/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 | } |