Annotation of OpenXM_contrib/gmp/mpz/perfpow.c, Revision 1.1
1.1 ! maekawa 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to perfpow.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz