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