Initial commit.

1 files changed, 342 insertions, 0 deletions

diff --git a/gmp-6.3.0/mpn/generic/perfpow.c b/gmp-6.3.0/mpn/generic/perfpow.c
new file mode 100644
index 0000000..9d46477
--- /dev/null
+++ b/gmp-6.3.0/mpn/generic/perfpow.c
@@ -0,0 +1,342 @@
+/* mpn_perfect_power_p -- mpn perfect power detection.
+
+   Contributed to the GNU project by Martin Boij.
+
+Copyright 2009, 2010, 2012, 2014 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 either:
+
+  * the GNU Lesser General Public License as published by the Free
+    Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+or
+
+  * the GNU General Public License as published by the Free Software
+    Foundation; either version 2 of the License, or (at your option) any
+    later version.
+
+or both in parallel, as here.
+
+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 copies of the GNU General Public License and the
+GNU Lesser General Public License along with the GNU MP Library.  If not,
+see https://www.gnu.org/licenses/.  */
+
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#define SMALL 20
+#define MEDIUM 100
+
+/* Return non-zero if {np,nn} == {xp,xn} ^ k.
+   Algorithm:
+       For s = 1, 2, 4, ..., s_max, compute the s least significant limbs of
+       {xp,xn}^k. Stop if they don't match the s least significant limbs of
+       {np,nn}.
+
+   FIXME: Low xn limbs can be expected to always match, if computed as a mod
+   B^{xn} root. So instead of using mpn_powlo, compute an approximation of the
+   most significant (normalized) limb of {xp,xn} ^ k (and an error bound), and
+   compare to {np, nn}. Or use an even cruder approximation based on fix-point
+   base 2 logarithm.  */
+static int
+pow_equals (mp_srcptr np, mp_size_t n,
+	    mp_srcptr xp,mp_size_t xn,
+	    mp_limb_t k, mp_bitcnt_t f,
+	    mp_ptr tp)
+{
+  mp_bitcnt_t y, z;
+  mp_size_t bn;
+  mp_limb_t h, l;
+
+  ASSERT (n > 1 || (n == 1 && np[0] > 1));
+  ASSERT (np[n - 1] > 0);
+  ASSERT (xn > 0);
+
+  if (xn == 1 && xp[0] == 1)
+    return 0;
+
+  z = 1 + (n >> 1);
+  for (bn = 1; bn < z; bn <<= 1)
+    {
+      mpn_powlo (tp, xp, &k, 1, bn, tp + bn);
+      if (mpn_cmp (tp, np, bn) != 0)
+	return 0;
+    }
+
+  /* Final check. Estimate the size of {xp,xn}^k before computing the power
+     with full precision.  Optimization: It might pay off to make a more
+     accurate estimation of the logarithm of {xp,xn}, rather than using the
+     index of the MSB.  */
+
+  MPN_SIZEINBASE_2EXP(y, xp, xn, 1);
+  y -= 1;  /* msb_index (xp, xn) */
+
+  umul_ppmm (h, l, k, y);
+  h -= l == 0;  --l;	/* two-limb decrement */
+
+  z = f - 1; /* msb_index (np, n) */
+  if (h == 0 && l <= z)
+    {
+      mp_limb_t *tp2;
+      mp_size_t i;
+      int ans;
+      mp_limb_t size;
+      TMP_DECL;
+
+      size = l + k;
+      ASSERT_ALWAYS (size >= k);
+
+      TMP_MARK;
+      y = 2 + size / GMP_LIMB_BITS;
+      tp2 = TMP_ALLOC_LIMBS (y);
+
+      i = mpn_pow_1 (tp, xp, xn, k, tp2);
+      if (i == n && mpn_cmp (tp, np, n) == 0)
+	ans = 1;
+      else
+	ans = 0;
+      TMP_FREE;
+      return ans;
+    }
+
+  return 0;
+}
+
+
+/* Return non-zero if N = {np,n} is a kth power.
+   I = {ip,n} = N^(-1) mod B^n.  */
+static int
+is_kth_power (mp_ptr rp, mp_srcptr np,
+	      mp_limb_t k, mp_srcptr ip,
+	      mp_size_t n, mp_bitcnt_t f,
+	      mp_ptr tp)
+{
+  mp_bitcnt_t b;
+  mp_size_t rn, xn;
+
+  ASSERT (n > 0);
+  ASSERT ((k & 1) != 0 || k == 2);
+  ASSERT ((np[0] & 1) != 0);
+
+  if (k == 2)
+    {
+      b = (f + 1) >> 1;
+      rn = 1 + b / GMP_LIMB_BITS;
+      if (mpn_bsqrtinv (rp, ip, b, tp) != 0)
+	{
+	  rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
+	  xn = rn;
+	  MPN_NORMALIZE (rp, xn);
+	  if (pow_equals (np, n, rp, xn, k, f, tp) != 0)
+	    return 1;
+
+	  /* Check if (2^b - r)^2 == n */
+	  mpn_neg (rp, rp, rn);
+	  rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
+	  MPN_NORMALIZE (rp, rn);
+	  if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
+	    return 1;
+	}
+    }
+  else
+    {
+      b = 1 + (f - 1) / k;
+      rn = 1 + (b - 1) / GMP_LIMB_BITS;
+      mpn_brootinv (rp, ip, rn, k, tp);
+      if ((b % GMP_LIMB_BITS) != 0)
+	rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
+      MPN_NORMALIZE (rp, rn);
+      if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
+	return 1;
+    }
+  MPN_ZERO (rp, rn); /* Untrash rp */
+  return 0;
+}
+
+static int
+perfpow (mp_srcptr np, mp_size_t n,
+	 mp_limb_t ub, mp_limb_t g,
+	 mp_bitcnt_t f, int neg)
+{
+  mp_ptr ip, tp, rp;
+  mp_limb_t k;
+  int ans;
+  mp_bitcnt_t b;
+  gmp_primesieve_t ps;
+  TMP_DECL;
+
+  ASSERT (n > 0);
+  ASSERT ((np[0] & 1) != 0);
+  ASSERT (ub > 0);
+
+  TMP_MARK;
+  gmp_init_primesieve (&ps);
+  b = (f + 3) >> 1;
+
+  TMP_ALLOC_LIMBS_3 (ip, n, rp, n, tp, 5 * n);
+
+  MPN_ZERO (rp, n);
+
+  /* FIXME: It seems the inverse in ninv is needed only to get non-inverted
+     roots. I.e., is_kth_power computes n^{1/2} as (n^{-1})^{-1/2} and
+     similarly for nth roots. It should be more efficient to compute n^{1/2} as
+     n * n^{-1/2}, with a mullo instead of a binvert. And we can do something
+     similar for kth roots if we switch to an iteration converging to n^{1/k -
+     1}, and we can then eliminate this binvert call. */
+  mpn_binvert (ip, np, 1 + (b - 1) / GMP_LIMB_BITS, tp);
+  if (b % GMP_LIMB_BITS)
+    ip[(b - 1) / GMP_LIMB_BITS] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
+
+  if (neg)
+    gmp_nextprime (&ps);
+
+  ans = 0;
+  if (g > 0)
+    {
+      ub = MIN (ub, g + 1);
+      while ((k = gmp_nextprime (&ps)) < ub)
+	{
+	  if ((g % k) == 0)
+	    {
+	      if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
+		{
+		  ans = 1;
+		  goto ret;
+		}
+	    }
+	}
+    }
+  else
+    {
+      while ((k = gmp_nextprime (&ps)) < ub)
+	{
+	  if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
+	    {
+	      ans = 1;
+	      goto ret;
+	    }
+	}
+    }
+ ret:
+  TMP_FREE;
+  return ans;
+}
+
+static const unsigned short nrtrial[] = { 100, 500, 1000 };
+
+/* Table of (log_{p_i} 2) values, where p_i is the (nrtrial[i] + 1)'th prime
+   number.  */
+static const double logs[] =
+  { 0.1099457228193620, 0.0847016403115322, 0.0772048195144415 };
+
+int
+mpn_perfect_power_p (mp_srcptr np, mp_size_t n)
+{
+  mp_limb_t *nc, factor, g;
+  mp_limb_t exp, d;
+  mp_bitcnt_t twos, count;
+  int ans, where, neg, trial;
+  TMP_DECL;
+
+  neg = n < 0;
+  if (neg)
+    {
+      n = -n;
+    }
+
+  if (n == 0 || (n == 1 && np[0] == 1)) /* Valgrind doesn't like
+					   (n <= (np[0] == 1)) */
+    return 1;
+
+  TMP_MARK;
+
+  count = 0;
+
+  twos = mpn_scan1 (np, 0);
+  if (twos != 0)
+    {
+      mp_size_t s;
+      if (twos == 1)
+	{
+	  return 0;
+	}
+      s = twos / GMP_LIMB_BITS;
+      if (s + 1 == n && POW2_P (np[s]))
+	{
+	  return ! (neg && POW2_P (twos));
+	}
+      count = twos % GMP_LIMB_BITS;
+      n -= s;
+      np += s;
+      if (count > 0)
+	{
+	  nc = TMP_ALLOC_LIMBS (n);
+	  mpn_rshift (nc, np, n, count);
+	  n -= (nc[n - 1] == 0);
+	  np = nc;
+	}
+    }
+  g = twos;
+
+  trial = (n > SMALL) + (n > MEDIUM);
+
+  where = 0;
+  factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
+
+  if (factor != 0)
+    {
+      if (count == 0) /* We did not allocate nc yet. */
+	{
+	  nc = TMP_ALLOC_LIMBS (n);
+	}
+
+      /* Remove factors found by trialdiv.  Optimization: If remove
+	 define _itch, we can allocate its scratch just once */
+
+      do
+	{
+	  binvert_limb (d, factor);
+
+	  /* After the first round we always have nc == np */
+	  exp = mpn_remove (nc, &n, np, n, &d, 1, ~(mp_bitcnt_t)0);
+
+	  if (g == 0)
+	    g = exp;
+	  else
+	    g = mpn_gcd_1 (&g, 1, exp);
+
+	  if (g == 1)
+	    {
+	      ans = 0;
+	      goto ret;
+	    }
+
+	  if ((n == 1) & (nc[0] == 1))
+	    {
+	      ans = ! (neg && POW2_P (g));
+	      goto ret;
+	    }
+
+	  np = nc;
+	  factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
+	}
+      while (factor != 0);
+    }
+
+  MPN_SIZEINBASE_2EXP(count, np, n, 1);   /* log (np) + 1 */
+  d = (mp_limb_t) (count * logs[trial] + 1e-9) + 1;
+  ans = perfpow (np, n, d, g, count, neg);
+
+ ret:
+  TMP_FREE;
+  return ans;
+}