Annotation of OpenXM_contrib/gmp/mpz/legendre.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpz_legendre (op1, op2).
2: Contributed by Bennet Yee (bsy) at Carnegie-Mellon University
3:
4: Copyright (C) 1992, 1996 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
1.1.1.2 ! maekawa 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
1.1 maekawa 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
1.1.1.2 ! maekawa 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 ! maekawa 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 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: #include "gmp.h"
24:
25: #if defined (DEBUG)
26: #include <stdio.h>
27: #endif
28:
29: /* Precondition: both p and q are positive */
30:
31: int
32: #if __STDC__
33: mpz_legendre (mpz_srcptr pi, mpz_srcptr qi)
34: #else
35: mpz_legendre (pi, qi)
36: mpz_srcptr pi, qi;
37: #endif
38: {
39: mpz_t p, q, qdiv2;
40: #ifdef Q_MINUS_1
41: mpz_t q_minus_1;
42: #endif
43: mpz_ptr mtmp;
44: register mpz_ptr pptr, qptr;
45: register int retval = 1;
46: register unsigned long int s;
47:
48: pptr = p;
49: mpz_init_set (pptr, pi);
50: qptr = q;
51: mpz_init_set (qptr, qi);
52:
53: #ifdef Q_MINUS_1
54: mpz_init (q_minus_1);
55: #endif
56: mpz_init (qdiv2);
57:
58: tail_recurse2:
59: #ifdef DEBUG
60: printf ("tail_recurse2: p=");
61: mpz_out_str (stdout, 10, pptr);
62: printf ("\nq=");
63: mpz_out_str (stdout, 10, qptr);
64: putchar ('\n');
65: #endif
66: s = mpz_scan1 (qptr, 0);
67: if (s) mpz_tdiv_q_2exp (qptr, qptr, s); /* J(a,2) = 1 */
68: #ifdef DEBUG
69: printf ("2 factor decomposition: p=");
70: mpz_out_str (stdout, 10, pptr);
71: printf ("\nq=");
72: mpz_out_str (stdout, 10, qptr);
73: putchar ('\n');
74: #endif
75: /* postcondition q odd */
76: if (!mpz_cmp_ui (qptr, 1L)) /* J(a,1) = 1 */
77: goto done;
78: mpz_mod (pptr, pptr, qptr); /* J(a,q) = J(b,q) when a == b mod q */
79: #ifdef DEBUG
80: printf ("mod out by q: p=");
81: mpz_out_str (stdout, 10, pptr);
82: printf ("\nq=");
83: mpz_out_str (stdout, 10, qptr);
84: putchar ('\n');
85: #endif
86: /* quick calculation to get approximate size first */
87: /* precondition: p < q */
88: if ((mpz_sizeinbase (pptr, 2) + 1 >= mpz_sizeinbase (qptr,2))
89: && (mpz_tdiv_q_2exp (qdiv2, qptr, 1L), mpz_cmp (pptr, qdiv2) > 0))
90: {
91: /* p > q/2 */
92: mpz_sub (pptr, qptr, pptr);
93: /* J(-1,q) = (-1)^((q-1)/2), q odd */
94: if (mpz_get_ui (qptr) & 2)
95: retval = -retval;
96: }
97: /* p < q/2 */
98: #ifdef Q_MINUS_1
99: mpz_sub_ui (q_minus_q, qptr, 1L);
100: #endif
101: tail_recurse: /* we use tail_recurse only if q has not changed */
102: #ifdef DEBUG
103: printf ("tail_recurse1: p=");
104: mpz_out_str (stdout, 10, pptr);
105: printf ("\nq=");
106: mpz_out_str (stdout, 10, qptr);
107: putchar ('\n');
108: #endif
109: /*
110: * J(0,q) = 0
111: * this occurs only if gcd(p,q) != 1 which is never true for
112: * Legendre function.
113: */
114: if (!mpz_cmp_ui (pptr, 0L))
115: {
116: retval = 0;
117: goto done;
118: }
119:
120: if (!mpz_cmp_ui (pptr, 1L))
121: {
122: /* J(1,q) = 1 */
123: /* retval *= 1; */
124: goto done;
125: }
126: #ifdef Q_MINUS_1
127: if (!mpz_cmp (pptr, q_minus_1))
128: {
129: /* J(-1,q) = (-1)^((q-1)/2) */
130: if (mpz_get_ui (qptr) & 2)
131: retval = -retval;
132: /* else retval *= 1; */
133: goto done;
134: }
135: #endif
136: /*
137: * we do not handle J(xy,q) except for x==2
138: * since we do not want to factor
139: */
140: if ((s = mpz_scan1 (pptr, 0)) != 0)
141: {
142: /*
143: * J(2,q) = (-1)^((q^2-1)/8)
144: *
145: * Note that q odd guarantees that q^2-1 is divisible by 8:
146: * Let a: q=2a+1. q^2 = 4a^2+4a+1, (q^2-1)/8 = a(a+1)/2, qed
147: *
148: * Now, note that this means that the low two bits of _a_
149: * (or the low bits of q shifted over by 1 determines
150: * the factor).
151: */
152: mpz_tdiv_q_2exp (pptr, pptr, s);
153:
154: /* even powers of 2 gives J(2,q)^{2n} = 1 */
155: if (s & 1)
156: {
157: s = mpz_get_ui (qptr) >> 1;
158: s = s * (s + 1);
159: if (s & 2)
160: retval = -retval;
161: }
162: goto tail_recurse;
163: }
164: /*
165: * we know p is odd since we have cast out 2s
166: * precondition that q is odd guarantees both odd.
167: *
168: * quadratic reciprocity
169: * J(p,q) = (-1)^((p-1)(q-1)/4) * J(q,p)
170: */
171: if ((s = mpz_scan1 (pptr, 1)) <= 2 && (s + mpz_scan1 (qptr, 1)) <= 2)
172: retval = -retval;
173:
174: mtmp = pptr; pptr = qptr; qptr = mtmp;
175: goto tail_recurse2;
176: done:
177: mpz_clear (p);
178: mpz_clear (q);
179: mpz_clear (qdiv2);
180: #ifdef Q_MINUS_1
181: mpz_clear (q_minus_1);
182: #endif
183: return retval;
184: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to legendre.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz