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