refmpz.c revision 1.1.1.2
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 test suite. 6 7The GNU MP Library test suite is free software; you can redistribute it 8and/or modify it under the terms of the GNU General Public License as 9published by the Free Software Foundation; either version 3 of the License, 10or (at your option) any later version. 11 12The GNU MP Library test suite is distributed in the hope that it will be 13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15Public License for more details. 16 17You should have received a copy of the GNU General Public License along with 18the GNU MP Library test suite. 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 ASSERT_ALWAYS (mpz_sgn (b) > 0); 188 ASSERT_ALWAYS (mpz_odd_p (b)); 189 190 return refmpz_kronecker (a, b); 191} 192 193/* Legendre symbol via powm. p must be an odd prime. */ 194int 195refmpz_legendre (mpz_srcptr a, mpz_srcptr p) 196{ 197 int res; 198 199 mpz_t r; 200 mpz_t e; 201 202 ASSERT_ALWAYS (mpz_sgn (p) > 0); 203 ASSERT_ALWAYS (mpz_odd_p (p)); 204 205 mpz_init (r); 206 mpz_init (e); 207 208 mpz_fdiv_r (r, a, p); 209 210 mpz_set (e, p); 211 mpz_sub_ui (e, e, 1); 212 mpz_fdiv_q_2exp (e, e, 1); 213 mpz_powm (r, r, e, p); 214 215 /* Normalize to a more or less symmetric range around zero */ 216 if (mpz_cmp (r, e) > 0) 217 mpz_sub (r, r, p); 218 219 ASSERT_ALWAYS (mpz_cmpabs_ui (r, 1) <= 0); 220 221 res = mpz_sgn (r); 222 223 mpz_clear (r); 224 mpz_clear (e); 225 226 return res; 227} 228 229 230int 231refmpz_kronecker_ui (mpz_srcptr a, unsigned long b) 232{ 233 mpz_t bz; 234 int ret; 235 mpz_init_set_ui (bz, b); 236 ret = refmpz_kronecker (a, bz); 237 mpz_clear (bz); 238 return ret; 239} 240 241int 242refmpz_kronecker_si (mpz_srcptr a, long b) 243{ 244 mpz_t bz; 245 int ret; 246 mpz_init_set_si (bz, b); 247 ret = refmpz_kronecker (a, bz); 248 mpz_clear (bz); 249 return ret; 250} 251 252int 253refmpz_ui_kronecker (unsigned long a, mpz_srcptr b) 254{ 255 mpz_t az; 256 int ret; 257 mpz_init_set_ui (az, a); 258 ret = refmpz_kronecker (az, b); 259 mpz_clear (az); 260 return ret; 261} 262 263int 264refmpz_si_kronecker (long a, mpz_srcptr b) 265{ 266 mpz_t az; 267 int ret; 268 mpz_init_set_si (az, a); 269 ret = refmpz_kronecker (az, b); 270 mpz_clear (az); 271 return ret; 272} 273 274 275void 276refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e) 277{ 278 mpz_t s, t; 279 unsigned long i; 280 281 mpz_init_set_ui (t, 1L); 282 mpz_init_set (s, b); 283 284 if ((e & 1) != 0) 285 mpz_mul (t, t, s); 286 287 for (i = 2; i <= e; i <<= 1) 288 { 289 mpz_mul (s, s, s); 290 if ((i & e) != 0) 291 mpz_mul (t, t, s); 292 } 293 294 mpz_set (w, t); 295 296 mpz_clear (s); 297 mpz_clear (t); 298} 299