contrib/gmp/mpz/perfpow.c

   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 }