1 /* mpn_divrem -- Divide natural numbers, producing both remainder and 2 quotient. 3 4 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc. 5 6 This file is part of the GNU MP Library. 7 8 The GNU MP Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 2.1 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MP Library; see the file COPYING.LIB. If not, write to 20 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, 21 MA 02111-1307, USA. */ 22 23 #include "gmp.h" 24 #include "gmp-impl.h" 25 #include "longlong.h" 26 27 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write 28 the NSIZE-DSIZE least significant quotient limbs at QP 29 and the DSIZE long remainder at NP. If QEXTRA_LIMBS is 30 non-zero, generate that many fraction bits and append them after the 31 other quotient limbs. 32 Return the most significant limb of the quotient, this is always 0 or 1. 33 34 Preconditions: 35 0. NSIZE >= DSIZE. 36 1. The most significant bit of the divisor must be set. 37 2. QP must either not overlap with the input operands at all, or 38 QP + DSIZE >= NP must hold true. (This means that it's 39 possible to put the quotient in the high part of NUM, right after the 40 remainder in NUM. 41 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. */ 42 43 mp_limb_t 44 #if __STDC__ 45 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs, 46 mp_ptr np, mp_size_t nsize, 47 mp_srcptr dp, mp_size_t dsize) 48 #else 49 mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize) 50 mp_ptr qp; 51 mp_size_t qextra_limbs; 52 mp_ptr np; 53 mp_size_t nsize; 54 mp_srcptr dp; 55 mp_size_t dsize; 56 #endif 57 { 58 mp_limb_t most_significant_q_limb = 0; 59 60 switch (dsize) 61 { 62 case 0: 63 /* We are asked to divide by zero, so go ahead and do it! (To make 64 the compiler not remove this statement, return the value.) */ 65 return 1 / dsize; 66 67 case 1: 68 { 69 mp_size_t i; 70 mp_limb_t n1; 71 mp_limb_t d; 72 73 d = dp[0]; 74 n1 = np[nsize - 1]; 75 76 if (n1 >= d) 77 { 78 n1 -= d; 79 most_significant_q_limb = 1; 80 } 81 82 qp += qextra_limbs; 83 for (i = nsize - 2; i >= 0; i--) 84 udiv_qrnnd (qp[i], n1, n1, np[i], d); 85 qp -= qextra_limbs; 86 87 for (i = qextra_limbs - 1; i >= 0; i--) 88 udiv_qrnnd (qp[i], n1, n1, 0, d); 89 90 np[0] = n1; 91 } 92 break; 93 94 case 2: 95 { 96 mp_size_t i; 97 mp_limb_t n1, n0, n2; 98 mp_limb_t d1, d0; 99 100 np += nsize - 2; 101 d1 = dp[1]; 102 d0 = dp[0]; 103 n1 = np[1]; 104 n0 = np[0]; 105 106 if (n1 >= d1 && (n1 > d1 || n0 >= d0)) 107 { 108 sub_ddmmss (n1, n0, n1, n0, d1, d0); 109 most_significant_q_limb = 1; 110 } 111 112 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) 113 { 114 mp_limb_t q; 115 mp_limb_t r; 116 117 if (i >= qextra_limbs) 118 np--; 119 else 120 np[0] = 0; 121 122 if (n1 == d1) 123 { 124 /* Q should be either 111..111 or 111..110. Need special 125 treatment of this rare case as normal division would 126 give overflow. */ 127 q = ~(mp_limb_t) 0; 128 129 r = n0 + d1; 130 if (r < d1) /* Carry in the addition? */ 131 { 132 add_ssaaaa (n1, n0, r - d0, np[0], 0, d0); 133 qp[i] = q; 134 continue; 135 } 136 n1 = d0 - (d0 != 0); 137 n0 = -d0; 138 } 139 else 140 { 141 udiv_qrnnd (q, r, n1, n0, d1); 142 umul_ppmm (n1, n0, d0, q); 143 } 144 145 n2 = np[0]; 146 q_test: 147 if (n1 > r || (n1 == r && n0 > n2)) 148 { 149 /* The estimated Q was too large. */ 150 q--; 151 152 sub_ddmmss (n1, n0, n1, n0, 0, d0); 153 r += d1; 154 if (r >= d1) /* If not carry, test Q again. */ 155 goto q_test; 156 } 157 158 qp[i] = q; 159 sub_ddmmss (n1, n0, r, n2, n1, n0); 160 } 161 np[1] = n1; 162 np[0] = n0; 163 } 164 break; 165 166 default: 167 { 168 mp_size_t i; 169 mp_limb_t dX, d1, n0; 170 171 np += nsize - dsize; 172 dX = dp[dsize - 1]; 173 d1 = dp[dsize - 2]; 174 n0 = np[dsize - 1]; 175 176 if (n0 >= dX) 177 { 178 if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0) 179 { 180 mpn_sub_n (np, np, dp, dsize); 181 n0 = np[dsize - 1]; 182 most_significant_q_limb = 1; 183 } 184 } 185 186 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) 187 { 188 mp_limb_t q; 189 mp_limb_t n1, n2; 190 mp_limb_t cy_limb; 191 192 if (i >= qextra_limbs) 193 { 194 np--; 195 n2 = np[dsize]; 196 } 197 else 198 { 199 n2 = np[dsize - 1]; 200 MPN_COPY_DECR (np + 1, np, dsize); 201 np[0] = 0; 202 } 203 204 if (n0 == dX) 205 /* This might over-estimate q, but it's probably not worth 206 the extra code here to find out. */ 207 q = ~(mp_limb_t) 0; 208 else 209 { 210 mp_limb_t r; 211 212 udiv_qrnnd (q, r, n0, np[dsize - 1], dX); 213 umul_ppmm (n1, n0, d1, q); 214 215 while (n1 > r || (n1 == r && n0 > np[dsize - 2])) 216 { 217 q--; 218 r += dX; 219 if (r < dX) /* I.e. "carry in previous addition?" */ 220 break; 221 n1 -= n0 < d1; 222 n0 -= d1; 223 } 224 } 225 226 /* Possible optimization: We already have (q * n0) and (1 * n1) 227 after the calculation of q. Taking advantage of that, we 228 could make this loop make two iterations less. */ 229 230 cy_limb = mpn_submul_1 (np, dp, dsize, q); 231 232 if (n2 != cy_limb) 233 { 234 mpn_add_n (np, np, dp, dsize); 235 q--; 236 } 237 238 qp[i] = q; 239 n0 = np[dsize - 1]; 240 } 241 } 242 } 243 244 return most_significant_q_limb; 245 } 246