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