Ciao, Il Dom, 19 Febbraio 2017 9:21 am, Niels Möller ha scritto: > But shifting, as you suggest, might be simpler, using > > 2 B^623 = 20023 (mod p)
A possible generalization follows: #include <assert.h> #include <stdio.h> #include <stdlib.h> #include <gmp.h> #define GSIZE (19937 / GMP_NUMB_BITS) #define GSHIFT (19937 % GMP_NUMB_BITS) #define SIZE (GSIZE + 1) #define K 20023 /* B^312 mod p */ #define B312MODP ((mp_limb_t) K << (GMP_NUMB_BITS-GSHIFT)) #define USE_NIELS_SHORTCUT (((mp_limb_t) K >> GSHIFT) == 0) /* Reduces r modulo p in place, result is SIZE limbs (non-canonical representation). */ static void reduce (mp_ptr rp, mp_size_t n) { mp_limb_t cy, hi; assert (n > GSIZE); assert (GSHIFT != 0); /* 19937 is prime */ if (USE_NIELS_SHORTCUT) { if (n == SIZE) return; while (n >= 2*SIZE) { rp[n - SIZE] = mpn_addmul_1 (rp + n - 2*SIZE, rp + n - SIZE, SIZE, B312MODP); n -= GSIZE; } cy = mpn_addmul_1 (rp, rp + SIZE, n - SIZE, B312MODP); if (n < 2*SIZE - 1) cy = mpn_add_1 (rp + n - SIZE, rp + n - SIZE, 2*SIZE - n - 1, cy); hi = rp[GSIZE]; rp[GSIZE] = cy + (hi & (((mp_limb_t)1 << GSHIFT) - 1)) + mpn_add_1 (rp, rp, SIZE - 1, (hi >> GSHIFT) * K); } else { while (n >= 2 * GSIZE) { hi = mpn_rshift (rp + n - GSIZE, rp + n - GSIZE, GSIZE, GSHIFT); cy = mpn_addmul_1 (rp + n - 2*GSIZE, rp + n - GSIZE, GSIZE, K); rp[n - GSIZE] = cy + (hi >> (GMP_NUMB_BITS - GSHIFT)); n -= GSIZE - 1; } hi = mpn_rshift (rp + GSIZE, rp + GSIZE, n - GSIZE, GSHIFT); cy = mpn_addmul_1 (rp, rp + GSIZE, n - GSIZE, K); cy = mpn_add_1 (rp + n - GSIZE, rp + n - GSIZE, 2*GSIZE - n, cy); rp[GSIZE] = cy + (hi >> (GMP_NUMB_BITS - GSHIFT)); } } int main (int argc, char **argv) { mpz_t p, x, r, ref; int i; gmp_randstate_t rands; mpz_inits (p, x, r, ref, NULL); gmp_randinit_default (rands); mpz_setbit (p, 19937); mpz_sub_ui (p, p, K); for (i = 0; i < 64000; i++) { mpz_urandomb (x, rands, 16); mpz_rrandomb (x, rands, mpz_get_ui (x) + 20000); mpz_tdiv_r (ref, x, p); mpz_set (r, x); reduce (mpz_limbs_modify (r, mpz_size (r)), mpz_size(r)); mpz_limbs_finish (r, SIZE); if (!mpz_congruent_p (r, ref, p)) { gmp_printf ("Failed (%i) for size %d: %Zx\n" " got: %Zx\n" " exp: %Zx\n", i, mpz_size (x), x, r, ref); exit (1); } } mpz_clears (p, x, r, ref, NULL); gmp_randclear(rands); } Best regards, m -- http://bodrato.it/ _______________________________________________ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs