Annotation of OpenXM_contrib/gmp/mpz/legendre.c, Revision 1.1
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
! 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.
! 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 Library General Public
! 16: License for more details.
! 17:
! 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. */
! 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