#
source:
trunk/third/gmp/tests/refmpz.c
@
18191

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

Line | |
---|---|

1 | /* Reference mpz functions. |

2 | |

3 | Copyright 1997, 1999, 2000, 2001 Free Software Foundation, Inc. |

4 | |

5 | This file is part of the GNU MP Library. |

6 | |

7 | The GNU MP Library is free software; you can redistribute it and/or modify |

8 | it under the terms of the GNU Lesser General Public License as published by |

9 | the Free Software Foundation; either version 2.1 of the License, or (at your |

10 | option) any later version. |

11 | |

12 | The GNU MP Library is distributed in the hope that it will be useful, but |

13 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

14 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

15 | License for more details. |

16 | |

17 | You should have received a copy of the GNU Lesser General Public License |

18 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

19 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

20 | MA 02111-1307, USA. */ |

21 | |

22 | /* always do assertion checking */ |

23 | #define WANT_ASSERT 1 |

24 | |

25 | #include <stdio.h> |

26 | #include <stdlib.h> /* for free */ |

27 | #include "gmp.h" |

28 | #include "gmp-impl.h" |

29 | #include "longlong.h" |

30 | #include "tests.h" |

31 | |

32 | |

33 | /* Change this to "#define TRACE(x) x" for some traces. */ |

34 | #define TRACE(x) |

35 | |

36 | |

37 | unsigned long |

38 | refmpz_hamdist (mpz_srcptr x, mpz_srcptr y) |

39 | { |

40 | mp_size_t tsize; |

41 | mp_ptr xp, yp; |

42 | unsigned long ret; |

43 | |

44 | if ((SIZ(x) < 0 && SIZ(y) >= 0) |

45 | || (SIZ(y) < 0 && SIZ(x) >= 0)) |

46 | return ULONG_MAX; |

47 | |

48 | tsize = MAX (ABSIZ(x), ABSIZ(y)); |

49 | |

50 | xp = refmpn_malloc_limbs (tsize); |

51 | refmpn_zero (xp, tsize); |

52 | refmpn_copy (xp, PTR(x), ABSIZ(x)); |

53 | |

54 | yp = refmpn_malloc_limbs (tsize); |

55 | refmpn_zero (yp, tsize); |

56 | refmpn_copy (yp, PTR(y), ABSIZ(y)); |

57 | |

58 | if (SIZ(x) < 0) |

59 | refmpn_neg_n (xp, xp, tsize); |

60 | |

61 | if (SIZ(x) < 0) |

62 | refmpn_neg_n (yp, yp, tsize); |

63 | |

64 | ret = refmpn_hamdist (xp, yp, tsize); |

65 | |

66 | free (xp); |

67 | free (yp); |

68 | return ret; |

69 | } |

70 | |

71 | |

72 | /* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */ |

73 | #define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b)) |

74 | |

75 | /* (a/b) effect due to sign of b: mpz/mpz */ |

76 | #define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b)) |

77 | |

78 | /* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd; |

79 | is (-1/b) if a<0, or +1 if a>=0 */ |

80 | #define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0]) |

81 | |

82 | int |

83 | refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig) |

84 | { |

85 | unsigned long twos; |

86 | mpz_t a, b; |

87 | int result_bit1 = 0; |

88 | |

89 | if (mpz_sgn (b_orig) == 0) |

90 | return JACOBI_Z0 (a_orig); /* (a/0) */ |

91 | |

92 | if (mpz_sgn (a_orig) == 0) |

93 | return JACOBI_0Z (b_orig); /* (0/b) */ |

94 | |

95 | if (mpz_even_p (a_orig) && mpz_even_p (b_orig)) |

96 | return 0; |

97 | |

98 | if (mpz_cmp_ui (b_orig, 1) == 0) |

99 | return 1; |

100 | |

101 | mpz_init_set (a, a_orig); |

102 | mpz_init_set (b, b_orig); |

103 | |

104 | if (mpz_sgn (b) < 0) |

105 | { |

106 | result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b); |

107 | mpz_neg (b, b); |

108 | } |

109 | if (mpz_even_p (b)) |

110 | { |

111 | twos = mpz_scan1 (b, 0L); |

112 | mpz_tdiv_q_2exp (b, b, twos); |

113 | result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); |

114 | } |

115 | |

116 | if (mpz_sgn (a) < 0) |

117 | { |

118 | result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]); |

119 | mpz_neg (a, a); |

120 | } |

121 | if (mpz_even_p (a)) |

122 | { |

123 | twos = mpz_scan1 (a, 0L); |

124 | mpz_tdiv_q_2exp (a, a, twos); |

125 | result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); |

126 | } |

127 | |

128 | for (;;) |

129 | { |

130 | ASSERT (mpz_odd_p (a)); |

131 | ASSERT (mpz_odd_p (b)); |

132 | ASSERT (mpz_sgn (a) > 0); |

133 | ASSERT (mpz_sgn (b) > 0); |

134 | |

135 | TRACE (printf ("top\n"); |

136 | mpz_trace (" a", a); |

137 | mpz_trace (" b", b)); |

138 | |

139 | if (mpz_cmp (a, b) < 0) |

140 | { |

141 | TRACE (printf ("swap\n")); |

142 | mpz_swap (a, b); |

143 | result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]); |

144 | } |

145 | |

146 | if (mpz_cmp_ui (b, 1) == 0) |

147 | break; |

148 | |

149 | mpz_sub (a, a, b); |

150 | TRACE (printf ("sub\n"); |

151 | mpz_trace (" a", a)); |

152 | if (mpz_sgn (a) == 0) |

153 | goto zero; |

154 | |

155 | twos = mpz_scan1 (a, 0L); |

156 | mpz_fdiv_q_2exp (a, a, twos); |

157 | TRACE (printf ("twos %lu\n", twos); |

158 | mpz_trace (" a", a)); |

159 | result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); |

160 | } |

161 | |

162 | mpz_clear (a); |

163 | mpz_clear (b); |

164 | return JACOBI_BIT1_TO_PN (result_bit1); |

165 | |

166 | zero: |

167 | mpz_clear (a); |

168 | mpz_clear (b); |

169 | return 0; |

170 | } |

171 | |

172 | /* Same as mpz_kronecker, but ignoring factors of 2 on b */ |

173 | int |

174 | refmpz_jacobi (mpz_srcptr a, mpz_srcptr b) |

175 | { |

176 | mpz_t b_odd; |

177 | mpz_init_set (b_odd, b); |

178 | if (mpz_sgn (b_odd) != 0) |

179 | mpz_fdiv_q_2exp (b_odd, b_odd, mpz_scan1 (b_odd, 0L)); |

180 | return refmpz_kronecker (a, b_odd); |

181 | } |

182 | |

183 | int |

184 | refmpz_legendre (mpz_srcptr a, mpz_srcptr b) |

185 | { |

186 | return refmpz_jacobi (a, b); |

187 | } |

188 | |

189 | |

190 | int |

191 | refmpz_kronecker_ui (mpz_srcptr a, unsigned long b) |

192 | { |

193 | mpz_t bz; |

194 | int ret; |

195 | mpz_init_set_ui (bz, b); |

196 | ret = refmpz_kronecker (a, bz); |

197 | mpz_clear (bz); |

198 | return ret; |

199 | } |

200 | |

201 | int |

202 | refmpz_kronecker_si (mpz_srcptr a, long b) |

203 | { |

204 | mpz_t bz; |

205 | int ret; |

206 | mpz_init_set_si (bz, b); |

207 | ret = refmpz_kronecker (a, bz); |

208 | mpz_clear (bz); |

209 | return ret; |

210 | } |

211 | |

212 | int |

213 | refmpz_ui_kronecker (unsigned long a, mpz_srcptr b) |

214 | { |

215 | mpz_t az; |

216 | int ret; |

217 | mpz_init_set_ui (az, a); |

218 | ret = refmpz_kronecker (az, b); |

219 | mpz_clear (az); |

220 | return ret; |

221 | } |

222 | |

223 | int |

224 | refmpz_si_kronecker (long a, mpz_srcptr b) |

225 | { |

226 | mpz_t az; |

227 | int ret; |

228 | mpz_init_set_si (az, a); |

229 | ret = refmpz_kronecker (az, b); |

230 | mpz_clear (az); |

231 | return ret; |

232 | } |

233 | |

234 | |

235 | void |

236 | refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e) |

237 | { |

238 | mpz_t s, t; |

239 | unsigned long i; |

240 | |

241 | mpz_init_set_ui (t, 1L); |

242 | mpz_init_set (s, b); |

243 | |

244 | if ((e & 1) != 0) |

245 | mpz_mul (t, t, s); |

246 | |

247 | for (i = 2; i <= e; i <<= 1) |

248 | { |

249 | mpz_mul (s, s, s); |

250 | if ((i & e) != 0) |

251 | mpz_mul (t, t, s); |

252 | } |

253 | |

254 | mpz_set (w, t); |

255 | |

256 | mpz_clear (s); |

257 | mpz_clear (t); |

258 | } |

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