#
source:
trunk/third/gmp/mpz/perfpow.c
@
18191

Revision 18191, 6.1 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

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 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 | |

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; |

74 | mp_size_t usize = SIZ (u); |

75 | TMP_DECL (marker); |

76 | |

77 | if (usize == 0) |

78 | return 1; /* consider 0 a perfect power */ |

79 | |

80 | n2 = mpz_scan1 (u, 0); |

81 | if (n2 == 1) |

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 */ |

86 | |

87 | TMP_MARK (marker); |

88 | |

89 | uns = ABS (usize) - n2 / BITS_PER_MP_LIMB; |

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); |

102 | if (rem == 0) /* divisable by this prime? */ |

103 | { |

104 | rem = mpz_tdiv_q_ui (q, u2, prime * prime); |

105 | if (rem != 0) |

106 | { |

107 | TMP_FREE (marker); |

108 | return 0; /* prime divides exactly once, reject */ |

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 | |

120 | if ((n & 1) == 0 && usize < 0) |

121 | { |

122 | TMP_FREE (marker); |

123 | return 0; /* even multiplicity with negative U, reject */ |

124 | } |

125 | |

126 | n2 = gcd (n2, n); |

127 | if (n2 == 1) |

128 | { |

129 | TMP_FREE (marker); |

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 */ |

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 | { |

148 | /* We found no factors above; have to check all values of n. */ |

149 | unsigned long int nth; |

150 | for (nth = usize < 0 ? 3 : 2;; nth++) |

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 | } |

**Note:**See TracBrowser for help on using the repository browser.