--- /dev/null
+/* 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;
+}
projects / ghc-hetmet.git / blobdiff