1 files changed, 118 insertions, 0 deletions
diff --git a/gnu/lib/libgmp/mpz_perfsqr.c b/gnu/lib/libgmp/mpz_perfsqr.c
new file mode 100644
index 0000000..4f14c06
--- /dev/null
+++ b/gnu/lib/libgmp/mpz_perfsqr.c
@@ -0,0 +1,118 @@
+/* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a pefect square,
+   zero otherwise.
+
+Copyright (C) 1991 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the GNU MP Library; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#if BITS_PER_MP_LIMB == 32
+static unsigned int primes[] = {3, 5, 7, 11, 13, 17, 19, 23, 29};
+static unsigned long int residue_map[] =
+{0x3, 0x13, 0x17, 0x23b, 0x161b, 0x1a317, 0x30af3, 0x5335f, 0x13d122f3};
+
+#define PP 0xC0CFD797L		/* 3 x 5 x 7 x 11 x 13 x ... x 29 */
+#endif
+
+/* sq_res_0x100[x mod 0x100] == 1 iff x mod 0x100 is a quadratic residue
+   modulo 0x100.  */
+static char sq_res_0x100[0x100] =
+{
+  1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+  0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
+};
+
+int
+#ifdef __STDC__
+mpz_perfect_square_p (const MP_INT *a)
+#else
+mpz_perfect_square_p (a)
+     const MP_INT *a;
+#endif
+{
+  mp_limb n1, n0;
+  mp_size i;
+  mp_size asize = a->size;
+  mp_srcptr aptr = a->d;
+  mp_limb rem;
+  mp_ptr root_ptr;
+
+  /* No negative numbers are perfect squares.  */
+  if (asize < 0)
+    return 0;
+
+  /* The first test excludes 55/64 (85.9%) of the perfect square candidates
+     in O(1) time.  */
+  if (sq_res_0x100[aptr[0] % 0x100] == 0)
+    return 0;
+
+#if BITS_PER_MP_LIMB == 32
+  /* The second test excludes 30652543/30808063 (99.5%) of the remaining
+     perfect square candidates in O(n) time.  */
+
+  /* Firstly, compute REM = A mod PP.  */
+  n1 = aptr[asize - 1];
+  if (n1 >= PP)
+    {
+      n1 = 0;
+      i = asize - 1;
+    }
+  else
+    i = asize - 2;
+
+  for (; i >= 0; i--)
+    {
+      mp_limb dummy;
+
+      n0 = aptr[i];
+      udiv_qrnnd (dummy, n1, n1, n0, PP);
+    }
+  rem = n1;
+
+  /* We have A mod PP in REM.  Now decide if REM is a quadratic residue
+     modulo the factors in PP.  */
+  for (i = 0; i < (sizeof primes) / sizeof (int); i++)
+    {
+      unsigned int p;
+
+      p = primes[i];
+      rem %= p;
+      if ((residue_map[i] & (1L << rem)) == 0)
+	return 0;
+    }
+#endif
+
+  /* For the third and last test, we finally compute the square root,
+     to make sure we've really got a perfect square.  */
+  root_ptr = (mp_ptr) alloca ((asize + 1) / 2 * BYTES_PER_MP_LIMB);
+
+  /* Iff mpn_sqrt returns zero, the square is perfect.  */
+  {
+    int retval = !mpn_sqrt (root_ptr, NULL, aptr, asize);
+    alloca (0);
+    return retval;
+  }
+}