#include #include #include #include #include #include #include "ptypes.h" #include "random_prime.h" #include "utility.h" #include "primality.h" #include "gmp_main.h" #include "isaac.h" #include "prime_iterator.h" static char pr[31] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127}; void mpz_random_nbit_prime(mpz_t p, UV n) { switch (n) { case 0: case 1: mpz_set_ui(p, 0); return; case 2: mpz_set_ui(p, pr[ 0+isaac_rand(2)]); return; case 3: mpz_set_ui(p, pr[ 2+isaac_rand(2)]); return; case 4: mpz_set_ui(p, pr[ 4+isaac_rand(2)]); return; case 5: mpz_set_ui(p, pr[ 6+isaac_rand(5)]); return; case 6: mpz_set_ui(p, pr[11+isaac_rand(7)]); return; case 7: mpz_set_ui(p, pr[18+isaac_rand(13)]); return; default: break; } /* For 32-bit inputs, use fast trivial method */ if (n <= 32) { uint32_t mask = (0xFFFFFFFFU >> (34-n)) << 1, base = mask+3; do { mpz_set_ui(p, base | (isaac_rand32() & mask)); } while (!_GMP_is_prob_prime(p)); } else { #if 0 do { /* Trivial method. */ mpz_isaac_urandomb(p, n); mpz_setbit(p, n-1); mpz_setbit(p, 0); } while (!_GMP_is_prob_prime(p)); #else mpz_t base; /* Fouque+Tibouchi Alg 1, without modulo checks */ mpz_init(base); if (n > 33) { mpz_isaac_urandomb(base, n-33); mpz_mul_2exp(base,base,1); } mpz_setbit(base, n-1); mpz_setbit(base, 0); do { mpz_set_ui(p, isaac_rand32()); mpz_mul_2exp(p, p, n-32); mpz_ior(p, p, base); } while (!_GMP_is_prob_prime(p)); mpz_clear(base); #endif } } /* PRIMEINC: pick random value, select next prime. */ /* Fast but bad distribution. */ static int _random_prime_primeinc(mpz_t p, mpz_t lo, mpz_t hi) { mpz_t r, t; mpz_init(t); mpz_init(r); mpz_sub(r, hi, lo); mpz_isaac_urandomm(t, r); mpz_clear(r); mpz_add(t, t, lo); mpz_sub_ui(t, t, 1); _GMP_next_prime(t); if (mpz_cmp(t,hi) > 0) { mpz_sub_ui(t, lo, 1); _GMP_next_prime(t); if (mpz_cmp(t,hi) > 0) { mpz_clear(t); return 0; } } mpz_set(p, t); mpz_clear(t); return 1; } /* TRIVIAL: pick random values until one is prime */ /* Perfect distribution. */ static int _random_prime_trivial(mpz_t p, mpz_t lo_in, mpz_t hi_in) { mpz_t r, t, lo, hi; int res = 0, tries = 10000; if (mpz_cmp_ui(hi_in,2) < 0 || mpz_cmp(lo_in,hi_in) > 0) return 0; mpz_init_set(lo, lo_in); mpz_init_set(hi, hi_in); if (mpz_cmp_ui(lo,2) <= 0) { mpz_set_ui(lo,1); } else if (mpz_even_p(lo)) { mpz_add_ui(lo,lo,1); } if (mpz_cmp_ui(hi,2) <= 0) { mpz_set_ui(hi,1); } else if (mpz_even_p(hi)) { mpz_sub_ui(hi,hi,1); } /* lo and hi are now odd */ if (mpz_cmp(lo,hi) >= 0) { if (mpz_cmp(lo,hi) > 0) { /* null range */ } else if (mpz_cmp_ui(lo,1) == 0) { mpz_set_ui(p,2); res = 1; } else if (_GMP_is_prob_prime(lo)) { mpz_set(p,lo); res = 1; } mpz_clear(hi); mpz_clear(lo); return res; } /* lo and hi are now odd and at least one odd between them */ mpz_init(t); mpz_init(r); mpz_sub(r, hi, lo); mpz_tdiv_q_2exp(r, r, 1); mpz_add_ui(r,r,1); do { mpz_isaac_urandomm(t, r); mpz_mul_2exp(t, t, 1); mpz_add(t, t, lo); if (mpz_cmp_ui(t,1) == 0) mpz_set_ui(t,2); /* map 1 back to 2 */ } while (!_GMP_is_prob_prime(t) && --tries > 0); if (tries > 0) { mpz_set(p, t); res = 1; } else { /* We couldn't find anything. Perhaps no primes in range. */ res = _random_prime_primeinc(p, lo, hi); } mpz_clear(r); mpz_clear(t); mpz_clear(hi); mpz_clear(lo); return res; } /* Set p to a random prime between lo and hi inclusive */ int mpz_random_prime(mpz_t p, mpz_t lo, mpz_t hi) { return _random_prime_trivial(p,lo,hi); } void mpz_random_ndigit_prime(mpz_t p, UV n) { mpz_t lo, hi; switch (n) { case 0: mpz_set_ui(p,0); return; case 1: mpz_set_ui(p, pr[isaac_rand(4)]); return; case 2: mpz_set_ui(p, pr[4+isaac_rand(21)]); return; default: break; } mpz_init_set_ui(lo,10); mpz_pow_ui(lo, lo, n-1); mpz_init(hi); mpz_mul_ui(hi, lo, 10); if (!mpz_random_prime(p, lo, hi)) croak("Failed to find %"UVuf" digit prime\n", n); mpz_clear(lo); mpz_clear(hi); } /* Random number rop such that 2*mult*rop+1 has nbits bits. */ static void _rand_in_bit_interval(mpz_t rop, UV nbits, mpz_t mult) { mpz_t t, lo, hi; mpz_init(t); mpz_init(lo); mpz_init(hi); mpz_mul_ui(t, mult, 2); mpz_setbit(lo, nbits-1); mpz_sub_ui(lo, lo, 1); mpz_cdiv_q(lo, lo, t); /* lo = ceil(2^(nbits-1)-1 / (2*mult)) */ mpz_setbit(hi, nbits); mpz_sub_ui(hi, hi, 2); mpz_fdiv_q(hi, hi, t); /* hi = floor(2^nbits-2 / (2*mult)) */ mpz_sub(t, hi, lo); mpz_isaac_urandomm(rop, t); mpz_add(rop, rop, lo); mpz_clear(t); mpz_clear(lo); mpz_clear(hi); } #define _SAFE_REJECT(q, p1, p2, p3, p4, p5) \ { uint32_t qm = mpz_fdiv_ui(q, p1*p2*p3*p4*p5); \ if ((qm % p1) == 0 || (qm % p1) == (p1>>1)) continue; \ if ((qm % p2) == 0 || (qm % p2) == (p2>>1)) continue; \ if ((qm % p3) == 0 || (qm % p3) == (p3>>1)) continue; \ if ((qm % p4) == 0 || (qm % p4) == (p4>>1)) continue; \ if ((qm % p5) == 0 || (qm % p5) == (p5>>1)) continue; \ } void mpz_random_safe_prime(mpz_t p, UV nbits) { static const unsigned char small_safe[] = {5,7,11,23,47,59,83,107}; mpz_t q, base; uint32_t qmod, PR, tlimit; int verbose; PRIME_ITERATOR(iter); switch (nbits) { case 0: case 1: case 2: mpz_set_ui(p, 0); return; case 3: mpz_set_ui(p, small_safe[ 0+isaac_rand(2)]); return; case 4: mpz_set_ui(p, 11); return; case 5: mpz_set_ui(p, 23); return; case 6: mpz_set_ui(p, small_safe[ 4+isaac_rand(2)]); return; case 7: mpz_set_ui(p, small_safe[ 6+isaac_rand(2)]); return; default: break; } mpz_init(q); mpz_init(base); if (nbits > 35) { mpz_isaac_urandomb(base, nbits-35); mpz_mul_2exp(base,base,2); } mpz_setbit(base, nbits-1); mpz_setbit(base, 1); mpz_setbit(base, 0); verbose = get_verbose_level(); tlimit = (nbits <= 512000) ? (nbits*(nbits/64.0) + 0.5) : 4000000000U; while (1) { /* 1. Generate random nbit p */ if (nbits > 35) { mpz_set_ui(p, isaac_rand32()); mpz_mul_2exp(p, p, nbits-33); } else { mpz_isaac_urandomb(p, nbits-3); mpz_mul_2exp(p, p, 2); } mpz_ior(p, p, base); /* 2. p = 2q+1 => q = p >> 1 */ mpz_div_2exp(q, p, 1); /* 3. Fast compositeness pretests for both q and p at the same time. */ qmod = mpz_fdiv_ui(q, 1155UL); if ( (qmod % 3) != 2 || (qmod % 5) == 0 || (qmod % 7) == 0 || (qmod % 11) == 0 || (qmod % 5) == 2 || (qmod % 7) == 3 || (qmod % 11) == 5) continue; if (nbits < 16) { /* 4. Pretest that p isn't easily composite */ if (!primality_pretest(p)) continue; /* 5. Pretest that q isn't easily composite */ if (!primality_pretest(q)) continue; } else { _SAFE_REJECT(q, 13U, 17U, 19U, 23U, 29U); _SAFE_REJECT(q, 31U, 37U, 41U, 43U, 47U); _SAFE_REJECT(q, 53U, 59U, 61U, 67U, 71U); _SAFE_REJECT(q, 73U, 79U, 83U, 89U, 97U); if (tlimit >= 101) { prime_iterator_setprime(&iter, PR = 101); while (PR <= tlimit) { uint32_t qm = mpz_fdiv_ui(q, PR); if (qm == 0 || qm == (PR>>1)) break; PR = prime_iterator_next(&iter); } if (PR <= tlimit) continue; } } if (verbose > 2) { printf("."); fflush(stdout); } /* 6. BPSW on q and M-R base 2 on p. */ if (!is_euler_plumb_pseudoprime(q)) continue; /* Start faster */ if (verbose > 2) { printf("+"); fflush(stdout); } if (!miller_rabin_ui(p, 2)) continue; if (verbose > 2) { printf("*"); fflush(stdout); } if (!_GMP_is_lucas_pseudoprime(q, 2)) continue; if (nbits > 64 && !miller_rabin_ui(q, 2)) continue; /* Verify fast test */ break; } mpz_clear(base); mpz_clear(q); prime_iterator_destroy(&iter); } /* Gordon's algorithm */ void mpz_random_strong_prime(mpz_t p, UV nbits) { mpz_t S, T, R, P0, t, i, j; UV rbits, sbits, tbits; if (nbits < 128) croak("random_strong_prime, bits must be >= 128"); if (nbits < 256) { rbits = ((nbits+1) >> 1) - 2; sbits = (nbits >> 1) - 20; tbits = rbits - 20; } else { UV N1, N2; { /* Calculate FIPS 186-4 C.10 recommended parameter */ UV t_, l2_; for (l2_ = 1, t_ = nbits; t_ >>= 1; ) l2_++; N1 = (nbits/2)-l2_-7; N2 = N1/2; } if (N1 > 200) N1 = 201; if (N2 > 100) N2 = 101; if (N2 < 100) N2 += N1/4; rbits = sbits = N1; tbits = N2; } mpz_init(S); mpz_init(T); mpz_init(R); mpz_init(P0); mpz_init(t); mpz_init(i); mpz_init(j); while (1) { mpz_random_nbit_prime(S, sbits); mpz_random_nbit_prime(T, tbits); _rand_in_bit_interval(i, rbits, T); while (1) { mpz_mul(t, i, T); mpz_mul_ui(t, t, 2); mpz_add_ui(R, t, 1); /* R = 2*i*T+1 */ if (_GMP_is_prob_prime(R)) break; mpz_add_ui(i,i,1); } mpz_sub_ui(t, R, 2); mpz_powm(P0, S, t, R); mpz_mul_ui(P0, P0, 2); mpz_mul(P0, P0, S); mpz_sub_ui(P0, P0, 1); mpz_mul(i, R, S); mpz_mul_ui(t, i, 2); _rand_in_bit_interval(j, nbits, i); while (1) { mpz_mul(p, j, t); mpz_add(p, p, P0); /* p = 2*j*R*S+p0 */ if (mpz_sizeinbase(p,2) > nbits) break; if (_GMP_is_prob_prime(p)) { mpz_clear(t); mpz_clear(i); mpz_clear(j); mpz_clear(S); mpz_clear(T); mpz_clear(R); mpz_clear(P0); /* p-1 has factor R. p+1 has factor S. r-1 has factor T. */ return; } mpz_add_ui(j,j,1); } } } /*===========================================================================*/ /* Proven primes (Maurer and Shawe-Taylor */ /*===========================================================================*/ #define MAKE_PROOF_START(proofptr, n, nums) \ if (proofptr) { \ char* thisproof, *thisptr; \ int prevlen = (*proofptr == 0) ? 0 : strlen(*proofptr); \ int thislen = (5 + mpz_sizeinbase(n,10)) * nums + 200; \ New(0, thisproof, thislen + prevlen + 1, char); \ thisptr = thisproof; \ thisptr += gmp_sprintf(thisptr, #define MAKE_PROOF_END(proofptr) \ ); \ if (*proofptr) { \ thisptr += gmp_sprintf(thisptr,"\n"); \ strcat(thisptr, *proofptr); \ Safefree(*proofptr); \ } \ *proofptr = thisproof; \ } #define USE_THEOREM5 0 void mpz_random_maurer_prime(mpz_t n, UV k, char** proofptr) { mpz_t t, a, q, I, R; double m, r, minr = USE_THEOREM5 ? 0.334 : 0.5; int i, verbose = get_verbose_level(); /* We could use safely use k <= 64. */ if (k <= 32) return mpz_random_nbit_prime(n, k); r = minr; /* size of q relative to size of n */ m = 20; /* always use at least this many bits of randomness */ if (k > 2*m) { do { double s = ((double)isaac_rand32()) / ((double)4294967295.0); /* [0,1] */ r = pow(2,s-1); /* exp2 is C99 */ #if USE_THEOREM5 r = 0.334 + 1.332 * (r-0.5); /* Stretch r to cover 0.334 - 1 */ #endif } while ((k-r*k) <= m); } #if 0 /* Improve efficiency for less than ideal distribution */ r -= 0.25; if (r < minr) r = minr; #endif mpz_init(t); mpz_init(a); mpz_init(q); mpz_init(I); mpz_init(R); mpz_random_maurer_prime(q, (UV)(r*k)+1, proofptr); mpz_setbit(I, k-1); mpz_mul_ui(t, q, 2); mpz_fdiv_q(I, I, t); /* I = floor(2^(k-1) / 2q) */ if (verbose && verbose != 3) { gmp_printf("r = %lf k = %lu q = %Zd I = %Zd\n",r,k,q,I); fflush(stdout); } while (1) { if (verbose > 2) { printf("."); fflush(stdout); } mpz_isaac_urandomm(R, I); /* [0, I-1] */ mpz_add(R, R, I); /* [I, 2I-I] */ mpz_add_ui(R, R, 1); /* [I+1,2I] */ #if USE_THEOREM5 mpz_setbit(R, 0); /* We need R to be odd */ #endif mpz_mul(n, R, q); mpz_mul_ui(n,n,2); mpz_add_ui(n,n,1); /* n = 2Rq+1 */ if (!primality_pretest(n)) continue; if (verbose > 2) { printf("+"); fflush(stdout); } /* if (!is_euler_plumb_pseudoprime(n)) continue; */ if (!miller_rabin_ui(n,2)) continue; /* n is a base-2 psp and probably prime */ if (verbose > 2) { printf("*"); fflush(stdout); } /* See if we can use BLS75 theorem 3 */ mpz_mul_ui(t, q, 2); mpz_add_ui(t, t, 1); mpz_mul(t, t, t); if (mpz_cmp(t, n) > 0) { for (i = 0; i < 20; i++) { mpz_set_ui(a, pr[i]); /* Check A^R mod N != N-1 */ mpz_powm(a, a, R, n); mpz_add_ui(t,a,1); if (mpz_cmp(t, n) == 0) continue; /* Check A^{Rq} mod N == N-1 */ mpz_powm(a, a, q, n); mpz_add_ui(t,a,1); if (mpz_cmp(t, n) != 0) continue; if (verbose > 2) { printf("(%"UVuf")",k); fflush(stdout); } /* Ensure all results passed BPSW. ~20% speed penalty. */ if (!_GMP_is_lucas_pseudoprime(n,2)) croak("Maurer internal failure"); MAKE_PROOF_START(proofptr, n, 3) "Type BLS3\nN %Zd\nQ %Zd\nA %u\n", n, q, pr[i] MAKE_PROOF_END(proofptr) mpz_clear(t); mpz_clear(a); mpz_clear(q); mpz_clear(I); mpz_clear(R); return; } /* Blast, we couldn't find the right 'a' value fast enough. Try a new n. */ continue; } /* Our q is smaller than sqrt(n)/2-1, so use BLS75 theorem 5. */ #if !USE_THEOREM5 croak("random_maurer_prime: internal bit size error"); #else /* Check for obvious generation problems. */ if (mpz_even_p(R)) continue; if (mpz_cmp_ui(R, 1) <= 0) continue; mpz_gcd(t, q, R); if (mpz_cmp_ui(t, 1) != 0) continue; /* Theorem 5 with m = 2, assuming (I) which we'll check after this. */ { mpz_t ts, tr, F; mpz_init(ts); mpz_init(tr); mpz_init(F); mpz_mul_ui(F, q, 2); /* Calculate r,s from page 624 of BLS75 */ mpz_mul_ui(t, F, 2); mpz_tdiv_qr(ts, tr, R, t); /* Verify the r,s condition */ mpz_mul(t, tr, tr); mpz_submul_ui(t, ts, 8); /* t = r^2-8s */ if (mpz_sgn(ts) != 0 && mpz_perfect_square_p(t)) { /* printf("fail r/s check\n"); */ mpz_clear(ts); mpz_clear(tr); mpz_clear(F); continue; } /* Verify size of N with m=2. a,t are temps. Should not fail. */ mpz_mul(t, F, tr); mpz_add_ui(a, t, 1); /* a = rF + 1 */ mpz_sub_ui(tr, tr, 1); mpz_mul(t, F, F); mpz_mul_ui(t, t, 2); mpz_mul(t, t, tr); mpz_add(a, a, t); /* a = (r-1)2FF + rF + 1 */ mpz_mul(t, F, F); mpz_mul(t, t, F); mpz_mul_ui(t, t, 4); mpz_add(a, a, t); /* a = 4FFF + (r-1)2FF + rF + 1 */ mpz_add_ui(t, F, 1); mpz_clear(tr); mpz_clear(ts); mpz_clear(F); if (mpz_cmp(n,a) >= 0) { /* printf("fail N size check\n"); */ continue; } /* Check divisibility required to use m=2 */ if (mpz_divisible_p(n,t)) { /* printf("fail N divisiblity check\n"); */ continue; } } #define SET_A_CHECK_PSP(i) \ mpz_set_ui(a, pr[i]); \ if (apsp[i] == -1) \ { mpz_sub_ui(t,n,1); mpz_powm(t,a,t,n); apsp[i] = (mpz_cmp_ui(t,1) == 0); } \ if (apsp[i] == 0) continue; #define CHECK_GCD(t) \ mpz_powm(t, a, t, n); mpz_sub_ui(t, t, 1); mpz_gcd(t, t, n); \ if (mpz_cmp_ui(t,1) != 0) continue; { int j, apsp[20]; /* apsp caches psp check. Init all to -1. */ for (i = 0; i < 20; i++) apsp[i] = -1; apsp[0] = 1; /* We passed a base 2 psp test to get here */ /* Find an a that works for p=2 */ for (i = 0; i < 20; i++) { SET_A_CHECK_PSP(i); mpz_mul(t, q, R); CHECK_GCD(t); /* We are good for p=2. Find an a for p=q. */ for (j = 0; j < 20; j++) { SET_A_CHECK_PSP(j); mpz_mul_ui(t, R, 2); CHECK_GCD(t); /* Success */ if (verbose > 2) { printf("(%lu)",k); fflush(stdout); } if (i == 0 && j == 0) { MAKE_PROOF_START(proofptr, n, 2) "Type BLS5\nN %Zd\nQ[1] %Zd\n----\n", n, q MAKE_PROOF_END(proofptr) } else { MAKE_PROOF_START(proofptr, n, 2) "Type BLS5\nN %Zd\nQ[1] %Zd\nA[0] %lu\nA[1] %lu\n----\n", n, q, pr[i], pr[j] MAKE_PROOF_END(proofptr) } /* Ensure all results passed BPSW. ~20% speed penalty. */ if (!_GMP_is_lucas_pseudoprime(n,2)) croak("Maurer internal failure"); mpz_clear(t); mpz_clear(a); mpz_clear(q); mpz_clear(I); mpz_clear(R); return; } break; /* Failed for p=q */ } } /* Blast, we couldn't find the right 'a' value fast enough. Try a new n. */ #endif } } /* FIPS 186-4 algorithm but using our CSPRNG (ISAAC) instead of SHA-256 */ void mpz_random_shawe_taylor_prime(mpz_t c, UV k, char** proofptr) { mpz_t c0, t, u, a, z; if (k <= 32) return mpz_random_nbit_prime(c, k); mpz_init(c0); mpz_init(t); mpz_init(u); mpz_init(a); mpz_init(z); mpz_random_shawe_taylor_prime(c0, 1 + (k+1)/2, proofptr); mpz_isaac_urandomb(t, k-1); mpz_setbit(t,k-1); /* Steps 18-21: t a random k-bit integer */ mpz_mul_ui(u, c0, 2); mpz_cdiv_q(t, t, u); /* Step 22: set t based on random integer */ while (1) { /* Steps 23-24 */ mpz_mul_ui(u, c0, 2); /* u = 2c0 */ mpz_mul(c, u, t); mpz_add_ui(c, c, 1); /* c = 2tc0+1 */ if (mpz_sizeinbase(c,2) > k) { mpz_set_ui(t,0); mpz_setbit(t,k-1); mpz_cdiv_q(t, t, u); mpz_mul(c, u, t); mpz_add_ui(c, c, 1); } /* Don't bother with Steps 26-31 for obvious composites */ if (primality_pretest(c) && miller_rabin_ui(c,2)) { /* Steps 26-29 */ mpz_sub_ui(u, c, 3); mpz_isaac_urandomm(a, u); mpz_add_ui(a, a, 2); /* Step 30 */ mpz_mul_ui(u, t, 2); mpz_powm(z, a, u, c); /* Step 31 */ mpz_sub_ui(u, z, 1); mpz_gcd(u, u, c); if (mpz_cmp_ui(u, 1) == 0) { mpz_powm(u, z, c0, c); if (mpz_cmp_ui(u, 1) == 0) { /* Ensure all results passed BPSW. ~20% speed penalty. */ if (!_GMP_is_lucas_pseudoprime(c,2)) croak("ST internal failure"); MAKE_PROOF_START(proofptr, c, 3) "Type Pocklington\nN %Zd\nQ %Zd\nA %Zd\n", c, c0, a MAKE_PROOF_END(proofptr) mpz_clear(c0); mpz_clear(t); mpz_clear(u); mpz_clear(a); mpz_clear(z); return; } } } mpz_add_ui(t,t,1); } } /*===========================================================================*/