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