1/* mpn_perfect_power_p -- mpn perfect power detection.
2
3   Contributed to the GNU project by Martin Boij.
4
5Copyright 2009, 2010, 2012, 2014 Free Software Foundation, Inc.
6
7This file is part of the GNU MP Library.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of either:
11
12  * the GNU Lesser General Public License as published by the Free
13    Software Foundation; either version 3 of the License, or (at your
14    option) any later version.
15
16or
17
18  * the GNU General Public License as published by the Free Software
19    Foundation; either version 2 of the License, or (at your option) any
20    later version.
21
22or both in parallel, as here.
23
24The GNU MP Library is distributed in the hope that it will be useful, but
25WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
27for more details.
28
29You should have received copies of the GNU General Public License and the
30GNU Lesser General Public License along with the GNU MP Library.  If not,
31see https://www.gnu.org/licenses/.  */
32
33#include "gmp-impl.h"
34#include "longlong.h"
35
36#define SMALL 20
37#define MEDIUM 100
38
39/* Return non-zero if {np,nn} == {xp,xn} ^ k.
40   Algorithm:
41       For s = 1, 2, 4, ..., s_max, compute the s least significant limbs of
42       {xp,xn}^k. Stop if they don't match the s least significant limbs of
43       {np,nn}.
44
45   FIXME: Low xn limbs can be expected to always match, if computed as a mod
46   B^{xn} root. So instead of using mpn_powlo, compute an approximation of the
47   most significant (normalized) limb of {xp,xn} ^ k (and an error bound), and
48   compare to {np, nn}. Or use an even cruder approximation based on fix-point
49   base 2 logarithm.  */
50static int
51pow_equals (mp_srcptr np, mp_size_t n,
52	    mp_srcptr xp,mp_size_t xn,
53	    mp_limb_t k, mp_bitcnt_t f,
54	    mp_ptr tp)
55{
56  mp_bitcnt_t y, z;
57  mp_size_t bn;
58  mp_limb_t h, l;
59
60  ASSERT (n > 1 || (n == 1 && np[0] > 1));
61  ASSERT (np[n - 1] > 0);
62  ASSERT (xn > 0);
63
64  if (xn == 1 && xp[0] == 1)
65    return 0;
66
67  z = 1 + (n >> 1);
68  for (bn = 1; bn < z; bn <<= 1)
69    {
70      mpn_powlo (tp, xp, &k, 1, bn, tp + bn);
71      if (mpn_cmp (tp, np, bn) != 0)
72	return 0;
73    }
74
75  /* Final check. Estimate the size of {xp,xn}^k before computing the power
76     with full precision.  Optimization: It might pay off to make a more
77     accurate estimation of the logarithm of {xp,xn}, rather than using the
78     index of the MSB.  */
79
80  MPN_SIZEINBASE_2EXP(y, xp, xn, 1);
81  y -= 1;  /* msb_index (xp, xn) */
82
83  umul_ppmm (h, l, k, y);
84  h -= l == 0;  --l;	/* two-limb decrement */
85
86  z = f - 1; /* msb_index (np, n) */
87  if (h == 0 && l <= z)
88    {
89      mp_limb_t *tp2;
90      mp_size_t i;
91      int ans;
92      mp_limb_t size;
93      TMP_DECL;
94
95      size = l + k;
96      ASSERT_ALWAYS (size >= k);
97
98      TMP_MARK;
99      y = 2 + size / GMP_LIMB_BITS;
100      tp2 = TMP_ALLOC_LIMBS (y);
101
102      i = mpn_pow_1 (tp, xp, xn, k, tp2);
103      if (i == n && mpn_cmp (tp, np, n) == 0)
104	ans = 1;
105      else
106	ans = 0;
107      TMP_FREE;
108      return ans;
109    }
110
111  return 0;
112}
113
114
115/* Return non-zero if N = {np,n} is a kth power.
116   I = {ip,n} = N^(-1) mod B^n.  */
117static int
118is_kth_power (mp_ptr rp, mp_srcptr np,
119	      mp_limb_t k, mp_srcptr ip,
120	      mp_size_t n, mp_bitcnt_t f,
121	      mp_ptr tp)
122{
123  mp_bitcnt_t b;
124  mp_size_t rn, xn;
125
126  ASSERT (n > 0);
127  ASSERT ((k & 1) != 0 || k == 2);
128  ASSERT ((np[0] & 1) != 0);
129
130  if (k == 2)
131    {
132      b = (f + 1) >> 1;
133      rn = 1 + b / GMP_LIMB_BITS;
134      if (mpn_bsqrtinv (rp, ip, b, tp) != 0)
135	{
136	  rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
137	  xn = rn;
138	  MPN_NORMALIZE (rp, xn);
139	  if (pow_equals (np, n, rp, xn, k, f, tp) != 0)
140	    return 1;
141
142	  /* Check if (2^b - r)^2 == n */
143	  mpn_neg (rp, rp, rn);
144	  rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
145	  MPN_NORMALIZE (rp, rn);
146	  if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
147	    return 1;
148	}
149    }
150  else
151    {
152      b = 1 + (f - 1) / k;
153      rn = 1 + (b - 1) / GMP_LIMB_BITS;
154      mpn_brootinv (rp, ip, rn, k, tp);
155      if ((b % GMP_LIMB_BITS) != 0)
156	rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
157      MPN_NORMALIZE (rp, rn);
158      if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
159	return 1;
160    }
161  MPN_ZERO (rp, rn); /* Untrash rp */
162  return 0;
163}
164
165static int
166perfpow (mp_srcptr np, mp_size_t n,
167	 mp_limb_t ub, mp_limb_t g,
168	 mp_bitcnt_t f, int neg)
169{
170  mp_ptr ip, tp, rp;
171  mp_limb_t k;
172  int ans;
173  mp_bitcnt_t b;
174  gmp_primesieve_t ps;
175  TMP_DECL;
176
177  ASSERT (n > 0);
178  ASSERT ((np[0] & 1) != 0);
179  ASSERT (ub > 0);
180
181  TMP_MARK;
182  gmp_init_primesieve (&ps);
183  b = (f + 3) >> 1;
184
185  TMP_ALLOC_LIMBS_3 (ip, n, rp, n, tp, 5 * n);
186
187  MPN_ZERO (rp, n);
188
189  /* FIXME: It seems the inverse in ninv is needed only to get non-inverted
190     roots. I.e., is_kth_power computes n^{1/2} as (n^{-1})^{-1/2} and
191     similarly for nth roots. It should be more efficient to compute n^{1/2} as
192     n * n^{-1/2}, with a mullo instead of a binvert. And we can do something
193     similar for kth roots if we switch to an iteration converging to n^{1/k -
194     1}, and we can then eliminate this binvert call. */
195  mpn_binvert (ip, np, 1 + (b - 1) / GMP_LIMB_BITS, tp);
196  if (b % GMP_LIMB_BITS)
197    ip[(b - 1) / GMP_LIMB_BITS] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
198
199  if (neg)
200    gmp_nextprime (&ps);
201
202  ans = 0;
203  if (g > 0)
204    {
205      ub = MIN (ub, g + 1);
206      while ((k = gmp_nextprime (&ps)) < ub)
207	{
208	  if ((g % k) == 0)
209	    {
210	      if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
211		{
212		  ans = 1;
213		  goto ret;
214		}
215	    }
216	}
217    }
218  else
219    {
220      while ((k = gmp_nextprime (&ps)) < ub)
221	{
222	  if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
223	    {
224	      ans = 1;
225	      goto ret;
226	    }
227	}
228    }
229 ret:
230  TMP_FREE;
231  return ans;
232}
233
234static const unsigned short nrtrial[] = { 100, 500, 1000 };
235
236/* Table of (log_{p_i} 2) values, where p_i is the (nrtrial[i] + 1)'th prime
237   number.  */
238static const double logs[] =
239  { 0.1099457228193620, 0.0847016403115322, 0.0772048195144415 };
240
241int
242mpn_perfect_power_p (mp_srcptr np, mp_size_t n)
243{
244  mp_limb_t *nc, factor, g;
245  mp_limb_t exp, d;
246  mp_bitcnt_t twos, count;
247  int ans, where, neg, trial;
248  TMP_DECL;
249
250  neg = n < 0;
251  if (neg)
252    {
253      n = -n;
254    }
255
256  if (n == 0 || (n == 1 && np[0] == 1)) /* Valgrind doesn't like
257					   (n <= (np[0] == 1)) */
258    return 1;
259
260  TMP_MARK;
261
262  count = 0;
263
264  twos = mpn_scan1 (np, 0);
265  if (twos != 0)
266    {
267      mp_size_t s;
268      if (twos == 1)
269	{
270	  return 0;
271	}
272      s = twos / GMP_LIMB_BITS;
273      if (s + 1 == n && POW2_P (np[s]))
274	{
275	  return ! (neg && POW2_P (twos));
276	}
277      count = twos % GMP_LIMB_BITS;
278      n -= s;
279      np += s;
280      if (count > 0)
281	{
282	  nc = TMP_ALLOC_LIMBS (n);
283	  mpn_rshift (nc, np, n, count);
284	  n -= (nc[n - 1] == 0);
285	  np = nc;
286	}
287    }
288  g = twos;
289
290  trial = (n > SMALL) + (n > MEDIUM);
291
292  where = 0;
293  factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
294
295  if (factor != 0)
296    {
297      if (count == 0) /* We did not allocate nc yet. */
298	{
299	  nc = TMP_ALLOC_LIMBS (n);
300	}
301
302      /* Remove factors found by trialdiv.  Optimization: If remove
303	 define _itch, we can allocate its scratch just once */
304
305      do
306	{
307	  binvert_limb (d, factor);
308
309	  /* After the first round we always have nc == np */
310	  exp = mpn_remove (nc, &n, np, n, &d, 1, ~(mp_bitcnt_t)0);
311
312	  if (g == 0)
313	    g = exp;
314	  else
315	    g = mpn_gcd_1 (&g, 1, exp);
316
317	  if (g == 1)
318	    {
319	      ans = 0;
320	      goto ret;
321	    }
322
323	  if ((n == 1) & (nc[0] == 1))
324	    {
325	      ans = ! (neg && POW2_P (g));
326	      goto ret;
327	    }
328
329	  np = nc;
330	  factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
331	}
332      while (factor != 0);
333    }
334
335  MPN_SIZEINBASE_2EXP(count, np, n, 1);   /* log (np) + 1 */
336  d = (mp_limb_t) (count * logs[trial] + 1e-9) + 1;
337  ans = perfpow (np, n, d, g, count, neg);
338
339 ret:
340  TMP_FREE;
341  return ans;
342}
343