From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/mpn/generic/perfpow.c | 342 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 342 insertions(+) create mode 100644 gmp-6.3.0/mpn/generic/perfpow.c (limited to 'gmp-6.3.0/mpn/generic/perfpow.c') 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; +} -- cgit v1.2.3