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