Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/perfpow.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_perfect_power_p -- mpn perfect power detection.
2
3
   Contributed to the GNU project by Martin Boij.
4
5
Copyright 2009, 2010, 2012, 2014 Free Software Foundation, Inc.
6
7
This file is part of the GNU MP Library.
8
9
The GNU MP Library is free software; you can redistribute it and/or modify
10
it 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
16
or
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
22
or both in parallel, as here.
23
24
The GNU MP Library is distributed in the hope that it will be useful, but
25
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
27
for more details.
28
29
You should have received copies of the GNU General Public License and the
30
GNU Lesser General Public License along with the GNU MP Library.  If not,
31
see https://www.gnu.org/licenses/.  */
32
33
#include "gmp-impl.h"
34
#include "longlong.h"
35
36
203
#define SMALL 20
37
203
#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.  */
50
static int
51
pow_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
1.11k
{
56
1.11k
  mp_bitcnt_t y, z;
57
1.11k
  mp_size_t bn;
58
1.11k
  mp_limb_t h, l;
59
60
1.11k
  ASSERT (n > 1 || (n == 1 && np[0] > 1));
61
1.11k
  ASSERT (np[n - 1] > 0);
62
1.11k
  ASSERT (xn > 0);
63
64
1.11k
  if (xn == 1 && xp[0] == 1)
65
71
    return 0;
66
67
1.04k
  z = 1 + (n >> 1);
68
1.86k
  for (bn = 1; bn < z; bn <<= 1)
69
1.73k
    {
70
1.73k
      mpn_powlo (tp, xp, &k, 1, bn, tp + bn);
71
1.73k
      if (mpn_cmp (tp, np, bn) != 0)
72
908
  return 0;
73
1.73k
    }
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
134
  MPN_SIZEINBASE_2EXP(y, xp, xn, 1);
81
134
  y -= 1;  /* msb_index (xp, xn) */
82
83
134
  umul_ppmm (h, l, k, y);
84
134
  h -= l == 0;  --l;  /* two-limb decrement */
85
86
134
  z = f - 1; /* msb_index (np, n) */
87
134
  if (h == 0 && l <= z)
88
134
    {
89
134
      mp_limb_t *tp2;
90
134
      mp_size_t i;
91
134
      int ans;
92
134
      mp_limb_t size;
93
134
      TMP_DECL;
94
95
134
      size = l + k;
96
134
      ASSERT_ALWAYS (size >= k);
97
98
134
      TMP_MARK;
99
134
      y = 2 + size / GMP_LIMB_BITS;
100
134
      tp2 = TMP_ALLOC_LIMBS (y);
101
102
134
      i = mpn_pow_1 (tp, xp, xn, k, tp2);
103
134
      if (i == n && mpn_cmp (tp, np, n) == 0)
104
0
  ans = 1;
105
134
      else
106
134
  ans = 0;
107
134
      TMP_FREE;
108
134
      return ans;
109
134
    }
110
111
0
  return 0;
112
134
}
113
114
115
/* Return non-zero if N = {np,n} is a kth power.
116
   I = {ip,n} = N^(-1) mod B^n.  */
117
static int
118
is_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
1.06k
{
123
1.06k
  mp_bitcnt_t b;
124
1.06k
  mp_size_t rn, xn;
125
126
1.06k
  ASSERT (n > 0);
127
1.06k
  ASSERT ((k & 1) != 0 || k == 2);
128
1.06k
  ASSERT ((np[0] & 1) != 0);
129
130
1.06k
  if (k == 2)
131
85
    {
132
85
      b = (f + 1) >> 1;
133
85
      rn = 1 + b / GMP_LIMB_BITS;
134
85
      if (mpn_bsqrtinv (rp, ip, b, tp) != 0)
135
65
  {
136
65
    rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
137
65
    xn = rn;
138
65
    MPN_NORMALIZE (rp, xn);
139
65
    if (pow_equals (np, n, rp, xn, k, f, tp) != 0)
140
0
      return 1;
141
142
    /* Check if (2^b - r)^2 == n */
143
65
    mpn_neg (rp, rp, rn);
144
65
    rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
145
65
    MPN_NORMALIZE (rp, rn);
146
65
    if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
147
0
      return 1;
148
65
  }
149
85
    }
150
983
  else
151
983
    {
152
983
      b = 1 + (f - 1) / k;
153
983
      rn = 1 + (b - 1) / GMP_LIMB_BITS;
154
983
      mpn_brootinv (rp, ip, rn, k, tp);
155
983
      if ((b % GMP_LIMB_BITS) != 0)
156
975
  rp[rn - 1] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
157
983
      MPN_NORMALIZE (rp, rn);
158
983
      if (pow_equals (np, n, rp, rn, k, f, tp) != 0)
159
0
  return 1;
160
983
    }
161
1.06k
  MPN_ZERO (rp, rn); /* Untrash rp */
162
1.06k
  return 0;
163
1.06k
}
164
165
static int
166
perfpow (mp_srcptr np, mp_size_t n,
167
   mp_limb_t ub, mp_limb_t g,
168
   mp_bitcnt_t f, int neg)
169
119
{
170
119
  mp_ptr ip, tp, rp;
171
119
  mp_limb_t k;
172
119
  int ans;
173
119
  mp_bitcnt_t b;
174
119
  gmp_primesieve_t ps;
175
119
  TMP_DECL;
176
177
119
  ASSERT (n > 0);
178
119
  ASSERT ((np[0] & 1) != 0);
179
119
  ASSERT (ub > 0);
180
181
119
  TMP_MARK;
182
119
  gmp_init_primesieve (&ps);
183
119
  b = (f + 3) >> 1;
184
185
119
  TMP_ALLOC_LIMBS_3 (ip, n, rp, n, tp, 5 * n);
186
187
119
  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
119
  mpn_binvert (ip, np, 1 + (b - 1) / GMP_LIMB_BITS, tp);
196
119
  if (b % GMP_LIMB_BITS)
197
114
    ip[(b - 1) / GMP_LIMB_BITS] &= (CNST_LIMB(1) << (b % GMP_LIMB_BITS)) - 1;
198
199
119
  if (neg)
200
0
    gmp_nextprime (&ps);
201
202
119
  ans = 0;
203
119
  if (g > 0)
204
86
    {
205
86
      ub = MIN (ub, g + 1);
206
687
      while ((k = gmp_nextprime (&ps)) < ub)
207
601
  {
208
601
    if ((g % k) == 0)
209
112
      {
210
112
        if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
211
0
    {
212
0
      ans = 1;
213
0
      goto ret;
214
0
    }
215
112
      }
216
601
  }
217
86
    }
218
33
  else
219
33
    {
220
989
      while ((k = gmp_nextprime (&ps)) < ub)
221
956
  {
222
956
    if (is_kth_power (rp, np, k, ip, n, f, tp) != 0)
223
0
      {
224
0
        ans = 1;
225
0
        goto ret;
226
0
      }
227
956
  }
228
33
    }
229
119
 ret:
230
119
  TMP_FREE;
231
119
  return ans;
232
119
}
233
234
static 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.  */
238
static const double logs[] =
239
  { 0.1099457228193620, 0.0847016403115322, 0.0772048195144415 };
240
241
int
242
mpn_perfect_power_p (mp_srcptr np, mp_size_t n)
243
213
{
244
213
  mp_limb_t *nc, factor, g;
245
213
  mp_limb_t exp, d;
246
213
  mp_bitcnt_t twos, count;
247
213
  int ans, where, neg, trial;
248
213
  TMP_DECL;
249
250
213
  neg = n < 0;
251
213
  if (neg)
252
0
    {
253
0
      n = -n;
254
0
    }
255
256
213
  if (n == 0 || (n == 1 && np[0] == 1)) /* Valgrind doesn't like
257
             (n <= (np[0] == 1)) */
258
5
    return 1;
259
260
208
  TMP_MARK;
261
262
208
  count = 0;
263
264
208
  twos = mpn_scan1 (np, 0);
265
208
  if (twos != 0)
266
153
    {
267
153
      mp_size_t s;
268
153
      if (twos == 1)
269
5
  {
270
5
    return 0;
271
5
  }
272
148
      s = twos / GMP_LIMB_BITS;
273
148
      if (s + 1 == n && POW2_P (np[s]))
274
0
  {
275
0
    return ! (neg && POW2_P (twos));
276
0
  }
277
148
      count = twos % GMP_LIMB_BITS;
278
148
      n -= s;
279
148
      np += s;
280
148
      if (count > 0)
281
147
  {
282
147
    nc = TMP_ALLOC_LIMBS (n);
283
147
    mpn_rshift (nc, np, n, count);
284
147
    n -= (nc[n - 1] == 0);
285
147
    np = nc;
286
147
  }
287
148
    }
288
203
  g = twos;
289
290
203
  trial = (n > SMALL) + (n > MEDIUM);
291
292
203
  where = 0;
293
203
  factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
294
295
203
  if (factor != 0)
296
166
    {
297
166
      if (count == 0) /* We did not allocate nc yet. */
298
23
  {
299
23
    nc = TMP_ALLOC_LIMBS (n);
300
23
  }
301
302
      /* Remove factors found by trialdiv.  Optimization: If remove
303
   define _itch, we can allocate its scratch just once */
304
305
166
      do
306
240
  {
307
240
    binvert_limb (d, factor);
308
309
    /* After the first round we always have nc == np */
310
240
    exp = mpn_remove (nc, &n, np, n, &d, 1, ~(mp_bitcnt_t)0);
311
312
240
    if (g == 0)
313
22
      g = exp;
314
218
    else
315
218
      g = mpn_gcd_1 (&g, 1, exp);
316
317
240
    if (g == 1)
318
81
      {
319
81
        ans = 0;
320
81
        goto ret;
321
81
      }
322
323
159
    if ((n == 1) & (nc[0] == 1))
324
3
      {
325
3
        ans = ! (neg && POW2_P (g));
326
3
        goto ret;
327
3
      }
328
329
156
    np = nc;
330
156
    factor = mpn_trialdiv (np, n, nrtrial[trial], &where);
331
156
  }
332
166
      while (factor != 0);
333
166
    }
334
335
119
  MPN_SIZEINBASE_2EXP(count, np, n, 1);   /* log (np) + 1 */
336
119
  d = (mp_limb_t) (count * logs[trial] + 1e-9) + 1;
337
119
  ans = perfpow (np, n, d, g, count, neg);
338
339
203
 ret:
340
203
  TMP_FREE;
341
203
  return ans;
342
119
}