Clean import of libgmp 2.0.2, with only the non-x86 bits removed.

BMakefiles and other bits will follow.

Requested by:	Andrey Chernov
Made world by:	Chuck Robey

10 files changed, 1941 insertions, 0 deletions

diff --git a/contrib/libgmp/mpn/tests/add_n.c b/contrib/libgmp/mpn/tests/add_n.c
new file mode 100644
index 0000000..c27d347
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/add_n.c
@@ -0,0 +1,211 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 10000000
+#endif
+#ifndef SIZE
+#define SIZE 328
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+
+mp_limb_t
+#if __STDC__
+refmpn_add_n (mp_ptr res_ptr,
+	       mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+refmpn_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_t 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;
+}
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t s2[SIZE];
+  mp_limb_t dx[SIZE+1];
+  mp_limb_t dy[SIZE+1];
+  int cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+
+      mpn_random2 (s1, size);
+      mpn_random2 (s2, size);
+
+      dx[size] = 0x12345678;
+      dy[size] = 0x12345678;
+
+#ifdef PRINT
+      mpn_print (s1, size);
+      mpn_print (s2, size);
+#endif
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_add_n (dx, s1, s2, size);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_add_n:   %ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%d ", cyx); mpn_print (dx, size);
+#endif
+
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = mpn_add_n (dx, s1, s2, size);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_add_n:   %ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%d ", cyx); mpn_print (dx, size);
+#endif
+
+#ifndef NOCHECK
+      /* Put garbage in the destination.  */
+      for (i = 0; i < size; i++)
+	{
+	  dx[i] = 0x7654321;
+	  dy[i] = 0x1234567;
+	}
+
+      cyx = refmpn_add_n (dx, s1, s2, size);
+      cyy = mpn_add_n (dy, s1, s2, size);
+      if (cyx != cyy || mpn_cmp (dx, dy, size) != 0
+	  || dx[size] != 0x12345678 || dy[size] != 0x12345678)
+	{
+#ifndef PRINT
+	  printf ("%d ", cyx); mpn_print (dx, size);
+	  printf ("%d ", cyy); mpn_print (dy, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/addmul_1.c b/contrib/libgmp/mpn/tests/addmul_1.c
new file mode 100644
index 0000000..23952a1
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/addmul_1.c
@@ -0,0 +1,223 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 20000000
+#endif
+#ifndef SIZE
+#define SIZE 496
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+mp_limb_t
+refmpn_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_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t prod_high, prod_low;
+  register mp_limb_t 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;
+}
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t dx[SIZE+2];
+  mp_limb_t dy[SIZE+2];
+  mp_limb_t cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  mp_limb_t xlimb;
+  mp_size_t size;
+  double cyc;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+
+      mpn_random2 (s1, size);
+      mpn_random2 (dy+1, size);
+
+      if (random () % 0x100 == 0)
+	xlimb = 0;
+      else
+	mpn_random2 (&xlimb, 1);
+
+      dy[size+1] = 0x12345678;
+      dy[0] = 0x87654321;
+
+#if defined (PRINT) || defined (XPRINT)
+      printf ("xlimb=%*lX\n", (int) (2 * sizeof(mp_limb_t)), xlimb);
+#endif
+#ifdef PRINT
+      mpn_print (dy+1, size);
+      mpn_print (s1, size);
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_addmul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      cyc = ((double) t * CLOCK) / (OPS * 1000.0);
+      printf ("refmpn_addmul_1: %5ldms (%.2f cycles/limb) [%.2f Gb/s]\n",
+	      t,
+	      cyc,
+	      CLOCK/cyc*BITS_PER_MP_LIMB*BITS_PER_MP_LIMB);
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyy = mpn_addmul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      cyc = ((double) t * CLOCK) / (OPS * 1000.0);
+      printf ("mpn_addmul_1:    %5ldms (%.2f cycles/limb) [%.2f Gb/s]\n",
+	      t,
+	      cyc,
+	      CLOCK/cyc*BITS_PER_MP_LIMB*BITS_PER_MP_LIMB);
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      cyx = refmpn_addmul_1 (dx+1, s1, size, xlimb);
+      cyy = mpn_addmul_1 (dy+1, s1, size, xlimb);
+
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+      mpn_print (dx+1, size);
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+      mpn_print (dy+1, size);
+#endif
+
+#ifndef NOCHECK
+      if (cyx != cyy || mpn_cmp (dx, dy, size+2) != 0
+	  || dx[size+1] != 0x12345678 || dx[0] != 0x87654321)
+	{
+#ifndef PRINT
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+	  mpn_print (dx+1, size);
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+	  mpn_print (dy+1, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/divmod_1.c b/contrib/libgmp/mpn/tests/divmod_1.c
new file mode 100644
index 0000000..f6b541e
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/divmod_1.c
@@ -0,0 +1,120 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 20000000
+#endif
+#ifndef SIZE
+#define SIZE 1000
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+main ()
+{
+  mp_limb_t nptr[SIZE];
+  mp_limb_t qptr[SIZE];
+  mp_limb_t pptr[SIZE];
+  mp_limb_t dlimb, rlimb, plimb;
+  mp_size_t nsize, qsize, psize;
+  int test;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      nsize = random () % SIZE + 1;
+#else
+      nsize = SIZE;
+#endif
+
+      mpn_random2 (nptr, nsize);
+
+      mpn_random2 (&dlimb, 1);
+      if (dlimb == 0)
+	abort ();
+
+      rlimb = mpn_divmod_1 (qptr, nptr, nsize, dlimb);
+      qsize = nsize - (qptr[nsize - 1] == 0);
+      if (qsize == 0)
+	{
+	  plimb = rlimb;
+	  psize = qsize;
+	}
+      else
+	{
+	  plimb = mpn_mul_1 (pptr, qptr, qsize, dlimb);
+	  psize = qsize;
+	  plimb += mpn_add_1 (pptr, pptr, psize, rlimb);
+	}
+      if (plimb != 0)
+	pptr[psize++] = plimb;
+
+
+      if (nsize != psize || mpn_cmp (nptr, pptr, nsize) != 0)
+	abort ();
+    }
+}
diff --git a/contrib/libgmp/mpn/tests/divrem.c b/contrib/libgmp/mpn/tests/divrem.c
new file mode 100644
index 0000000..6eafc99
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/divrem.c
@@ -0,0 +1,129 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 20000000
+#endif
+#ifndef SIZE
+#define SIZE 100
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+main ()
+{
+  mp_limb_t nptr[2 * SIZE];
+  mp_limb_t dptr[SIZE];
+  mp_limb_t qptr[2 * SIZE];
+  mp_limb_t pptr[2 * SIZE];
+  mp_limb_t rptr[2 * SIZE];
+  mp_size_t nsize, dsize, qsize, rsize, psize;
+  int test;
+  mp_limb_t qlimb;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      nsize = random () % (2 * SIZE) + 1;
+      dsize = random () % nsize + 1;
+#else
+      nsize = 2 * SIZE;
+      dsize = SIZE;
+#endif
+
+      mpn_random2 (nptr, nsize);
+      mpn_random2 (dptr, dsize);
+      dptr[dsize - 1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
+
+      MPN_COPY (rptr, nptr, nsize);
+      qlimb = mpn_divrem (qptr, (mp_size_t) 0, rptr, nsize, dptr, dsize);
+      rsize = dsize;
+      qsize = nsize - dsize;
+      qptr[qsize] = qlimb;
+      qsize += qlimb;
+      if (qsize == 0 || qsize > 2 * SIZE)
+	{
+	  continue;		/* bogus */
+	}
+      else
+	{
+	  mp_limb_t cy;
+	  if (qsize > dsize)
+	    mpn_mul (pptr, qptr, qsize, dptr, dsize);
+	  else
+	    mpn_mul (pptr, dptr, dsize, qptr, qsize);
+	  psize = qsize + dsize;
+	  psize -= pptr[psize - 1] == 0;
+	  cy = mpn_add (pptr, pptr, psize, rptr, rsize);
+	  pptr[psize] = cy;
+	  psize += cy;
+	}
+
+      if (nsize != psize || mpn_cmp (nptr, pptr, nsize) != 0)
+	abort ();
+    }
+}
diff --git a/contrib/libgmp/mpn/tests/lshift.c b/contrib/libgmp/mpn/tests/lshift.c
new file mode 100644
index 0000000..f50c5dc
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/lshift.c
@@ -0,0 +1,226 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 10000000
+#endif
+#ifndef SIZE
+#define SIZE 496
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+mp_limb_t
+refmpn_lshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+{
+  register mp_limb_t high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb_t 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;
+}
+
+#ifndef CNT
+#define CNT 4
+#endif
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t dx[SIZE+2];
+  mp_limb_t dy[SIZE+2];
+  mp_limb_t cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  int cnt = CNT;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+      mpn_random2 (s1, size);
+
+      dx[size+1] = 0x12345678;
+      dy[size+1] = 0x12345678;
+      dx[0] = 0x87654321;
+      dy[0] = 0x87654321;
+
+#ifdef PRINT
+      mpn_print (s1, size);
+#endif
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_lshift (dx+1, s1, size, cnt);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_lshift: %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx); mpn_print (dx+1, size);
+#endif
+
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyy = mpn_lshift (dx+1, s1, size, cnt);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_lshift:  %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy); mpn_print (dx+1, size);
+#endif
+
+#ifndef NOCHECK
+      /* Put garbage in the destination.  */
+      for (i = 1; i <= size; i++)
+	{
+	  dx[i] = 0x7654321;
+	  dy[i] = 0x1234567;
+	}
+
+      cyx = refmpn_lshift (dx+1, s1, size, cnt);
+      cyy = mpn_lshift (dy+1, s1, size, cnt);
+
+      if (cyx != cyy || mpn_cmp (dx, dy, size+2) != 0
+	  || dx[size+1] != 0x12345678 || dx[0] != 0x87654321)
+	{
+#ifndef PRINT
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+	  mpn_print (dx+1, size);
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+	  mpn_print (dy+1, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/mul_1.c b/contrib/libgmp/mpn/tests/mul_1.c
new file mode 100644
index 0000000..2b522fa
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/mul_1.c
@@ -0,0 +1,212 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 20000000
+#endif
+#ifndef SIZE
+#define SIZE 496
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+mp_limb_t
+refmpn_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_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t 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;
+}
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t dx[SIZE+2];
+  mp_limb_t dy[SIZE+2];
+  mp_limb_t cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  mp_limb_t xlimb;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+
+      mpn_random2 (s1, size);
+      mpn_random2 (dy+1, size);
+
+      if (random () % 0x100 == 0)
+	xlimb = 0;
+      else
+	mpn_random2 (&xlimb, 1);
+
+      dy[size+1] = 0x12345678;
+      dy[0] = 0x87654321;
+
+#if defined (PRINT) || defined (XPRINT)
+      printf ("xlimb=%*lX\n", (int) (2 * sizeof(mp_limb_t)), xlimb);
+#endif
+#ifdef PRINT
+      mpn_print (s1, size);
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_mul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_mul_1: %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyy = mpn_mul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_mul_1:    %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+
+      cyx = refmpn_mul_1 (dx+1, s1, size, xlimb);
+      cyy = mpn_mul_1 (dy+1, s1, size, xlimb);
+
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+      mpn_print (dx+1, size);
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+      mpn_print (dy+1, size);
+#endif
+
+#ifndef NOCHECK
+      if (cyx != cyy || mpn_cmp (dx, dy, size+2) != 0
+	  || dx[size+1] != 0x12345678 || dx[0] != 0x87654321)
+	{
+#ifndef PRINT
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+	  mpn_print (dx+1, size);
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+	  mpn_print (dy+1, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/rshift.c b/contrib/libgmp/mpn/tests/rshift.c
new file mode 100644
index 0000000..2482bf3
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/rshift.c
@@ -0,0 +1,227 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 10000000
+#endif
+#ifndef SIZE
+#define SIZE 496
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+mp_limb_t
+refmpn_rshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+{
+  register mp_limb_t high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb_t 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 0;
+    }
+#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;
+    }
+  low_limb >>= sh_1;
+  wp[i] = low_limb;
+
+  return retval;
+}
+
+#ifndef CNT
+#define CNT 4
+#endif
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t dx[SIZE+2];
+  mp_limb_t dy[SIZE+2];
+  mp_limb_t cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  int cnt = CNT;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+      mpn_random2 (s1, size);
+
+      dx[size+1] = 0x12345678;
+      dy[size+1] = 0x12345678;
+      dx[0] = 0x87654321;
+      dy[0] = 0x87654321;
+
+#ifdef PRINT
+      mpn_print (s1, size);
+#endif
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_rshift (dx+1, s1, size, cnt);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_rshift: %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx); mpn_print (dx+1, size);
+#endif
+
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyy = mpn_rshift (dx+1, s1, size, cnt);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_rshift:  %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy); mpn_print (dx+1, size);
+#endif
+
+#ifndef NOCHECK
+      /* Put garbage in the destination.  */
+      for (i = 1; i <= size; i++)
+	{
+	  dx[i] = 0x7654321;
+	  dy[i] = 0x1234567;
+	}
+
+      cyx = refmpn_rshift (dx+1, s1, size, cnt);
+      cyy = mpn_rshift (dy+1, s1, size, cnt);
+
+      if (cyx != cyy || mpn_cmp (dx, dy, size+2) != 0
+	  || dx[size+1] != 0x12345678 || dx[0] != 0x87654321)
+	{
+#ifndef PRINT
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+	  mpn_print (dx+1, size);
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+	  mpn_print (dy+1, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/sub_n.c b/contrib/libgmp/mpn/tests/sub_n.c
new file mode 100644
index 0000000..2b9031b
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/sub_n.c
@@ -0,0 +1,211 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 10000000
+#endif
+#ifndef SIZE
+#define SIZE 328
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+
+mp_limb_t
+#if __STDC__
+refmpn_sub_n (mp_ptr res_ptr,
+	       mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+refmpn_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_t 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;
+}
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t s2[SIZE];
+  mp_limb_t dx[SIZE+1];
+  mp_limb_t dy[SIZE+1];
+  int cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+
+      mpn_random2 (s1, size);
+      mpn_random2 (s2, size);
+
+      dx[size] = 0x12345678;
+      dy[size] = 0x12345678;
+
+#ifdef PRINT
+      mpn_print (s1, size);
+      mpn_print (s2, size);
+#endif
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_sub_n (dx, s1, s2, size);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_sub_n:   %ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%d ", cyx); mpn_print (dx, size);
+#endif
+
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = mpn_sub_n (dx, s1, s2, size);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_sub_n:   %ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+#ifdef PRINT
+      printf ("%d ", cyx); mpn_print (dx, size);
+#endif
+
+#ifndef NOCHECK
+      /* Put garbage in the destination.  */
+      for (i = 0; i < size; i++)
+	{
+	  dx[i] = 0x7654321;
+	  dy[i] = 0x1234567;
+	}
+
+      cyx = refmpn_sub_n (dx, s1, s2, size);
+      cyy = mpn_sub_n (dy, s1, s2, size);
+      if (cyx != cyy || mpn_cmp (dx, dy, size) != 0
+	  || dx[size] != 0x12345678 || dy[size] != 0x12345678)
+	{
+#ifndef PRINT
+	  printf ("%d ", cyx); mpn_print (dx, size);
+	  printf ("%d ", cyy); mpn_print (dy, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/submul_1.c b/contrib/libgmp/mpn/tests/submul_1.c
new file mode 100644
index 0000000..0e464e6
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/submul_1.c
@@ -0,0 +1,218 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef USG
+#include <sys/time.h>
+#include <sys/resource.h>
+
+unsigned long
+cputime ()
+{
+    struct rusage rus;
+
+    getrusage (0, &rus);
+    return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#include <time.h>
+
+#ifndef CLOCKS_PER_SEC
+#define CLOCKS_PER_SEC 1000000
+#endif
+
+#if CLOCKS_PER_SEC >= 10000
+#define CLOCK_TO_MILLISEC(cl) ((cl) / (CLOCKS_PER_SEC / 1000))
+#else
+#define CLOCK_TO_MILLISEC(cl) ((cl) * 1000 / CLOCKS_PER_SEC)
+#endif
+
+unsigned long
+cputime ()
+{
+  return CLOCK_TO_MILLISEC (clock ());
+}
+#endif
+
+#define M * 1000000
+
+#ifndef CLOCK
+#if defined (__m88k__)
+#define CLOCK 20 M
+#elif defined (__i386__)
+#define CLOCK (16.666667 M)
+#elif defined (__m68k__)
+#define CLOCK (20 M)
+#elif defined (_IBMR2)
+#define CLOCK (25 M)
+#elif defined (__sparc__)
+#define CLOCK (20 M)
+#elif defined (__sun__)
+#define CLOCK (20 M)
+#elif defined (__mips)
+#define CLOCK (40 M)
+#elif defined (__hppa__)
+#define CLOCK (50 M)
+#elif defined (__alpha)
+#define CLOCK (133 M)
+#else
+#error "Don't know CLOCK of your machine"
+#endif
+#endif
+
+#ifndef OPS
+#define OPS 20000000
+#endif
+#ifndef SIZE
+#define SIZE 496
+#endif
+#ifndef TIMES
+#define TIMES OPS/SIZE
+#else
+#undef OPS
+#define OPS (SIZE*TIMES)
+#endif
+
+mp_limb_t
+refmpn_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_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t prod_high, prod_low;
+  register mp_limb_t 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;
+}
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_limb_t s1[SIZE];
+  mp_limb_t dx[SIZE+2];
+  mp_limb_t dy[SIZE+2];
+  mp_limb_t cyx, cyy;
+  int i;
+  long t0, t;
+  int test;
+  mp_limb_t xlimb;
+  mp_size_t size;
+
+  for (test = 0; ; test++)
+    {
+#ifdef RANDOM
+      size = (random () % SIZE + 1);
+#else
+      size = SIZE;
+#endif
+
+      mpn_random2 (s1, size);
+      mpn_random2 (dy+1, size);
+
+      if (random () % 0x100 == 0)
+	xlimb = 0;
+      else
+	mpn_random2 (&xlimb, 1);
+
+      dy[size+1] = 0x12345678;
+      dy[0] = 0x87654321;
+
+#if defined (PRINT) || defined (XPRINT)
+      printf ("xlimb=%*lX\n", (int) (2 * sizeof(mp_limb_t)), xlimb);
+#endif
+#ifdef PRINT
+      mpn_print (dy+1, size);
+      mpn_print (s1, size);
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyx = refmpn_submul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("refmpn_submul_1: %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      t0 = cputime();
+      for (i = 0; i < TIMES; i++)
+	cyy = mpn_submul_1 (dx+1, s1, size, xlimb);
+      t = cputime() - t0;
+#if TIMES != 1
+      printf ("mpn_submul_1:    %5ldms (%.2f cycles/limb)\n",
+	      t,
+	      ((double) t * CLOCK) / (OPS * 1000.0));
+#endif
+
+      MPN_COPY (dx, dy, size+2);
+      cyx = refmpn_submul_1 (dx+1, s1, size, xlimb);
+      cyy = mpn_submul_1 (dy+1, s1, size, xlimb);
+
+#ifdef PRINT
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+      mpn_print (dx+1, size);
+      printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+      mpn_print (dy+1, size);
+#endif
+
+#ifndef NOCHECK
+      if (cyx != cyy || mpn_cmp (dx, dy, size+2) != 0
+	  || dx[size+1] != 0x12345678 || dx[0] != 0x87654321)
+	{
+#ifndef PRINT
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyx);
+	  mpn_print (dx+1, size);
+	  printf ("%*lX ", (int) (2 * sizeof(mp_limb_t)), cyy);
+	  mpn_print (dy+1, size);
+#endif
+	  abort();
+	}
+#endif
+    }
+}
+
+mpn_print (mp_ptr p, mp_size_t size)
+{
+  mp_size_t i;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      printf ("%0*lX", (int) (2 * sizeof(mp_limb_t)), p[i]);
+#ifdef SPACE
+      if (i != 0)
+	printf (" ");
+#endif
+    }
+  puts ("");
+}
diff --git a/contrib/libgmp/mpn/tests/tst-addsub.c b/contrib/libgmp/mpn/tests/tst-addsub.c
new file mode 100644
index 0000000..e02b9d5
--- /dev/null
+++ b/contrib/libgmp/mpn/tests/tst-addsub.c
@@ -0,0 +1,164 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+
+#define ADD 1
+#define SUB 2
+
+#ifndef METHOD
+#define METHOD ADD
+#endif
+
+#if METHOD == ADD
+#define REFCALL refmpn_add_n
+#define TESTCALL mpn_add_n
+#endif
+
+#if METHOD == SUB
+#define REFCALL refmpn_sub_n
+#define TESTCALL mpn_sub_n
+#endif
+
+mp_limb_t refmpn_add_n ();
+mp_limb_t refmpn_sub_n ();
+
+#define SIZE 100
+
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  mp_size_t alloc_size, max_size, size, i, cumul_size;
+  mp_ptr s1, s2, dx, dy;
+  int s1_align, s2_align, d_align;
+  long pass, n_passes;
+  mp_limb_t cx, cy;
+
+  max_size = SIZE;
+  n_passes = 1000000;
+
+  argc--; argv++;
+  if (argc)
+    {
+      max_size = atol (*argv);
+      argc--; argv++;
+    }
+
+  alloc_size = max_size + 32;
+  s1 = malloc (alloc_size * BYTES_PER_MP_LIMB);
+  s2 = malloc (alloc_size * BYTES_PER_MP_LIMB);
+  dx = malloc (alloc_size * BYTES_PER_MP_LIMB);
+  dy = malloc (alloc_size * BYTES_PER_MP_LIMB);
+
+  cumul_size = 0;
+  for (pass = 0; pass < n_passes; pass++)
+    {
+      cumul_size += size;
+      if (cumul_size >= 1000000)
+	{
+	  cumul_size -= 1000000;
+	  printf ("%d ", pass); fflush (stdout);
+	}
+      s1_align = random () % 32;
+      s2_align = random () % 32;
+      d_align = random () % 32;
+
+      size = random () % max_size + 1;
+
+      mpn_random2 (s1 + s1_align, size);
+      mpn_random2 (s2 + s2_align, size);
+
+      for (i = 0; i < alloc_size; i++)
+	dx[i] = dy[i] = i + 0x9876500;
+
+      cx = TESTCALL (dx + d_align, s1 + s1_align, s2 + s2_align, size);
+      cy = REFCALL (dy + d_align, s1 + s1_align, s2 + s2_align, size);
+
+      if (cx != cy || mpn_cmp (dx, dy, alloc_size) != 0)
+	abort ();
+    }
+
+  printf ("%d passes OK\n", n_passes);
+  exit (0);
+}
+
+mp_limb_t
+#if __STDC__
+refmpn_add_n (mp_ptr res_ptr,
+	      mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+refmpn_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_t 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;
+}
+
+mp_limb_t
+#if __STDC__
+refmpn_sub_n (mp_ptr res_ptr,
+	       mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+refmpn_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_t 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;
+}