1 /* mpz_perfect_power_p(arg) -- Return non-zero if ARG is a perfect power, 2 zero otherwise. 3 4 Copyright 1998, 1999, 2000, 2001, 2005, 2008 Free Software Foundation, Inc. 5 6 This file is part of the GNU MP Library. 7 8 The GNU MP Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 20 21 /* 22 We are to determine if c is a perfect power, c = a ^ b. 23 Assume c is divisible by 2^n and that codd = c/2^n is odd. 24 Assume a is divisible by 2^m and that aodd = a/2^m is odd. 25 It is always true that m divides n. 26 27 * If n is prime, either 1) a is 2*aodd and b = n 28 or 2) a = c and b = 1. 29 So for n prime, we readily have a solution. 30 * If n is factorable into the non-trivial factors p1,p2,... 31 Since m divides n, m has a subset of n's factors and b = n / m. 32 */ 33 34 /* This is a naive approach to recognizing perfect powers. 35 Many things can be improved. In particular, we should use p-adic 36 arithmetic for computing possible roots. */ 37 38 #include <stdio.h> /* for NULL */ 39 #include "gmp.h" 40 #include "gmp-impl.h" 41 #include "longlong.h" 42 43 static unsigned long int gcd __GMP_PROTO ((unsigned long int, unsigned long int)); 44 static int isprime __GMP_PROTO ((unsigned long int)); 45 46 static const unsigned short primes[] = 47 { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 48 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131, 49 137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223, 50 227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311, 51 313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409, 52 419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503, 53 509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613, 54 617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719, 55 727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827, 56 829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941, 57 947,953,967,971,977,983,991,997,0 58 }; 59 #define SMALLEST_OMITTED_PRIME 1009 60 61 62 #define POW2P(a) (((a) & ((a) - 1)) == 0) 63 64 int 65 mpz_perfect_power_p (mpz_srcptr u) 66 { 67 unsigned long int prime; 68 unsigned long int n, n2; 69 int i; 70 unsigned long int rem; 71 mpz_t u2, q; 72 int exact; 73 mp_size_t uns; 74 mp_size_t usize = SIZ (u); 75 TMP_DECL; 76 77 if (mpz_cmpabs_ui (u, 1) <= 0) 78 return 1; /* -1, 0, and +1 are perfect powers */ 79 80 n2 = mpz_scan1 (u, 0); 81 if (n2 == 1) 82 return 0; /* 2 divides exactly once. */ 83 84 TMP_MARK; 85 86 uns = ABS (usize) - n2 / BITS_PER_MP_LIMB; 87 MPZ_TMP_INIT (q, uns); 88 MPZ_TMP_INIT (u2, uns); 89 90 mpz_tdiv_q_2exp (u2, u, n2); 91 mpz_abs (u2, u2); 92 93 if (mpz_cmp_ui (u2, 1) == 0) 94 { 95 TMP_FREE; 96 /* factoring completed; consistent power */ 97 return ! (usize < 0 && POW2P(n2)); 98 } 99 100 if (isprime (n2)) 101 goto n2prime; 102 103 for (i = 1; primes[i] != 0; i++) 104 { 105 prime = primes[i]; 106 107 if (mpz_cmp_ui (u2, prime) < 0) 108 break; 109 110 if (mpz_divisible_ui_p (u2, prime)) /* divisible by this prime? */ 111 { 112 rem = mpz_tdiv_q_ui (q, u2, prime * prime); 113 if (rem != 0) 114 { 115 TMP_FREE; 116 return 0; /* prime divides exactly once, reject */ 117 } 118 mpz_swap (q, u2); 119 for (n = 2;;) 120 { 121 rem = mpz_tdiv_q_ui (q, u2, prime); 122 if (rem != 0) 123 break; 124 mpz_swap (q, u2); 125 n++; 126 } 127 128 n2 = gcd (n2, n); 129 if (n2 == 1) 130 { 131 TMP_FREE; 132 return 0; /* we have multiplicity 1 of some factor */ 133 } 134 135 if (mpz_cmp_ui (u2, 1) == 0) 136 { 137 TMP_FREE; 138 /* factoring completed; consistent power */ 139 return ! (usize < 0 && POW2P(n2)); 140 } 141 142 /* As soon as n2 becomes a prime number, stop factoring. 143 Either we have u=x^n2 or u is not a perfect power. */ 144 if (isprime (n2)) 145 goto n2prime; 146 } 147 } 148 149 if (n2 == 0) 150 { 151 /* We found no factors above; have to check all values of n. */ 152 unsigned long int nth; 153 for (nth = usize < 0 ? 3 : 2;; nth++) 154 { 155 if (! isprime (nth)) 156 continue; 157 #if 0 158 exact = mpz_padic_root (q, u2, nth, PTH); 159 if (exact) 160 #endif 161 exact = mpz_root (q, u2, nth); 162 if (exact) 163 { 164 TMP_FREE; 165 return 1; 166 } 167 if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0) 168 { 169 TMP_FREE; 170 return 0; 171 } 172 } 173 } 174 else 175 { 176 unsigned long int nth; 177 178 if (usize < 0 && POW2P(n2)) 179 { 180 TMP_FREE; 181 return 0; 182 } 183 184 /* We found some factors above. We just need to consider values of n 185 that divides n2. */ 186 for (nth = 2; nth <= n2; nth++) 187 { 188 if (! isprime (nth)) 189 continue; 190 if (n2 % nth != 0) 191 continue; 192 #if 0 193 exact = mpz_padic_root (q, u2, nth, PTH); 194 if (exact) 195 #endif 196 exact = mpz_root (q, u2, nth); 197 if (exact) 198 { 199 if (! (usize < 0 && POW2P(nth))) 200 { 201 TMP_FREE; 202 return 1; 203 } 204 } 205 if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0) 206 { 207 TMP_FREE; 208 return 0; 209 } 210 } 211 212 TMP_FREE; 213 return 0; 214 } 215 216 n2prime: 217 if (usize < 0 && POW2P(n2)) 218 { 219 TMP_FREE; 220 return 0; 221 } 222 223 exact = mpz_root (NULL, u2, n2); 224 TMP_FREE; 225 return exact; 226 } 227 228 static unsigned long int 229 gcd (unsigned long int a, unsigned long int b) 230 { 231 int an2, bn2, n2; 232 233 if (a == 0) 234 return b; 235 if (b == 0) 236 return a; 237 238 count_trailing_zeros (an2, a); 239 a >>= an2; 240 241 count_trailing_zeros (bn2, b); 242 b >>= bn2; 243 244 n2 = MIN (an2, bn2); 245 246 while (a != b) 247 { 248 if (a > b) 249 { 250 a -= b; 251 do 252 a >>= 1; 253 while ((a & 1) == 0); 254 } 255 else /* b > a. */ 256 { 257 b -= a; 258 do 259 b >>= 1; 260 while ((b & 1) == 0); 261 } 262 } 263 264 return a << n2; 265 } 266 267 static int 268 isprime (unsigned long int t) 269 { 270 unsigned long int q, r, d; 271 272 if (t < 3 || (t & 1) == 0) 273 return t == 2; 274 275 for (d = 3, r = 1; r != 0; d += 2) 276 { 277 q = t / d; 278 r = t - q * d; 279 if (q < d) 280 return 1; 281 } 282 return 0; 283 } 284