1 /* mpz_legendre (op1, op2).
2 Contributed by Bennet Yee (bsy) at Carnegie-Mellon University
4 Copyright (C) 1992, 1996 Free Software Foundation, Inc.
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Library General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or (at your
11 option) any later version.
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 Library General Public
16 License for more details.
18 You should have received a copy of the GNU Library 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. */
29 /* Precondition: both p and q are positive */
33 mpz_legendre (mpz_srcptr pi, mpz_srcptr qi)
44 register mpz_ptr pptr, qptr;
45 register int retval = 1;
46 register unsigned long int s;
49 mpz_init_set (pptr, pi);
51 mpz_init_set (qptr, qi);
60 printf ("tail_recurse2: p=");
61 mpz_out_str (stdout, 10, pptr);
63 mpz_out_str (stdout, 10, qptr);
66 s = mpz_scan1 (qptr, 0);
67 if (s) mpz_tdiv_q_2exp (qptr, qptr, s); /* J(a,2) = 1 */
69 printf ("2 factor decomposition: p=");
70 mpz_out_str (stdout, 10, pptr);
72 mpz_out_str (stdout, 10, qptr);
75 /* postcondition q odd */
76 if (!mpz_cmp_ui (qptr, 1L)) /* J(a,1) = 1 */
78 mpz_mod (pptr, pptr, qptr); /* J(a,q) = J(b,q) when a == b mod q */
80 printf ("mod out by q: p=");
81 mpz_out_str (stdout, 10, pptr);
83 mpz_out_str (stdout, 10, qptr);
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))
92 mpz_sub (pptr, qptr, pptr);
93 /* J(-1,q) = (-1)^((q-1)/2), q odd */
94 if (mpz_get_ui (qptr) & 2)
99 mpz_sub_ui (q_minus_q, qptr, 1L);
101 tail_recurse: /* we use tail_recurse only if q has not changed */
103 printf ("tail_recurse1: p=");
104 mpz_out_str (stdout, 10, pptr);
106 mpz_out_str (stdout, 10, qptr);
111 * this occurs only if gcd(p,q) != 1 which is never true for
114 if (!mpz_cmp_ui (pptr, 0L))
120 if (!mpz_cmp_ui (pptr, 1L))
127 if (!mpz_cmp (pptr, q_minus_1))
129 /* J(-1,q) = (-1)^((q-1)/2) */
130 if (mpz_get_ui (qptr) & 2)
132 /* else retval *= 1; */
137 * we do not handle J(xy,q) except for x==2
138 * since we do not want to factor
140 if ((s = mpz_scan1 (pptr, 0)) != 0)
143 * J(2,q) = (-1)^((q^2-1)/8)
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
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
152 mpz_tdiv_q_2exp (pptr, pptr, s);
154 /* even powers of 2 gives J(2,q)^{2n} = 1 */
157 s = mpz_get_ui (qptr) >> 1;
165 * we know p is odd since we have cast out 2s
166 * precondition that q is odd guarantees both odd.
168 * quadratic reciprocity
169 * J(p,q) = (-1)^((p-1)(q-1)/4) * J(q,p)
171 if ((s = mpz_scan1 (pptr, 1)) <= 2 && (s + mpz_scan1 (qptr, 1)) <= 2)
174 mtmp = pptr; pptr = qptr; qptr = mtmp;
181 mpz_clear (q_minus_1);