diff options
Diffstat (limited to 'gnu/lib/libgmp/mpn_div.c')
-rw-r--r-- | gnu/lib/libgmp/mpn_div.c | 321 |
1 files changed, 321 insertions, 0 deletions
diff --git a/gnu/lib/libgmp/mpn_div.c b/gnu/lib/libgmp/mpn_div.c new file mode 100644 index 0000000..8609206 --- /dev/null +++ b/gnu/lib/libgmp/mpn_div.c @@ -0,0 +1,321 @@ +/* mpn_div -- Divide natural numbers, producing both remainder and + quotient. + +Copyright (C) 1991 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2, or (at your option) +any later version. + +The GNU MP Library is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the GNU MP Library; see the file COPYING. If not, write to +the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +/* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write + the quotient at QUOT_PTR and the remainder at NUM_PTR. + + Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE) + or (NSIZE - DSIZE + 1). + + Argument constraints: + 1. The most significant bit of d must be set. + 2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage + area must be distinct from either input operands. + + The exact sizes of the quotient and remainder must be determined + by the caller, in spite of the return value. The return value just + informs the caller about if the highest digit is written or not, and + it may very well be 0. */ + +/* THIS WILL BE IMPROVED SOON. MORE COMMENTS AND FASTER CODE. */ + +mp_size +#ifdef __STDC__ +mpn_div (mp_ptr quot_ptr, + mp_ptr num_ptr, mp_size num_size, + mp_srcptr den_ptr, mp_size den_size) +#else +mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size) + mp_ptr quot_ptr; + mp_ptr num_ptr; + mp_size num_size; + mp_srcptr den_ptr; + mp_size den_size; +#endif +{ + mp_size q_is_long = 0; + + switch (den_size) + { + case 0: + /* We are asked to divide by zero, so go ahead and do it! + (To make the compiler not remove this statement, assign NUM_SIZE + and fall through.) */ + num_size = 1 / den_size; + + case 1: + { + mp_size i; + mp_limb n1, n0; + mp_limb d; + + d = den_ptr[0]; + i = num_size - 1; + n1 = num_ptr[i]; + i--; + + if (n1 >= d) + { + q_is_long = 1; + n1 = 0; + i++; + } + + for (; i >= 0; i--) + { + n0 = num_ptr[i]; + udiv_qrnnd (quot_ptr[i], n1, n1, n0, d); + } + + num_ptr[0] = n1; + } + break; + + case 2: + { + mp_size i; + mp_limb n0, n1, n2; + mp_limb d0, d1; + + num_ptr += num_size - 2; + d0 = den_ptr[1]; + d1 = den_ptr[0]; + n0 = num_ptr[1]; + n1 = num_ptr[0]; + + if (n0 >= d0) + { + q_is_long = 1; + n1 = n0; + n0 = 0; + num_ptr++; + num_size++; + } + + for (i = num_size - den_size - 1; i >= 0; i--) + { + mp_limb q; + mp_limb r; + + num_ptr--; + if (n0 == d0) + { + /* Q should be either 111..111 or 111..110. Need special + treatment of this rare case as normal division would + give overflow. */ + q = ~0; + + r = n1 + d0; + if (r < d0) /* Carry in the addition? */ + { + n2 = num_ptr[0]; + + add_ssaaaa (n0, n1, r - d1, n2, 0, d1); + quot_ptr[i] = q; + continue; + } + n0 = d1 - (d1 != 0); + n1 = -d1; + } + else + { + udiv_qrnnd (q, r, n0, n1, d0); + umul_ppmm (n0, n1, d1, q); + } + + n2 = num_ptr[0]; + q_test: + if (n0 > r || (n0 == r && n1 > n2)) + { + /* The estimated Q was too large. */ + q--; + + sub_ddmmss (n0, n1, n0, n1, 0, d1); + r += d0; + if (r >= d0) /* If not carry, test q again. */ + goto q_test; + } + + quot_ptr[i] = q; + sub_ddmmss (n0, n1, r, n2, n0, n1); + } + num_ptr[1] = n0; + num_ptr[0] = n1; + } + break; + + default: + { + mp_size i; + mp_limb d0 = den_ptr[den_size - 1]; + mp_limb d1 = den_ptr[den_size - 2]; + mp_limb n0 = num_ptr[num_size - 1]; + int ugly_hack_flag = 0; + + if (n0 >= d0) + { + + /* There's a problem with this case, which shows up later in the + code. q becomes 1 (and sometimes 0) the first time when + we've been here, and w_cy == 0 after the main do-loops below. + But c = num_ptr[j] reads rubbish outside the num_ptr vector! + Maybe I can solve this cleanly when I fix the early-end + optimization here in the default case. For now, I change the + add_back entering condition, to kludge. Leaving the stray + memref behind! + + HACK: Added ugly_hack_flag to make it work. */ + + q_is_long = 1; + n0 = 0; + num_size++; + ugly_hack_flag = 1; + } + + num_ptr += num_size; + den_ptr += den_size; + for (i = num_size - den_size - 1; i >= 0; i--) + { + mp_limb q; + mp_limb n1; + mp_limb w_cy; + mp_limb d, c; + mp_size j; + + num_ptr--; + if (n0 == d0) + /* This might over-estimate q, but it's probably not worth + the extra code here to find out. */ + q = ~0; + else + { + mp_limb r; + + udiv_qrnnd (q, r, n0, num_ptr[-1], d0); + umul_ppmm (n1, n0, d1, q); + + while (n1 > r || (n1 == r && n0 > num_ptr[-2])) + { + q--; + r += d0; + if (r < d0) /* I.e. "carry in previous addition?" */ + break; + n1 -= n0 < d1; + n0 -= d1; + } + } + + w_cy = 0; + j = -den_size; + do + { + d = den_ptr[j]; + c = num_ptr[j]; + umul_ppmm (n1, n0, d, q); + n0 += w_cy; + w_cy = (n0 < w_cy) + n1; + n0 = c - n0; + num_ptr[j] = n0; + if (n0 > c) + goto cy_loop; + ncy_loop: + j++; + } + while (j < 0); + + if (ugly_hack_flag) + { + c = 0; + ugly_hack_flag = 0; + } + else + c = num_ptr[j]; + if (c >= w_cy) + goto store_q; + goto add_back; + + do + { + d = den_ptr[j]; + c = num_ptr[j]; + umul_ppmm (n1, n0, d, q); + n0 += w_cy; + w_cy = (n0 < w_cy) + n1; + n0 = c - n0 - 1; + num_ptr[j] = n0; + if (n0 < c) + goto ncy_loop; + cy_loop: + j++; + } + while (j < 0); + + if (ugly_hack_flag) + { + c = 0; + ugly_hack_flag = 0; + } + else + c = num_ptr[j]; + w_cy++; + if (c >= w_cy) + goto store_q; + + add_back: + j = -den_size; + do + { + d = den_ptr[j]; + n0 = num_ptr[j] + d; + num_ptr[j] = n0; + if (n0 < d) + goto ab_cy_loop; + ab_ncy_loop: + j++; + } + while (j < 0); + abort (); /* We should always have a carry out! */ + + do + { + d = den_ptr[j]; + n0 = num_ptr[j] + d + 1; + num_ptr[j] = n0; + if (n0 > d) + goto ab_ncy_loop; + ab_cy_loop: + j++; + } + while (j < 0); + q--; + + store_q: + quot_ptr[i] = q; + } + } + } + + return q_is_long; +} |