ni...@lysator.liu.se (Niels Möller) writes: > ni...@lysator.liu.se (Niels Möller) writes: > > Or if we want to take advantage of the structure, we need an mpn > function to reduce numbers modulo 2^19937 - 20023.
Below is a sketch for the 64-bit case, not yet working. These things are a bit tricky to get right, but it's not very complex code either. Regards, /Niels #include <assert.h> #include <stdio.h> #include <stdlib.h> #include <gmp.h> #if GMP_NUMB_BITS != 64 #error unsupported limb size #endif #define SIZE 312 #define K 20023 /* B^312 mod p */ #define B312MODP ((mp_limb_t) K << 31) /* Reduces r modulo p in place, result is 312 limbs (non-canonical representation). */ static void reduce (mp_ptr rp, mp_size_t n) { mp_limb_t cy, hi; assert (n >= SIZE); if (n == SIZE) return; while (n >= 2*SIZE) { rp[n - SIZE] = mpn_addmul_1 (rp + n - 2*SIZE, rp + n - SIZE, SIZE, B312MODP); n -= (SIZE-1); } 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[SIZE - 1]; rp[SIZE - 1] = cy + (hi & (((mp_limb_t)1<<31) - 1)) + mpn_add_1 (rp, rp, SIZE - 1, (hi >> 31) * K); } 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_set_ui (p, 1); mpz_mul_2exp (p, p, 19937); mpz_sub_ui (p, p, K); for (i = 0; i < 1000; 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 for size %d: %Zx\n" " got: %Zx\n" " exp: %Zx\n", mpz_size (x), x, r, ref); exit (1); } } mpz_clears (p, x, r, ref, NULL); gmp_randclear(rands); } -- Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677. Internet email is subject to wholesale government surveillance. _______________________________________________ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs