refmpz.c revision 1.1.1.1
1/* Reference mpz functions. 2 3Copyright 1997, 1999, 2000, 2001, 2002 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library. 6 7The GNU MP Library is free software; you can redistribute it and/or modify 8it under the terms of the GNU Lesser General Public License as published by 9the Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12The GNU MP Library is distributed in the hope that it will be useful, but 13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15License for more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 19 20/* always do assertion checking */ 21#define WANT_ASSERT 1 22 23#include <stdio.h> 24#include <stdlib.h> /* for free */ 25#include "gmp.h" 26#include "gmp-impl.h" 27#include "longlong.h" 28#include "tests.h" 29 30 31/* Change this to "#define TRACE(x) x" for some traces. */ 32#define TRACE(x) 33 34 35/* FIXME: Shouldn't use plain mpz functions in a reference routine. */ 36void 37refmpz_combit (mpz_ptr r, unsigned long bit) 38{ 39 if (mpz_tstbit (r, bit)) 40 mpz_clrbit (r, bit); 41 else 42 mpz_setbit (r, bit); 43} 44 45 46unsigned long 47refmpz_hamdist (mpz_srcptr x, mpz_srcptr y) 48{ 49 mp_size_t xsize, ysize, tsize; 50 mp_ptr xp, yp; 51 unsigned long ret; 52 53 if ((SIZ(x) < 0 && SIZ(y) >= 0) 54 || (SIZ(y) < 0 && SIZ(x) >= 0)) 55 return ULONG_MAX; 56 57 xsize = ABSIZ(x); 58 ysize = ABSIZ(y); 59 tsize = MAX (xsize, ysize); 60 61 xp = refmpn_malloc_limbs (tsize); 62 refmpn_zero (xp, tsize); 63 refmpn_copy (xp, PTR(x), xsize); 64 65 yp = refmpn_malloc_limbs (tsize); 66 refmpn_zero (yp, tsize); 67 refmpn_copy (yp, PTR(y), ysize); 68 69 if (SIZ(x) < 0) 70 refmpn_neg (xp, xp, tsize); 71 72 if (SIZ(x) < 0) 73 refmpn_neg (yp, yp, tsize); 74 75 ret = refmpn_hamdist (xp, yp, tsize); 76 77 free (xp); 78 free (yp); 79 return ret; 80} 81 82 83/* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */ 84#define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b)) 85 86/* (a/b) effect due to sign of b: mpz/mpz */ 87#define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b)) 88 89/* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd; 90 is (-1/b) if a<0, or +1 if a>=0 */ 91#define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0]) 92 93int 94refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig) 95{ 96 unsigned long twos; 97 mpz_t a, b; 98 int result_bit1 = 0; 99 100 if (mpz_sgn (b_orig) == 0) 101 return JACOBI_Z0 (a_orig); /* (a/0) */ 102 103 if (mpz_sgn (a_orig) == 0) 104 return JACOBI_0Z (b_orig); /* (0/b) */ 105 106 if (mpz_even_p (a_orig) && mpz_even_p (b_orig)) 107 return 0; 108 109 if (mpz_cmp_ui (b_orig, 1) == 0) 110 return 1; 111 112 mpz_init_set (a, a_orig); 113 mpz_init_set (b, b_orig); 114 115 if (mpz_sgn (b) < 0) 116 { 117 result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b); 118 mpz_neg (b, b); 119 } 120 if (mpz_even_p (b)) 121 { 122 twos = mpz_scan1 (b, 0L); 123 mpz_tdiv_q_2exp (b, b, twos); 124 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); 125 } 126 127 if (mpz_sgn (a) < 0) 128 { 129 result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]); 130 mpz_neg (a, a); 131 } 132 if (mpz_even_p (a)) 133 { 134 twos = mpz_scan1 (a, 0L); 135 mpz_tdiv_q_2exp (a, a, twos); 136 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); 137 } 138 139 for (;;) 140 { 141 ASSERT (mpz_odd_p (a)); 142 ASSERT (mpz_odd_p (b)); 143 ASSERT (mpz_sgn (a) > 0); 144 ASSERT (mpz_sgn (b) > 0); 145 146 TRACE (printf ("top\n"); 147 mpz_trace (" a", a); 148 mpz_trace (" b", b)); 149 150 if (mpz_cmp (a, b) < 0) 151 { 152 TRACE (printf ("swap\n")); 153 mpz_swap (a, b); 154 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]); 155 } 156 157 if (mpz_cmp_ui (b, 1) == 0) 158 break; 159 160 mpz_sub (a, a, b); 161 TRACE (printf ("sub\n"); 162 mpz_trace (" a", a)); 163 if (mpz_sgn (a) == 0) 164 goto zero; 165 166 twos = mpz_scan1 (a, 0L); 167 mpz_fdiv_q_2exp (a, a, twos); 168 TRACE (printf ("twos %lu\n", twos); 169 mpz_trace (" a", a)); 170 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); 171 } 172 173 mpz_clear (a); 174 mpz_clear (b); 175 return JACOBI_BIT1_TO_PN (result_bit1); 176 177 zero: 178 mpz_clear (a); 179 mpz_clear (b); 180 return 0; 181} 182 183/* Same as mpz_kronecker, but ignoring factors of 2 on b */ 184int 185refmpz_jacobi (mpz_srcptr a, mpz_srcptr b) 186{ 187 mpz_t b_odd; 188 mpz_init_set (b_odd, b); 189 if (mpz_sgn (b_odd) != 0) 190 mpz_fdiv_q_2exp (b_odd, b_odd, mpz_scan1 (b_odd, 0L)); 191 return refmpz_kronecker (a, b_odd); 192} 193 194int 195refmpz_legendre (mpz_srcptr a, mpz_srcptr b) 196{ 197 return refmpz_jacobi (a, b); 198} 199 200 201int 202refmpz_kronecker_ui (mpz_srcptr a, unsigned long b) 203{ 204 mpz_t bz; 205 int ret; 206 mpz_init_set_ui (bz, b); 207 ret = refmpz_kronecker (a, bz); 208 mpz_clear (bz); 209 return ret; 210} 211 212int 213refmpz_kronecker_si (mpz_srcptr a, long b) 214{ 215 mpz_t bz; 216 int ret; 217 mpz_init_set_si (bz, b); 218 ret = refmpz_kronecker (a, bz); 219 mpz_clear (bz); 220 return ret; 221} 222 223int 224refmpz_ui_kronecker (unsigned long a, mpz_srcptr b) 225{ 226 mpz_t az; 227 int ret; 228 mpz_init_set_ui (az, a); 229 ret = refmpz_kronecker (az, b); 230 mpz_clear (az); 231 return ret; 232} 233 234int 235refmpz_si_kronecker (long a, mpz_srcptr b) 236{ 237 mpz_t az; 238 int ret; 239 mpz_init_set_si (az, a); 240 ret = refmpz_kronecker (az, b); 241 mpz_clear (az); 242 return ret; 243} 244 245 246void 247refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e) 248{ 249 mpz_t s, t; 250 unsigned long i; 251 252 mpz_init_set_ui (t, 1L); 253 mpz_init_set (s, b); 254 255 if ((e & 1) != 0) 256 mpz_mul (t, t, s); 257 258 for (i = 2; i <= e; i <<= 1) 259 { 260 mpz_mul (s, s, s); 261 if ((i & e) != 0) 262 mpz_mul (t, t, s); 263 } 264 265 mpz_set (w, t); 266 267 mpz_clear (s); 268 mpz_clear (t); 269} 270