Annotation of OpenXM_contrib/gmp/mpz/perfpow.c, Revision 1.1.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