1 /* Reference mpz functions. 2 3 Copyright 1997, 1999-2002 Free Software Foundation, Inc. 4 5 This file is part of the GNU MP Library test suite. 6 7 The GNU MP Library test suite is free software; you can redistribute it 8 and/or modify it under the terms of the GNU General Public License as 9 published by the Free Software Foundation; either version 3 of the License, 10 or (at your option) any later version. 11 12 The GNU MP Library test suite is distributed in the hope that it will be 13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15 Public License for more details. 16 17 You should have received a copy of the GNU General Public License along with 18 the GNU MP Library test suite. If not, see https://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-impl.h" 26 #include "longlong.h" 27 #include "tests.h" 28 29 30 /* Change this to "#define TRACE(x) x" for some traces. */ 31 #define TRACE(x) 32 33 34 /* FIXME: Shouldn't use plain mpz functions in a reference routine. */ 35 void 36 refmpz_combit (mpz_ptr r, unsigned long bit) 37 { 38 if (mpz_tstbit (r, bit)) 39 mpz_clrbit (r, bit); 40 else 41 mpz_setbit (r, bit); 42 } 43 44 45 unsigned long 46 refmpz_hamdist (mpz_srcptr x, mpz_srcptr y) 47 { 48 mp_size_t xsize, ysize, tsize; 49 mp_ptr xp, yp; 50 unsigned long ret; 51 52 if ((SIZ(x) < 0 && SIZ(y) >= 0) 53 || (SIZ(y) < 0 && SIZ(x) >= 0)) 54 return ULONG_MAX; 55 56 xsize = ABSIZ(x); 57 ysize = ABSIZ(y); 58 tsize = MAX (xsize, ysize); 59 60 xp = refmpn_malloc_limbs (tsize); 61 refmpn_zero (xp, tsize); 62 refmpn_copy (xp, PTR(x), xsize); 63 64 yp = refmpn_malloc_limbs (tsize); 65 refmpn_zero (yp, tsize); 66 refmpn_copy (yp, PTR(y), ysize); 67 68 if (SIZ(x) < 0) 69 refmpn_neg (xp, xp, tsize); 70 71 if (SIZ(x) < 0) 72 refmpn_neg (yp, yp, tsize); 73 74 ret = refmpn_hamdist (xp, yp, tsize); 75 76 free (xp); 77 free (yp); 78 return ret; 79 } 80 81 void 82 refmpz_gcd (mpz_ptr g, mpz_srcptr a_orig, mpz_srcptr b_orig) 83 { 84 mp_bitcnt_t a_twos, b_twos, common_twos; 85 mpz_t a; 86 mpz_t b; 87 mpz_init (a); 88 mpz_init (b); 89 mpz_abs (a, a_orig); 90 mpz_abs (b, b_orig); 91 92 if (mpz_sgn (a) == 0) 93 { 94 mpz_set (g, b); 95 return; 96 } 97 if (mpz_sgn (b) == 0) 98 { 99 mpz_set (g, a); 100 return; 101 } 102 a_twos = mpz_scan1 (a, 0); 103 mpz_tdiv_q_2exp (a, a, a_twos); 104 105 b_twos = mpz_scan1 (b, 0); 106 mpz_tdiv_q_2exp (b, b, b_twos); 107 108 common_twos = MIN(a_twos, b_twos); 109 for (;;) 110 { 111 int c; 112 mp_bitcnt_t twos; 113 c = mpz_cmp (a, b); 114 if (c == 0) 115 break; 116 if (c < 0) 117 mpz_swap (a, b); 118 mpz_sub (a, a, b); 119 twos = mpz_scan1 (a, 0); 120 mpz_tdiv_q_2exp (a, a, twos); 121 } 122 mpz_mul_2exp (g, a, common_twos); 123 124 mpz_clear (a); 125 mpz_clear (b); 126 } 127 128 /* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */ 129 #define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b)) 130 131 /* (a/b) effect due to sign of b: mpz/mpz */ 132 #define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b)) 133 134 /* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd; 135 is (-1/b) if a<0, or +1 if a>=0 */ 136 #define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0]) 137 138 int 139 refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig) 140 { 141 unsigned long twos; 142 mpz_t a, b; 143 int result_bit1 = 0; 144 145 if (mpz_sgn (b_orig) == 0) 146 return JACOBI_Z0 (a_orig); /* (a/0) */ 147 148 if (mpz_sgn (a_orig) == 0) 149 return JACOBI_0Z (b_orig); /* (0/b) */ 150 151 if (mpz_even_p (a_orig) && mpz_even_p (b_orig)) 152 return 0; 153 154 if (mpz_cmp_ui (b_orig, 1) == 0) 155 return 1; 156 157 mpz_init_set (a, a_orig); 158 mpz_init_set (b, b_orig); 159 160 if (mpz_sgn (b) < 0) 161 { 162 result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b); 163 mpz_neg (b, b); 164 } 165 if (mpz_even_p (b)) 166 { 167 twos = mpz_scan1 (b, 0L); 168 mpz_tdiv_q_2exp (b, b, twos); 169 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); 170 } 171 172 if (mpz_sgn (a) < 0) 173 { 174 result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]); 175 mpz_neg (a, a); 176 } 177 if (mpz_even_p (a)) 178 { 179 twos = mpz_scan1 (a, 0L); 180 mpz_tdiv_q_2exp (a, a, twos); 181 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); 182 } 183 184 for (;;) 185 { 186 ASSERT (mpz_odd_p (a)); 187 ASSERT (mpz_odd_p (b)); 188 ASSERT (mpz_sgn (a) > 0); 189 ASSERT (mpz_sgn (b) > 0); 190 191 TRACE (printf ("top\n"); 192 mpz_trace (" a", a); 193 mpz_trace (" b", b)); 194 195 if (mpz_cmp (a, b) < 0) 196 { 197 TRACE (printf ("swap\n")); 198 mpz_swap (a, b); 199 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]); 200 } 201 202 if (mpz_cmp_ui (b, 1) == 0) 203 break; 204 205 mpz_sub (a, a, b); 206 TRACE (printf ("sub\n"); 207 mpz_trace (" a", a)); 208 if (mpz_sgn (a) == 0) 209 goto zero; 210 211 twos = mpz_scan1 (a, 0L); 212 mpz_fdiv_q_2exp (a, a, twos); 213 TRACE (printf ("twos %lu\n", twos); 214 mpz_trace (" a", a)); 215 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); 216 } 217 218 mpz_clear (a); 219 mpz_clear (b); 220 return JACOBI_BIT1_TO_PN (result_bit1); 221 222 zero: 223 mpz_clear (a); 224 mpz_clear (b); 225 return 0; 226 } 227 228 /* Same as mpz_kronecker, but ignoring factors of 2 on b */ 229 int 230 refmpz_jacobi (mpz_srcptr a, mpz_srcptr b) 231 { 232 ASSERT_ALWAYS (mpz_sgn (b) > 0); 233 ASSERT_ALWAYS (mpz_odd_p (b)); 234 235 return refmpz_kronecker (a, b); 236 } 237 238 /* Legendre symbol via powm. p must be an odd prime. */ 239 int 240 refmpz_legendre (mpz_srcptr a, mpz_srcptr p) 241 { 242 int res; 243 244 mpz_t r; 245 mpz_t e; 246 247 ASSERT_ALWAYS (mpz_sgn (p) > 0); 248 ASSERT_ALWAYS (mpz_odd_p (p)); 249 250 mpz_init (r); 251 mpz_init (e); 252 253 mpz_fdiv_r (r, a, p); 254 255 mpz_set (e, p); 256 mpz_sub_ui (e, e, 1); 257 mpz_fdiv_q_2exp (e, e, 1); 258 mpz_powm (r, r, e, p); 259 260 /* Normalize to a more or less symmetric range around zero */ 261 if (mpz_cmp (r, e) > 0) 262 mpz_sub (r, r, p); 263 264 ASSERT_ALWAYS (mpz_cmpabs_ui (r, 1) <= 0); 265 266 res = mpz_sgn (r); 267 268 mpz_clear (r); 269 mpz_clear (e); 270 271 return res; 272 } 273 274 275 int 276 refmpz_kronecker_ui (mpz_srcptr a, unsigned long b) 277 { 278 mpz_t bz; 279 int ret; 280 mpz_init_set_ui (bz, b); 281 ret = refmpz_kronecker (a, bz); 282 mpz_clear (bz); 283 return ret; 284 } 285 286 int 287 refmpz_kronecker_si (mpz_srcptr a, long b) 288 { 289 mpz_t bz; 290 int ret; 291 mpz_init_set_si (bz, b); 292 ret = refmpz_kronecker (a, bz); 293 mpz_clear (bz); 294 return ret; 295 } 296 297 int 298 refmpz_ui_kronecker (unsigned long a, mpz_srcptr b) 299 { 300 mpz_t az; 301 int ret; 302 mpz_init_set_ui (az, a); 303 ret = refmpz_kronecker (az, b); 304 mpz_clear (az); 305 return ret; 306 } 307 308 int 309 refmpz_si_kronecker (long a, mpz_srcptr b) 310 { 311 mpz_t az; 312 int ret; 313 mpz_init_set_si (az, a); 314 ret = refmpz_kronecker (az, b); 315 mpz_clear (az); 316 return ret; 317 } 318 319 320 void 321 refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e) 322 { 323 mpz_t s, t; 324 unsigned long i; 325 326 mpz_init_set_ui (t, 1L); 327 mpz_init_set (s, b); 328 329 if ((e & 1) != 0) 330 mpz_mul (t, t, s); 331 332 for (i = 2; i <= e; i <<= 1) 333 { 334 mpz_mul (s, s, s); 335 if ((i & e) != 0) 336 mpz_mul (t, t, s); 337 } 338 339 mpz_set (w, t); 340 341 mpz_clear (s); 342 mpz_clear (t); 343 } 344