Annotation of OpenXM_contrib/gmp/tests/refmpz.c, Revision 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