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