Imported from gmp-1.900
authorroland <roland>
Fri, 17 Feb 1995 20:10:31 +0000 (20:10 +0000)
committerroland <roland>
Fri, 17 Feb 1995 20:10:31 +0000 (20:10 +0000)
15 files changed:
sysdeps/generic/add_1.c [new file with mode: 0644]
sysdeps/generic/add_n.c [new file with mode: 0644]
sysdeps/generic/addmul_1.c [new file with mode: 0644]
sysdeps/generic/cmp.c [new file with mode: 0644]
sysdeps/generic/divmod.c [new file with mode: 0644]
sysdeps/generic/divmod_1.c [new file with mode: 0644]
sysdeps/generic/gmp-mparam.h [new file with mode: 0644]
sysdeps/generic/lshift.c [new file with mode: 0644]
sysdeps/generic/mod_1.c [new file with mode: 0644]
sysdeps/generic/mul.c [new file with mode: 0644]
sysdeps/generic/mul_1.c [new file with mode: 0644]
sysdeps/generic/mul_n.c [new file with mode: 0644]
sysdeps/generic/rshift.c [new file with mode: 0644]
sysdeps/generic/sub_n.c [new file with mode: 0644]
sysdeps/generic/submul_1.c [new file with mode: 0644]

diff --git a/sysdeps/generic/add_1.c b/sysdeps/generic/add_1.c
new file mode 100644 (file)
index 0000000..7b1c697
--- /dev/null
@@ -0,0 +1,61 @@
+/* mpn_add_1 --
+
+Copyright (C) 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#define __mpn_add_1 __noname
+#include "gmp.h"
+#undef __mpn_add_1
+
+#include "gmp-impl.h"
+
+mp_limb
+__mpn_add_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     register mp_size_t s1_size;
+     register mp_limb s2_limb;
+{
+  register mp_limb x;
+
+  x = *s1_ptr++;
+  s2_limb = x + s2_limb;
+  *res_ptr++ = s2_limb;
+  if (s2_limb < x)
+    {
+      while (--s1_size != 0)
+       {
+         x = *s1_ptr++ + 1;
+         *res_ptr++ = x;
+         if (x != 0)
+           goto fin;
+       }
+
+      return 1;
+    }
+
+ fin:
+  if (res_ptr != s1_ptr)
+    {
+      mp_size_t i;
+      for (i = 0; i < s1_size - 1; i++)
+       res_ptr[i] = s1_ptr[i];
+    }
+
+  return 0;
+}
diff --git a/sysdeps/generic/add_n.c b/sysdeps/generic/add_n.c
new file mode 100644 (file)
index 0000000..6989ab0
--- /dev/null
@@ -0,0 +1,61 @@
+/* __mpn_add_n -- Add two limb vectors of equal, non-zero length.
+
+Copyright (C) 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mp_limb
+#if __STDC__
+__mpn_add_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+__mpn_add_n (res_ptr, s1_ptr, s2_ptr, size)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     register mp_srcptr s2_ptr;
+     mp_size_t size;
+#endif
+{
+  register mp_limb x, y, cy;
+  register mp_size_t j;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  s2_ptr -= j;
+  res_ptr -= j;
+
+  cy = 0;
+  do
+    {
+      y = s2_ptr[j];
+      x = s1_ptr[j];
+      y += cy;                 /* add previous carry to one addend */
+      cy = (y < cy);           /* get out carry from that addition */
+      y = x + y;               /* add other addend */
+      cy = (y < x) + cy;       /* get out carry from that add, combine */
+      res_ptr[j] = y;
+    }
+  while (++j != 0);
+
+  return cy;
+}
diff --git a/sysdeps/generic/addmul_1.c b/sysdeps/generic/addmul_1.c
new file mode 100644 (file)
index 0000000..fdf3541
--- /dev/null
@@ -0,0 +1,64 @@
+/* __mpn_addmul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
+   by S2_LIMB, add the S1_SIZE least significant limbs of the product to the
+   limb vector pointed to by RES_PTR.  Return the most significant limb of
+   the product, adjusted for carry-out from the addition.
+
+Copyright (C) 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+mp_limb
+__mpn_addmul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb s2_limb;
+{
+  register mp_limb cy_limb;
+  register mp_size_t j;
+  register mp_limb prod_high, prod_low;
+  register mp_limb x;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  res_ptr -= j;
+  s1_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      x = res_ptr[j];
+      prod_low = x + prod_low;
+      cy_limb += (prod_low < x);
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}
diff --git a/sysdeps/generic/cmp.c b/sysdeps/generic/cmp.c
new file mode 100644 (file)
index 0000000..144c885
--- /dev/null
@@ -0,0 +1,55 @@
+/* __mpn_cmp -- Compare two low-level natural-number integers.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Compare OP1_PTR/OP1_SIZE with OP2_PTR/OP2_SIZE.
+   There are no restrictions on the relative sizes of
+   the two arguments.
+   Return 1 if OP1 > OP2, 0 if they are equal, and -1 if OP1 < OP2.  */
+
+int
+#if __STDC__
+__mpn_cmp (mp_srcptr op1_ptr, mp_srcptr op2_ptr, mp_size_t size)
+#else
+__mpn_cmp (op1_ptr, op2_ptr, size)
+     mp_srcptr op1_ptr;
+     mp_srcptr op2_ptr;
+     mp_size_t size;
+#endif
+{
+  mp_size_t i;
+  mp_limb op1_word, op2_word;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      op1_word = op1_ptr[i];
+      op2_word = op2_ptr[i];
+      if (op1_word != op2_word)
+       goto diff;
+    }
+  return 0;
+ diff:
+  /* This can *not* be simplified to
+       op2_word - op2_word
+     since that expression might give signed overflow.  */
+  return (op1_word > op2_word) ? 1 : -1;
+}
diff --git a/sysdeps/generic/divmod.c b/sysdeps/generic/divmod.c
new file mode 100644 (file)
index 0000000..76b9bca
--- /dev/null
@@ -0,0 +1,234 @@
+/* __mpn_divmod -- Divide natural numbers, producing both remainder and
+   quotient.
+
+Copyright (C) 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 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 NUM_SIZE-DEN_SIZE least significant quotient limbs at QUOT_PTR
+   and the DEN_SIZE long remainder at NUM_PTR.
+   Return the most significant limb of the quotient, this is always 0 or 1.
+
+   Argument constraints:
+   1. The most significant bit of the divisor must be set.
+   2. QUOT_PTR must either not overlap with the input operands at all, or
+      QUOT_PTR + DEN_SIZE >= NUM_PTR must hold true.  (This means that it's
+      possible to put the quotient in the high part of NUM, right after the
+      remainder in NUM.  */
+
+mp_limb
+#if __STDC__
+__mpn_divmod (mp_ptr quot_ptr,
+             mp_ptr num_ptr, mp_size_t num_size,
+             mp_srcptr den_ptr, mp_size_t den_size)
+#else
+__mpn_divmod (quot_ptr, num_ptr, num_size, den_ptr, den_size)
+     mp_ptr quot_ptr;
+     mp_ptr num_ptr;
+     mp_size_t num_size;
+     mp_srcptr den_ptr;
+     mp_size_t den_size;
+#endif
+{
+  mp_limb most_significant_q_limb = 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, return the value.)  */
+      return 1 / den_size;
+
+    case 1:
+      {
+       mp_size_t i;
+       mp_limb n1, n0;
+       mp_limb d;
+
+       d = den_ptr[0];
+       n1 = num_ptr[num_size - 1];
+
+       if (n1 >= d)
+         {
+           most_significant_q_limb = 1;
+           n1 -= d;
+         }
+
+       for (i = num_size - 2; 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_t i;
+       mp_limb n1, n0, n2;
+       mp_limb d1, d0;
+
+       num_ptr += num_size - 2;
+       d1 = den_ptr[1];
+       d0 = den_ptr[0];
+       n1 = num_ptr[1];
+       n0 = num_ptr[0];
+
+       if (n1 >= d1 && (n1 > d1 || n0 >= d0))
+         {
+           most_significant_q_limb = 1;
+           sub_ddmmss (n1, n0, n1, n0, d1, d0);
+         }
+
+       for (i = num_size - den_size - 1; i >= 0; i--)
+         {
+           mp_limb q;
+           mp_limb r;
+
+           num_ptr--;
+           if (n1 == d1)
+             {
+               /* Q should be either 111..111 or 111..110.  Need special
+                  treatment of this rare case as normal division would
+                  give overflow.  */
+               q = ~(mp_limb) 0;
+
+               r = n0 + d1;
+               if (r < d1)     /* Carry in the addition? */
+                 {
+                   add_ssaaaa (n1, n0, r - d0, num_ptr[0], 0, d0);
+                   quot_ptr[i] = q;
+                   continue;
+                 }
+               n1 = d0 - (d0 != 0);
+               n0 = -d0;
+             }
+           else
+             {
+               udiv_qrnnd (q, r, n1, n0, d1);
+               umul_ppmm (n1, n0, d0, q);
+             }
+
+           n2 = num_ptr[0];
+         q_test:
+           if (n1 > r || (n1 == r && n0 > n2))
+             {
+               /* The estimated Q was too large.  */
+               q--;
+
+               sub_ddmmss (n1, n0, n1, n0, 0, d0);
+               r += d1;
+               if (r >= d1)    /* If not carry, test Q again.  */
+                 goto q_test;
+             }
+
+           quot_ptr[i] = q;
+           sub_ddmmss (n1, n0, r, n2, n1, n0);
+         }
+       num_ptr[1] = n1;
+       num_ptr[0] = n0;
+      }
+      break;
+
+    default:
+      {
+       mp_size_t i;
+       mp_limb dX, d1, n0;
+
+       num_ptr += num_size;
+       den_ptr += den_size;
+       dX = den_ptr[-1];
+       d1 = den_ptr[-2];
+       n0 = num_ptr[-1];
+
+       if (n0 >= dX)
+         {
+           if (n0 > dX
+               || __mpn_cmp (num_ptr - den_size, den_ptr - den_size,
+                             den_size - 1) >= 0)
+             {
+               __mpn_sub_n (num_ptr - den_size,
+                            num_ptr - den_size, den_ptr - den_size,
+                            den_size);
+               most_significant_q_limb = 1;
+             }
+
+           n0 = num_ptr[-1];
+         }
+
+       for (i = num_size - den_size - 1; i >= 0; i--)
+         {
+           mp_limb q;
+           mp_limb n1;
+           mp_limb cy_limb;
+
+           num_ptr--;
+           if (n0 == dX)
+             /* This might over-estimate q, but it's probably not worth
+                the extra code here to find out.  */
+             q = ~(mp_limb) 0;
+           else
+             {
+               mp_limb r;
+
+               udiv_qrnnd (q, r, n0, num_ptr[-1], dX);
+               umul_ppmm (n1, n0, d1, q);
+
+               while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
+                 {
+                   q--;
+                   r += dX;
+                   if (r < dX) /* I.e. "carry in previous addition?"  */
+                     break;
+                   n1 -= n0 < d1;
+                   n0 -= d1;
+                 }
+             }
+
+           /* Possible optimization: We already have (q * n0) and (1 * n1)
+              after the calculation of q.  Taking advantage of that, we
+              could make this loop make two iterations less.  */
+
+           cy_limb = __mpn_submul_1 (num_ptr - den_size,
+                                     den_ptr - den_size, den_size, q);
+
+           if (num_ptr[0] != cy_limb)
+             {
+               mp_limb cy;
+               cy = __mpn_add_n (num_ptr - den_size,
+                                 num_ptr - den_size,
+                                 den_ptr - den_size, den_size);
+               if (cy == 0)
+                 abort ();
+               q--;
+             }
+
+           quot_ptr[i] = q;
+           n0 = num_ptr[-1];
+         }
+      }
+    }
+
+  return most_significant_q_limb;
+}
diff --git a/sysdeps/generic/divmod_1.c b/sysdeps/generic/divmod_1.c
new file mode 100644 (file)
index 0000000..d156eeb
--- /dev/null
@@ -0,0 +1,209 @@
+/* __mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
+   Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
+   Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
+   Return the single-limb remainder.
+   There are no constraints on the value of the divisor.
+
+   QUOT_PTR and DIVIDEND_PTR might point to the same limb.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef UMUL_TIME
+#define UMUL_TIME 1
+#endif
+
+#ifndef UDIV_TIME
+#define UDIV_TIME UMUL_TIME
+#endif
+
+/* FIXME: We should be using invert_limb (or invert_normalized_limb)
+   here (not udiv_qrnnd).  */
+
+mp_limb
+#if __STDC__
+__mpn_divmod_1 (mp_ptr quot_ptr,
+             mp_srcptr dividend_ptr, mp_size_t dividend_size,
+             mp_limb divisor_limb)
+#else
+__mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
+     mp_ptr quot_ptr;
+     mp_srcptr dividend_ptr;
+     mp_size_t dividend_size;
+     mp_limb divisor_limb;
+#endif
+{
+  mp_size_t i;
+  mp_limb n1, n0, r;
+  int dummy;
+
+  /* ??? Should this be handled at all?  Rely on callers?  */
+  if (dividend_size == 0)
+    return 0;
+
+  /* If multiplication is much faster than division, and the
+     dividend is large, pre-invert the divisor, and use
+     only multiplications in the inner loop.  */
+
+  /* This test should be read:
+       Does it ever help to use udiv_qrnnd_preinv?
+        && Does what we save compensate for the inversion overhead?  */
+  if (UDIV_TIME > (2 * UMUL_TIME + 6)
+      && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
+    {
+      int normalization_steps;
+
+      count_leading_zeros (normalization_steps, divisor_limb);
+      if (normalization_steps != 0)
+       {
+         mp_limb divisor_limb_inverted;
+
+         divisor_limb <<= normalization_steps;
+
+         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+            result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+            most significant bit (with weight 2**N) implicit.  */
+
+#if 0 /* This can't happen when normalization_steps != 0 */
+         /* Special case for DIVISOR_LIMB == 100...000.  */
+         if (divisor_limb << 1 == 0)
+           divisor_limb_inverted = ~(mp_limb) 0;
+         else
+#endif
+         udiv_qrnnd (divisor_limb_inverted, dummy,
+                     -divisor_limb, 0, divisor_limb);
+
+         n1 = dividend_ptr[dividend_size - 1];
+         r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+         /* Possible optimization:
+            if (r == 0
+            && divisor_limb > ((n1 << normalization_steps)
+                            | (dividend_ptr[dividend_size - 2] >> ...)))
+            ...one division less... */
+
+         for (i = dividend_size - 2; i >= 0; i--)
+           {
+             n0 = dividend_ptr[i];
+             udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
+                                ((n1 << normalization_steps)
+                                 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+                                divisor_limb, divisor_limb_inverted);
+             n1 = n0;
+           }
+         udiv_qrnnd_preinv (quot_ptr[0], r, r,
+                            n1 << normalization_steps,
+                            divisor_limb, divisor_limb_inverted);
+         return r >> normalization_steps;
+       }
+      else
+       {
+         mp_limb divisor_limb_inverted;
+
+         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+            result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+            most significant bit (with weight 2**N) implicit.  */
+
+         /* Special case for DIVISOR_LIMB == 100...000.  */
+         if (divisor_limb << 1 == 0)
+           divisor_limb_inverted = ~(mp_limb) 0;
+         else
+           udiv_qrnnd (divisor_limb_inverted, dummy,
+                       -divisor_limb, 0, divisor_limb);
+
+         i = dividend_size - 1;
+         r = dividend_ptr[i];
+
+         if (r >= divisor_limb)
+           r = 0;
+         else
+           {
+             quot_ptr[i] = 0;
+             i--;
+           }
+
+         for (; i >= 0; i--)
+           {
+             n0 = dividend_ptr[i];
+             udiv_qrnnd_preinv (quot_ptr[i], r, r,
+                                n0, divisor_limb, divisor_limb_inverted);
+           }
+         return r;
+       }
+    }
+  else
+    {
+      if (UDIV_NEEDS_NORMALIZATION)
+       {
+         int normalization_steps;
+
+         count_leading_zeros (normalization_steps, divisor_limb);
+         if (normalization_steps != 0)
+           {
+             divisor_limb <<= normalization_steps;
+
+             n1 = dividend_ptr[dividend_size - 1];
+             r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+             /* Possible optimization:
+                if (r == 0
+                && divisor_limb > ((n1 << normalization_steps)
+                                | (dividend_ptr[dividend_size - 2] >> ...)))
+                ...one division less... */
+
+             for (i = dividend_size - 2; i >= 0; i--)
+               {
+                 n0 = dividend_ptr[i];
+                 udiv_qrnnd (quot_ptr[i + 1], r, r,
+                             ((n1 << normalization_steps)
+                              | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+                             divisor_limb);
+                 n1 = n0;
+               }
+             udiv_qrnnd (quot_ptr[0], r, r,
+                         n1 << normalization_steps,
+                         divisor_limb);
+             return r >> normalization_steps;
+           }
+       }
+      /* No normalization needed, either because udiv_qrnnd doesn't require
+        it, or because DIVISOR_LIMB is already normalized.  */
+
+      i = dividend_size - 1;
+      r = dividend_ptr[i];
+
+      if (r >= divisor_limb)
+       r = 0;
+      else
+       {
+         quot_ptr[i] = 0;
+         i--;
+       }
+
+      for (; i >= 0; i--)
+       {
+         n0 = dividend_ptr[i];
+         udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
+       }
+      return r;
+    }
+}
diff --git a/sysdeps/generic/gmp-mparam.h b/sysdeps/generic/gmp-mparam.h
new file mode 100644 (file)
index 0000000..4286ebf
--- /dev/null
@@ -0,0 +1,26 @@
+/* gmp-mparam.h -- Compiler/machine parameter header file.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#define BITS_PER_MP_LIMB 32
+#define BYTES_PER_MP_LIMB 4
+#define BITS_PER_LONGINT 32
+#define BITS_PER_INT 32
+#define BITS_PER_SHORTINT 16
+#define BITS_PER_CHAR 8
diff --git a/sysdeps/generic/lshift.c b/sysdeps/generic/lshift.c
new file mode 100644 (file)
index 0000000..1ba0903
--- /dev/null
@@ -0,0 +1,86 @@
+/* __mpn_lshift -- Shift left low level.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the left
+   and store the USIZE least significant digits of the result at WP.
+   Return the bits shifted out from the most significant digit.
+
+   Argument constraints:
+   1. 0 < CNT < BITS_PER_MP_LIMB
+   2. If the result is to be written over the input, WP must be >= UP.
+*/
+
+mp_limb
+#if __STDC__
+__mpn_lshift (register mp_ptr wp,
+           register mp_srcptr up, mp_size_t usize,
+           register unsigned int cnt)
+#else
+__mpn_lshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+#endif
+{
+  register mp_limb high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb retval;
+
+#ifdef DEBUG
+  if (usize == 0 || cnt == 0)
+    abort ();
+#endif
+
+  sh_1 = cnt;
+#if 0
+  if (sh_1 == 0)
+    {
+      if (wp != up)
+       {
+         /* Copy from high end to low end, to allow specified input/output
+            overlapping.  */
+         for (i = usize - 1; i >= 0; i--)
+           wp[i] = up[i];
+       }
+      return 0;
+    }
+#endif
+
+  wp += 1;
+  sh_2 = BITS_PER_MP_LIMB - sh_1;
+  i = usize - 1;
+  low_limb = up[i];
+  retval = low_limb >> sh_2;
+  high_limb = low_limb;
+  while (--i >= 0)
+    {
+      low_limb = up[i];
+      wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
+      high_limb = low_limb;
+    }
+  wp[i] = high_limb << sh_1;
+
+  return retval;
+}
diff --git a/sysdeps/generic/mod_1.c b/sysdeps/generic/mod_1.c
new file mode 100644 (file)
index 0000000..ae4ed09
--- /dev/null
@@ -0,0 +1,198 @@
+/* __mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
+   Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
+   Return the single-limb remainder.
+   There are no constraints on the value of the divisor.
+
+   QUOT_PTR and DIVIDEND_PTR might point to the same limb.
+
+Copyright (C) 1991, 1993, 1994, 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef UMUL_TIME
+#define UMUL_TIME 1
+#endif
+
+#ifndef UDIV_TIME
+#define UDIV_TIME UMUL_TIME
+#endif
+
+/* FIXME: We should be using invert_limb (or invert_normalized_limb)
+   here (not udiv_qrnnd).  */
+
+mp_limb
+#if __STDC__
+__mpn_mod_1 (mp_srcptr dividend_ptr, mp_size_t dividend_size,
+            mp_limb divisor_limb)
+#else
+__mpn_mod_1 (dividend_ptr, dividend_size, divisor_limb)
+     mp_srcptr dividend_ptr;
+     mp_size_t dividend_size;
+     mp_limb divisor_limb;
+#endif
+{
+  mp_size_t i;
+  mp_limb n1, n0, r;
+  int dummy;
+
+  /* Botch: Should this be handled at all?  Rely on callers?  */
+  if (dividend_size == 0)
+    return 0;
+
+  /* If multiplication is much faster than division, and the
+     dividend is large, pre-invert the divisor, and use
+     only multiplications in the inner loop.  */
+
+  /* This test should be read:
+       Does it ever help to use udiv_qrnnd_preinv?
+        && Does what we save compensate for the inversion overhead?  */
+  if (UDIV_TIME > (2 * UMUL_TIME + 6)
+      && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
+    {
+      int normalization_steps;
+
+      count_leading_zeros (normalization_steps, divisor_limb);
+      if (normalization_steps != 0)
+       {
+         mp_limb divisor_limb_inverted;
+
+         divisor_limb <<= normalization_steps;
+
+         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+            result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+            most significant bit (with weight 2**N) implicit.  */
+
+         /* Special case for DIVISOR_LIMB == 100...000.  */
+         if (divisor_limb << 1 == 0)
+           divisor_limb_inverted = ~(mp_limb) 0;
+         else
+           udiv_qrnnd (divisor_limb_inverted, dummy,
+                       -divisor_limb, 0, divisor_limb);
+
+         n1 = dividend_ptr[dividend_size - 1];
+         r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+         /* Possible optimization:
+            if (r == 0
+            && divisor_limb > ((n1 << normalization_steps)
+                            | (dividend_ptr[dividend_size - 2] >> ...)))
+            ...one division less... */
+
+         for (i = dividend_size - 2; i >= 0; i--)
+           {
+             n0 = dividend_ptr[i];
+             udiv_qrnnd_preinv (dummy, r, r,
+                                ((n1 << normalization_steps)
+                                 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+                                divisor_limb, divisor_limb_inverted);
+             n1 = n0;
+           }
+         udiv_qrnnd_preinv (dummy, r, r,
+                            n1 << normalization_steps,
+                            divisor_limb, divisor_limb_inverted);
+         return r >> normalization_steps;
+       }
+      else
+       {
+         mp_limb divisor_limb_inverted;
+
+         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+            result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+            most significant bit (with weight 2**N) implicit.  */
+
+         /* Special case for DIVISOR_LIMB == 100...000.  */
+         if (divisor_limb << 1 == 0)
+           divisor_limb_inverted = ~(mp_limb) 0;
+         else
+           udiv_qrnnd (divisor_limb_inverted, dummy,
+                       -divisor_limb, 0, divisor_limb);
+
+         i = dividend_size - 1;
+         r = dividend_ptr[i];
+
+         if (r >= divisor_limb)
+           r = 0;
+         else
+           i--;
+
+         for (; i >= 0; i--)
+           {
+             n0 = dividend_ptr[i];
+             udiv_qrnnd_preinv (dummy, r, r,
+                                n0, divisor_limb, divisor_limb_inverted);
+           }
+         return r;
+       }
+    }
+  else
+    {
+      if (UDIV_NEEDS_NORMALIZATION)
+       {
+         int normalization_steps;
+
+         count_leading_zeros (normalization_steps, divisor_limb);
+         if (normalization_steps != 0)
+           {
+             divisor_limb <<= normalization_steps;
+
+             n1 = dividend_ptr[dividend_size - 1];
+             r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+             /* Possible optimization:
+                if (r == 0
+                && divisor_limb > ((n1 << normalization_steps)
+                                | (dividend_ptr[dividend_size - 2] >> ...)))
+                ...one division less... */
+
+             for (i = dividend_size - 2; i >= 0; i--)
+               {
+                 n0 = dividend_ptr[i];
+                 udiv_qrnnd (dummy, r, r,
+                             ((n1 << normalization_steps)
+                              | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+                             divisor_limb);
+                 n1 = n0;
+               }
+             udiv_qrnnd (dummy, r, r,
+                         n1 << normalization_steps,
+                         divisor_limb);
+             return r >> normalization_steps;
+           }
+       }
+      /* No normalization needed, either because udiv_qrnnd doesn't require
+        it, or because DIVISOR_LIMB is already normalized.  */
+
+      i = dividend_size - 1;
+      r = dividend_ptr[i];
+
+      if (r >= divisor_limb)
+       r = 0;
+      else
+       i--;
+
+      for (; i >= 0; i--)
+       {
+         n0 = dividend_ptr[i];
+         udiv_qrnnd (dummy, r, r, n0, divisor_limb);
+       }
+      return r;
+    }
+}
diff --git a/sysdeps/generic/mul.c b/sysdeps/generic/mul.c
new file mode 100644 (file)
index 0000000..cd2acb5
--- /dev/null
@@ -0,0 +1,147 @@
+/* __mpn_mul -- Multiply two natural numbers.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
+   and v (pointed to by VP, with VSIZE limbs), and store the result at
+   PRODP.  USIZE + VSIZE limbs are always stored, but if the input
+   operands are normalized.  Return the most significant limb of the
+   result.
+
+   NOTE: The space pointed to by PRODP is overwritten before finished
+   with U and V, so overlap is an error.
+
+   Argument constraints:
+   1. USIZE >= VSIZE.
+   2. PRODP != UP and PRODP != VP, i.e. the destination
+      must be distinct from the multiplier and the multiplicand.  */
+
+/* If KARATSUBA_THRESHOLD is not already defined, define it to a
+   value which is good on most machines.  */
+#ifndef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 32
+#endif
+
+mp_limb
+#if __STDC__
+__mpn_mul (mp_ptr prodp,
+         mp_srcptr up, mp_size_t usize,
+         mp_srcptr vp, mp_size_t vsize)
+#else
+__mpn_mul (prodp, up, usize, vp, vsize)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_size_t usize;
+     mp_srcptr vp;
+     mp_size_t vsize;
+#endif
+{
+  mp_ptr prod_endp = prodp + usize + vsize - 1;
+  mp_limb cy;
+  mp_ptr tspace;
+
+  if (vsize < KARATSUBA_THRESHOLD)
+    {
+      /* Handle simple cases with traditional multiplication.
+
+        This is the most critical code of the entire function.  All
+        multiplies rely on this, both small and huge.  Small ones arrive
+        here immediately.  Huge ones arrive here as this is the base case
+        for Karatsuba's recursive algorithm below.  */
+      mp_size_t i;
+      mp_limb cy_limb;
+      mp_limb v_limb;
+
+      if (vsize == 0)
+       return 0;
+
+      /* Multiply by the first limb in V separately, as the result can be
+        stored (not added) to PROD.  We also avoid a loop for zeroing.  */
+      v_limb = vp[0];
+      if (v_limb <= 1)
+       {
+         if (v_limb == 1)
+           MPN_COPY (prodp, up, usize);
+         else
+           MPN_ZERO (prodp, usize);
+         cy_limb = 0;
+       }
+      else
+       cy_limb = __mpn_mul_1 (prodp, up, usize, v_limb);
+
+      prodp[usize] = cy_limb;
+      prodp++;
+
+      /* For each iteration in the outer loop, multiply one limb from
+        U with one limb from V, and add it to PROD.  */
+      for (i = 1; i < vsize; i++)
+       {
+         v_limb = vp[i];
+         if (v_limb <= 1)
+           {
+             cy_limb = 0;
+             if (v_limb == 1)
+               cy_limb = __mpn_add_n (prodp, prodp, up, usize);
+           }
+         else
+           cy_limb = __mpn_addmul_1 (prodp, up, usize, v_limb);
+
+         prodp[usize] = cy_limb;
+         prodp++;
+       }
+      return cy_limb;
+    }
+
+  tspace = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
+  MPN_MUL_N_RECURSE (prodp, up, vp, vsize, tspace);
+
+  prodp += vsize;
+  up += vsize;
+  usize -= vsize;
+  if (usize >= vsize)
+    {
+      mp_ptr tp = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
+      do
+       {
+         MPN_MUL_N_RECURSE (tp, up, vp, vsize, tspace);
+         cy = __mpn_add_n (prodp, prodp, tp, vsize);
+         __mpn_add_1 (prodp + vsize, tp + vsize, vsize, cy);
+         prodp += vsize;
+         up += vsize;
+         usize -= vsize;
+       }
+      while (usize >= vsize);
+    }
+
+  /* True: usize < vsize.  */
+
+  /* Make life simple: Recurse.  */
+
+  if (usize != 0)
+    {
+      __mpn_mul (tspace, vp, vsize, up, usize);
+      cy = __mpn_add_n (prodp, prodp, tspace, vsize);
+      __mpn_add_1 (prodp + vsize, tspace + vsize, usize, cy);
+    }
+
+  return *prod_endp;
+}
diff --git a/sysdeps/generic/mul_1.c b/sysdeps/generic/mul_1.c
new file mode 100644 (file)
index 0000000..37dbc33
--- /dev/null
@@ -0,0 +1,58 @@
+/* __mpn_mul_1 -- Multiply a limb vector with a single limb and
+   store the product in a second limb vector.
+
+Copyright (C) 1991, 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+mp_limb
+__mpn_mul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb s2_limb;
+{
+  register mp_limb cy_limb;
+  register mp_size_t j;
+  register mp_limb prod_high, prod_low;
+
+  /* The loop counter and index J goes from -S1_SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  res_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}
diff --git a/sysdeps/generic/mul_n.c b/sysdeps/generic/mul_n.c
new file mode 100644 (file)
index 0000000..7900988
--- /dev/null
@@ -0,0 +1,420 @@
+/* __mpn_mul_n -- Multiply two natural numbers of length n.
+
+Copyright (C) 1991, 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
+   both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
+   always stored.  Return the most significant limb.
+
+   Argument constraints:
+   1. PRODP != UP and PRODP != VP, i.e. the destination
+      must be distinct from the multiplier and the multiplicand.  */
+
+/* If KARATSUBA_THRESHOLD is not already defined, define it to a
+   value which is good on most machines.  */
+#ifndef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 32
+#endif
+
+/* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
+#if KARATSUBA_THRESHOLD < 2
+#undef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 2
+#endif
+
+void
+#if __STDC__
+____mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t, mp_ptr);
+#else
+____mpn_mul_n ();
+#endif
+
+/* Handle simple cases with traditional multiplication.
+
+   This is the most critical code of multiplication.  All multiplies rely
+   on this, both small and huge.  Small ones arrive here immediately.  Huge
+   ones arrive here as this is the base case for Karatsuba's recursive
+   algorithm below.  */
+
+void
+#if __STDC__
+____mpn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
+#else
+____mpn_mul_n_basecase (prodp, up, vp, size)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_srcptr vp;
+     mp_size_t size;
+#endif
+{
+  mp_size_t i;
+  mp_limb cy_limb;
+  mp_limb v_limb;
+
+  /* Multiply by the first limb in V separately, as the result can be
+     stored (not added) to PROD.  We also avoid a loop for zeroing.  */
+  v_limb = vp[0];
+  if (v_limb <= 1)
+    {
+      if (v_limb == 1)
+       MPN_COPY (prodp, up, size);
+      else
+       MPN_ZERO (prodp, size);
+      cy_limb = 0;
+    }
+  else
+    cy_limb = __mpn_mul_1 (prodp, up, size, v_limb);
+
+  prodp[size] = cy_limb;
+  prodp++;
+
+  /* For each iteration in the outer loop, multiply one limb from
+     U with one limb from V, and add it to PROD.  */
+  for (i = 1; i < size; i++)
+    {
+      v_limb = vp[i];
+      if (v_limb <= 1)
+       {
+         cy_limb = 0;
+         if (v_limb == 1)
+           cy_limb = __mpn_add_n (prodp, prodp, up, size);
+       }
+      else
+       cy_limb = __mpn_addmul_1 (prodp, up, size, v_limb);
+
+      prodp[size] = cy_limb;
+      prodp++;
+    }
+}
+
+void
+#if __STDC__
+____mpn_mul_n (mp_ptr prodp,
+            mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
+#else
+____mpn_mul_n (prodp, up, vp, size, tspace)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_srcptr vp;
+     mp_size_t size;
+     mp_ptr tspace;
+#endif
+{
+  if ((size & 1) != 0)
+    {
+      /* The size is odd, the code code below doesn't handle that.
+        Multiply the least significant (size - 1) limbs with a recursive
+        call, and handle the most significant limb of S1 and S2
+        separately.  */
+      /* A slightly faster way to do this would be to make the Karatsuba
+        code below behave as if the size were even, and let it check for
+        odd size in the end.  I.e., in essence move this code to the end.
+        Doing so would save us a recursive call, and potentially make the
+        stack grow a lot less.  */
+
+      mp_size_t esize = size - 1;      /* even size */
+      mp_limb cy_limb;
+
+      MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
+      cy_limb = __mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
+      prodp[esize + esize] = cy_limb;
+      cy_limb = __mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
+
+      prodp[esize + size] = cy_limb;
+    }
+  else
+    {
+      /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
+
+        Split U in two pieces, U1 and U0, such that
+        U = U0 + U1*(B**n),
+        and V in V1 and V0, such that
+        V = V0 + V1*(B**n).
+
+        UV is then computed recursively using the identity
+
+               2n   n          n                     n
+        UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
+                       1 1        1  0   0  1              0 0
+
+        Where B = 2**BITS_PER_MP_LIMB.  */
+
+      mp_size_t hsize = size >> 1;
+      mp_limb cy;
+      int negflg;
+
+      /*** Product H.   ________________  ________________
+                       |_____U1 x V1____||____U0 x V0_____|  */
+      /* Put result in upper part of PROD and pass low part of TSPACE
+        as new TSPACE.  */
+      MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
+
+      /*** Product M.   ________________
+                       |_(U1-U0)(V0-V1)_|  */
+      if (__mpn_cmp (up + hsize, up, hsize) >= 0)
+       {
+         __mpn_sub_n (prodp, up + hsize, up, hsize);
+         negflg = 0;
+       }
+      else
+       {
+         __mpn_sub_n (prodp, up, up + hsize, hsize);
+         negflg = 1;
+       }
+      if (__mpn_cmp (vp + hsize, vp, hsize) >= 0)
+       {
+         __mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
+         negflg ^= 1;
+       }
+      else
+       {
+         __mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
+         /* No change of NEGFLG.  */
+       }
+      /* Read temporary operands from low part of PROD.
+        Put result in low part of TSPACE using upper part of TSPACE
+        as new TSPACE.  */
+      MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
+
+      /*** Add/copy product H.  */
+      MPN_COPY (prodp + hsize, prodp + size, hsize);
+      cy = __mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
+
+      /*** Add product M (if NEGFLG M is a negative number).  */
+      if (negflg)
+       cy -= __mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
+      else
+       cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
+
+      /*** Product L.   ________________  ________________
+                       |________________||____U0 x V0_____|  */
+      /* Read temporary operands from low part of PROD.
+        Put result in low part of TSPACE using upper part of TSPACE
+        as new TSPACE.  */
+      MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
+
+      /*** Add/copy Product L (twice).  */
+
+      cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
+      if (cy)
+       {
+         if (cy > 0)
+           __mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
+         else
+           {
+             __mpn_sub_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
+             abort ();
+           }
+       }
+
+      MPN_COPY (prodp, tspace, hsize);
+      cy = __mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
+      if (cy)
+       __mpn_add_1 (prodp + size, prodp + size, size, 1);
+    }
+}
+
+void
+#if __STDC__
+____mpn_sqr_n_basecase (mp_ptr prodp, mp_srcptr up, mp_size_t size)
+#else
+____mpn_sqr_n_basecase (prodp, up, size)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_size_t size;
+#endif
+{
+  mp_size_t i;
+  mp_limb cy_limb;
+  mp_limb v_limb;
+
+  /* Multiply by the first limb in V separately, as the result can be
+     stored (not added) to PROD.  We also avoid a loop for zeroing.  */
+  v_limb = up[0];
+  if (v_limb <= 1)
+    {
+      if (v_limb == 1)
+       MPN_COPY (prodp, up, size);
+      else
+       MPN_ZERO (prodp, size);
+      cy_limb = 0;
+    }
+  else
+    cy_limb = __mpn_mul_1 (prodp, up, size, v_limb);
+
+  prodp[size] = cy_limb;
+  prodp++;
+
+  /* For each iteration in the outer loop, multiply one limb from
+     U with one limb from V, and add it to PROD.  */
+  for (i = 1; i < size; i++)
+    {
+      v_limb = up[i];
+      if (v_limb <= 1)
+       {
+         cy_limb = 0;
+         if (v_limb == 1)
+           cy_limb = __mpn_add_n (prodp, prodp, up, size);
+       }
+      else
+       cy_limb = __mpn_addmul_1 (prodp, up, size, v_limb);
+
+      prodp[size] = cy_limb;
+      prodp++;
+    }
+}
+
+void
+#if __STDC__
+____mpn_sqr_n (mp_ptr prodp,
+            mp_srcptr up, mp_size_t size, mp_ptr tspace)
+#else
+____mpn_sqr_n (prodp, up, size, tspace)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_size_t size;
+     mp_ptr tspace;
+#endif
+{
+  if ((size & 1) != 0)
+    {
+      /* The size is odd, the code code below doesn't handle that.
+        Multiply the least significant (size - 1) limbs with a recursive
+        call, and handle the most significant limb of S1 and S2
+        separately.  */
+      /* A slightly faster way to do this would be to make the Karatsuba
+        code below behave as if the size were even, and let it check for
+        odd size in the end.  I.e., in essence move this code to the end.
+        Doing so would save us a recursive call, and potentially make the
+        stack grow a lot less.  */
+
+      mp_size_t esize = size - 1;      /* even size */
+      mp_limb cy_limb;
+
+      MPN_SQR_N_RECURSE (prodp, up, esize, tspace);
+      cy_limb = __mpn_addmul_1 (prodp + esize, up, esize, up[esize]);
+      prodp[esize + esize] = cy_limb;
+      cy_limb = __mpn_addmul_1 (prodp + esize, up, size, up[esize]);
+
+      prodp[esize + size] = cy_limb;
+    }
+  else
+    {
+      mp_size_t hsize = size >> 1;
+      mp_limb cy;
+
+      /*** Product H.   ________________  ________________
+                       |_____U1 x U1____||____U0 x U0_____|  */
+      /* Put result in upper part of PROD and pass low part of TSPACE
+        as new TSPACE.  */
+      MPN_SQR_N_RECURSE (prodp + size, up + hsize, hsize, tspace);
+
+      /*** Product M.   ________________
+                       |_(U1-U0)(U0-U1)_|  */
+      if (__mpn_cmp (up + hsize, up, hsize) >= 0)
+       {
+         __mpn_sub_n (prodp, up + hsize, up, hsize);
+       }
+      else
+       {
+         __mpn_sub_n (prodp, up, up + hsize, hsize);
+       }
+
+      /* Read temporary operands from low part of PROD.
+        Put result in low part of TSPACE using upper part of TSPACE
+        as new TSPACE.  */
+      MPN_SQR_N_RECURSE (tspace, prodp, hsize, tspace + size);
+
+      /*** Add/copy product H.  */
+      MPN_COPY (prodp + hsize, prodp + size, hsize);
+      cy = __mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
+
+      /*** Add product M (if NEGFLG M is a negative number).  */
+      cy -= __mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
+
+      /*** Product L.   ________________  ________________
+                       |________________||____U0 x U0_____|  */
+      /* Read temporary operands from low part of PROD.
+        Put result in low part of TSPACE using upper part of TSPACE
+        as new TSPACE.  */
+      MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
+
+      /*** Add/copy Product L (twice).  */
+
+      cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
+      if (cy)
+       {
+         if (cy > 0)
+           __mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
+         else
+           {
+             __mpn_sub_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
+             abort ();
+           }
+       }
+
+      MPN_COPY (prodp, tspace, hsize);
+      cy = __mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
+      if (cy)
+       __mpn_add_1 (prodp + size, prodp + size, size, 1);
+    }
+}
+
+/* This should be made into an inline function in gmp.h.  */
+inline void
+#if __STDC__
+__mpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
+#else
+__mpn_mul_n (prodp, up, vp, size)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_srcptr vp;
+     mp_size_t size;
+#endif
+{
+  if (up == vp)
+    {
+      if (size < KARATSUBA_THRESHOLD)
+       {
+         ____mpn_sqr_n_basecase (prodp, up, size);
+       }
+      else
+       {
+         mp_ptr tspace;
+         tspace = (mp_ptr) alloca (2 * size * BYTES_PER_MP_LIMB);
+         ____mpn_sqr_n (prodp, up, size, tspace);
+       }
+    }
+  else
+    {
+      if (size < KARATSUBA_THRESHOLD)
+       {
+         ____mpn_mul_n_basecase (prodp, up, vp, size);
+       }
+      else
+       {
+         mp_ptr tspace;
+         tspace = (mp_ptr) alloca (2 * size * BYTES_PER_MP_LIMB);
+         ____mpn_mul_n (prodp, up, vp, size, tspace);
+       }
+    }
+}
diff --git a/sysdeps/generic/rshift.c b/sysdeps/generic/rshift.c
new file mode 100644 (file)
index 0000000..966cc7b
--- /dev/null
@@ -0,0 +1,87 @@
+/* __mpn_rshift -- Shift right a low-level natural-number integer.
+
+Copyright (C) 1991, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
+   and store the USIZE least significant limbs of the result at WP.
+   The bits shifted out to the right are returned.
+
+   Argument constraints:
+   1. 0 < CNT < BITS_PER_MP_LIMB
+   2. If the result is to be written over the input, WP must be <= UP.
+*/
+
+mp_limb
+#if __STDC__
+__mpn_rshift (register mp_ptr wp,
+           register mp_srcptr up, mp_size_t usize,
+           register unsigned int cnt)
+#else
+__mpn_rshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+#endif
+{
+  register mp_limb high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb retval;
+
+#ifdef DEBUG
+  if (usize == 0 || cnt == 0)
+    abort ();
+#endif
+
+  sh_1 = cnt;
+
+#if 0
+  if (sh_1 == 0)
+    {
+      if (wp != up)
+       {
+         /* Copy from low end to high end, to allow specified input/output
+            overlapping.  */
+         for (i = 0; i < usize; i++)
+           wp[i] = up[i];
+       }
+      return usize;
+    }
+#endif
+
+  wp -= 1;
+  sh_2 = BITS_PER_MP_LIMB - sh_1;
+  high_limb = up[0];
+  retval = high_limb << sh_2;
+  low_limb = high_limb;
+
+  for (i = 1; i < usize; i++)
+    {
+      high_limb = up[i];
+      wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
+      low_limb = high_limb;
+    }
+  wp[i] = low_limb >> sh_1;
+
+  return retval;
+}
diff --git a/sysdeps/generic/sub_n.c b/sysdeps/generic/sub_n.c
new file mode 100644 (file)
index 0000000..6b33e66
--- /dev/null
@@ -0,0 +1,61 @@
+/* __mpn_sub_n -- Subtract two limb vectors of equal, non-zero length.
+
+Copyright (C) 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mp_limb
+#if __STDC__
+__mpn_sub_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+__mpn_sub_n (res_ptr, s1_ptr, s2_ptr, size)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     register mp_srcptr s2_ptr;
+     mp_size_t size;
+#endif
+{
+  register mp_limb x, y, cy;
+  register mp_size_t j;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  s2_ptr -= j;
+  res_ptr -= j;
+
+  cy = 0;
+  do
+    {
+      y = s2_ptr[j];
+      x = s1_ptr[j];
+      y += cy;                 /* add previous carry to subtrahend */
+      cy = (y < cy);           /* get out carry from that addition */
+      y = x - y;               /* main subtract */
+      cy = (y > x) + cy;       /* get out carry from the subtract, combine */
+      res_ptr[j] = y;
+    }
+  while (++j != 0);
+
+  return cy;
+}
diff --git a/sysdeps/generic/submul_1.c b/sysdeps/generic/submul_1.c
new file mode 100644 (file)
index 0000000..855dd3f
--- /dev/null
@@ -0,0 +1,64 @@
+/* __mpn_submul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
+   by S2_LIMB, subtract the S1_SIZE least significant limbs of the product
+   from the limb vector pointed to by RES_PTR.  Return the most significant
+   limb of the product, adjusted for carry-out from the subtraction.
+
+Copyright (C) 1992, 1993, 1994 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 Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+mp_limb
+__mpn_submul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb s2_limb;
+{
+  register mp_limb cy_limb;
+  register mp_size_t j;
+  register mp_limb prod_high, prod_low;
+  register mp_limb x;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  res_ptr -= j;
+  s1_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      x = res_ptr[j];
+      prod_low = x - prod_low;
+      cy_limb += (prod_low > x);
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}