divexact.c revision 1.1.1.2
1/* mpn_divexact(qp,np,nn,dp,dn,tp) -- Divide N = {np,nn} by D = {dp,dn} storing 2 the result in Q = {qp,nn-dn+1} expecting no remainder. Overlap allowed 3 between Q and N; all other overlap disallowed. 4 5 Contributed to the GNU project by Torbjorn Granlund. 6 7 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY 8 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 9 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. 10 11Copyright 2006, 2007, 2009 Free Software Foundation, Inc. 12 13This file is part of the GNU MP Library. 14 15The GNU MP Library is free software; you can redistribute it and/or modify 16it under the terms of the GNU Lesser General Public License as published by 17the Free Software Foundation; either version 3 of the License, or (at your 18option) any later version. 19 20The GNU MP Library is distributed in the hope that it will be useful, but 21WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 22or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 23License for more details. 24 25You should have received a copy of the GNU Lesser General Public License 26along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 27 28 29#include "gmp.h" 30#include "gmp-impl.h" 31#include "longlong.h" 32 33#if 1 34void 35mpn_divexact (mp_ptr qp, 36 mp_srcptr np, mp_size_t nn, 37 mp_srcptr dp, mp_size_t dn) 38{ 39 unsigned shift; 40 mp_size_t qn; 41 mp_ptr tp; 42 TMP_DECL; 43 44 ASSERT (dn > 0); 45 ASSERT (nn >= dn); 46 ASSERT (dp[dn-1] > 0); 47 48 while (dp[0] == 0) 49 { 50 ASSERT (np[0] == 0); 51 dp++; 52 np++; 53 dn--; 54 nn--; 55 } 56 57 if (dn == 1) 58 { 59 MPN_DIVREM_OR_DIVEXACT_1 (qp, np, nn, dp[0]); 60 return; 61 } 62 63 TMP_MARK; 64 65 qn = nn + 1 - dn; 66 count_trailing_zeros (shift, dp[0]); 67 68 if (shift > 0) 69 { 70 mp_ptr wp; 71 mp_size_t ss; 72 ss = (dn > qn) ? qn + 1 : dn; 73 74 tp = TMP_ALLOC_LIMBS (ss); 75 mpn_rshift (tp, dp, ss, shift); 76 dp = tp; 77 78 /* Since we have excluded dn == 1, we have nn > qn, and we need 79 to shift one limb beyond qn. */ 80 wp = TMP_ALLOC_LIMBS (qn + 1); 81 mpn_rshift (wp, np, qn + 1, shift); 82 np = wp; 83 } 84 85 if (dn > qn) 86 dn = qn; 87 88 tp = TMP_ALLOC_LIMBS (mpn_bdiv_q_itch (qn, dn)); 89 mpn_bdiv_q (qp, np, qn, dp, dn, tp); 90 TMP_FREE; 91} 92 93#else 94 95/* We use the Jebelean's bidirectional exact division algorithm. This is 96 somewhat naively implemented, with equal quotient parts done by 2-adic 97 division and truncating division. Since 2-adic division is faster, it 98 should be used for a larger chunk. 99 100 This code is horrendously ugly, in all sorts of ways. 101 102 * It was hacked without much care or thought, but with a testing program. 103 * It handles scratch space frivolously, and furthermore the itch function 104 is broken. 105 * Doesn't provide any measures to deal with mu_divappr_q's +3 error. We 106 have yet to provoke an error due to this, though. 107 * Algorithm selection leaves a lot to be desired. In particular, the choice 108 between DC and MU isn't a point, but we treat it like one. 109 * It makes the msb part 1 or 2 limbs larger than the lsb part, in spite of 110 that the latter is faster. We should at least reverse this, but perhaps 111 we should make the lsb part considerably larger. (How do we tune this?) 112*/ 113 114mp_size_t 115mpn_divexact_itch (mp_size_t nn, mp_size_t dn) 116{ 117 return nn + dn; /* FIXME this is not right */ 118} 119 120void 121mpn_divexact (mp_ptr qp, 122 mp_srcptr np, mp_size_t nn, 123 mp_srcptr dp, mp_size_t dn, 124 mp_ptr scratch) 125{ 126 mp_size_t qn; 127 mp_size_t nn0, qn0; 128 mp_size_t nn1, qn1; 129 mp_ptr tp; 130 mp_limb_t qml; 131 mp_limb_t qh; 132 int cnt; 133 mp_ptr xdp; 134 mp_limb_t di; 135 mp_limb_t cy; 136 gmp_pi1_t dinv; 137 TMP_DECL; 138 139 TMP_MARK; 140 141 qn = nn - dn + 1; 142 143 /* For small divisors, and small quotients, don't use Jebelean's algorithm. */ 144 if (dn < DIVEXACT_JEB_THRESHOLD || qn < DIVEXACT_JEB_THRESHOLD) 145 { 146 tp = scratch; 147 MPN_COPY (tp, np, qn); 148 binvert_limb (di, dp[0]); di = -di; 149 dn = MIN (dn, qn); 150 mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di); 151 TMP_FREE; 152 return; 153 } 154 155 qn0 = ((nn - dn) >> 1) + 1; /* low quotient size */ 156 157 /* If quotient is much larger than the divisor, the bidirectional algorithm 158 does not work as currently implemented. Fall back to plain bdiv. */ 159 if (qn0 > dn) 160 { 161 if (BELOW_THRESHOLD (dn, DC_BDIV_Q_THRESHOLD)) 162 { 163 tp = scratch; 164 MPN_COPY (tp, np, qn); 165 binvert_limb (di, dp[0]); di = -di; 166 dn = MIN (dn, qn); 167 mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di); 168 } 169 else if (BELOW_THRESHOLD (dn, MU_BDIV_Q_THRESHOLD)) 170 { 171 tp = scratch; 172 MPN_COPY (tp, np, qn); 173 binvert_limb (di, dp[0]); di = -di; 174 mpn_dcpi1_bdiv_q (qp, tp, qn, dp, dn, di); 175 } 176 else 177 { 178 mpn_mu_bdiv_q (qp, np, qn, dp, dn, scratch); 179 } 180 TMP_FREE; 181 return; 182 } 183 184 nn0 = qn0 + qn0; 185 186 nn1 = nn0 - 1 + ((nn-dn) & 1); 187 qn1 = qn0; 188 if (LIKELY (qn0 != dn)) 189 { 190 nn1 = nn1 + 1; 191 qn1 = qn1 + 1; 192 if (UNLIKELY (dp[dn - 1] == 1 && qn1 != dn)) 193 { 194 /* If the leading divisor limb == 1, i.e. has just one bit, we have 195 to include an extra limb in order to get the needed overlap. */ 196 /* FIXME: Now with the mu_divappr_q function, we should really need 197 more overlap. That indicates one of two things: (1) The test code 198 is not good. (2) We actually overlap too much by default. */ 199 nn1 = nn1 + 1; 200 qn1 = qn1 + 1; 201 } 202 } 203 204 tp = TMP_ALLOC_LIMBS (nn1 + 1); 205 206 count_leading_zeros (cnt, dp[dn - 1]); 207 208 /* Normalize divisor, store into tmp area. */ 209 if (cnt != 0) 210 { 211 xdp = TMP_ALLOC_LIMBS (qn1); 212 mpn_lshift (xdp, dp + dn - qn1, qn1, cnt); 213 } 214 else 215 { 216 xdp = (mp_ptr) dp + dn - qn1; 217 } 218 219 /* Shift dividend according to the divisor normalization. */ 220 /* FIXME: We compute too much here for XX_divappr_q, but these functions' 221 interfaces want a pointer to the imaginative least significant limb, not 222 to the least significant *used* limb. Of course, we could leave nn1-qn1 223 rubbish limbs in the low part, to save some time. */ 224 if (cnt != 0) 225 { 226 cy = mpn_lshift (tp, np + nn - nn1, nn1, cnt); 227 if (cy != 0) 228 { 229 tp[nn1] = cy; 230 nn1++; 231 } 232 } 233 else 234 { 235 /* FIXME: This copy is not needed for mpn_mu_divappr_q, except when the 236 mpn_sub_n right before is executed. */ 237 MPN_COPY (tp, np + nn - nn1, nn1); 238 } 239 240 invert_pi1 (dinv, xdp[qn1 - 1], xdp[qn1 - 2]); 241 if (BELOW_THRESHOLD (qn1, DC_DIVAPPR_Q_THRESHOLD)) 242 { 243 qp[qn0 - 1 + nn1 - qn1] = mpn_sbpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, dinv.inv32); 244 } 245 else if (BELOW_THRESHOLD (qn1, MU_DIVAPPR_Q_THRESHOLD)) 246 { 247 qp[qn0 - 1 + nn1 - qn1] = mpn_dcpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, &dinv); 248 } 249 else 250 { 251 /* FIXME: mpn_mu_divappr_q doesn't handle qh != 0. Work around it with a 252 conditional subtraction here. */ 253 qh = mpn_cmp (tp + nn1 - qn1, xdp, qn1) >= 0; 254 if (qh) 255 mpn_sub_n (tp + nn1 - qn1, tp + nn1 - qn1, xdp, qn1); 256 mpn_mu_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, scratch); 257 qp[qn0 - 1 + nn1 - qn1] = qh; 258 } 259 qml = qp[qn0 - 1]; 260 261 binvert_limb (di, dp[0]); di = -di; 262 263 if (BELOW_THRESHOLD (qn0, DC_BDIV_Q_THRESHOLD)) 264 { 265 MPN_COPY (tp, np, qn0); 266 mpn_sbpi1_bdiv_q (qp, tp, qn0, dp, qn0, di); 267 } 268 else if (BELOW_THRESHOLD (qn0, MU_BDIV_Q_THRESHOLD)) 269 { 270 MPN_COPY (tp, np, qn0); 271 mpn_dcpi1_bdiv_q (qp, tp, qn0, dp, qn0, di); 272 } 273 else 274 { 275 mpn_mu_bdiv_q (qp, np, qn0, dp, qn0, scratch); 276 } 277 278 if (qml < qp[qn0 - 1]) 279 mpn_decr_u (qp + qn0, 1); 280 281 TMP_FREE; 282} 283#endif 284