diff options
Diffstat (limited to 'keygen/prime.c')
| -rw-r--r-- | keygen/prime.c | 762 |
1 files changed, 762 insertions, 0 deletions
diff --git a/keygen/prime.c b/keygen/prime.c new file mode 100644 index 00000000..d9bdebba --- /dev/null +++ b/keygen/prime.c @@ -0,0 +1,762 @@ +/* + * Prime generation. + */ + +#include <assert.h> +#include <math.h> + +#include "ssh.h" +#include "mpint.h" +#include "mpunsafe.h" +#include "sshkeygen.h" + +/* ---------------------------------------------------------------------- + * Standard probabilistic prime-generation algorithm: + * + * - get a number from our PrimeCandidateSource which will at least + * avoid being divisible by any prime under 2^16 + * + * - perform the Miller-Rabin primality test enough times to + * ensure the probability of it being composite is 2^-80 or + * less + * + * - go back to square one if any M-R test fails. + */ + +static PrimeGenerationContext *probprime_new_context( + const PrimeGenerationPolicy *policy) +{ + PrimeGenerationContext *ctx = snew(PrimeGenerationContext); + ctx->vt = policy; + return ctx; +} + +static void probprime_free_context(PrimeGenerationContext *ctx) +{ + sfree(ctx); +} + +static ProgressPhase probprime_add_progress_phase( + const PrimeGenerationPolicy *policy, + ProgressReceiver *prog, unsigned bits) +{ + /* + * The density of primes near x is 1/(log x). When x is about 2^b, + * that's 1/(b log 2). + * + * But we're only doing the expensive part of the process (the M-R + * checks) for a number that passes the initial winnowing test of + * having no factor less than 2^16 (at least, unless the prime is + * so small that PrimeCandidateSource gives up on that winnowing). + * The density of _those_ numbers is about 1/19.76. So the odds of + * hitting a prime per expensive attempt are boosted by a factor + * of 19.76. + */ + const double log_2 = 0.693147180559945309417232121458; + double winnow_factor = (bits < 32 ? 1.0 : 19.76); + double prob = winnow_factor / (bits * log_2); + + /* + * Estimate the cost of prime generation as the cost of the M-R + * modexps. + */ + double cost = (miller_rabin_checks_needed(bits) * + estimate_modexp_cost(bits)); + return progress_add_probabilistic(prog, cost, prob); +} + +static mp_int *probprime_generate( + PrimeGenerationContext *ctx, + PrimeCandidateSource *pcs, ProgressReceiver *prog) +{ + pcs_ready(pcs); + + while (true) { + progress_report_attempt(prog); + + mp_int *p = pcs_generate(pcs); + if (!p) { + pcs_free(pcs); + return NULL; + } + + MillerRabin *mr = miller_rabin_new(p); + bool known_bad = false; + unsigned nchecks = miller_rabin_checks_needed(mp_get_nbits(p)); + for (unsigned check = 0; check < nchecks; check++) { + if (!miller_rabin_test_random(mr)) { + known_bad = true; + break; + } + } + miller_rabin_free(mr); + + if (!known_bad) { + /* + * We have a prime! + */ + pcs_free(pcs); + return p; + } + + mp_free(p); + } +} + +static strbuf *null_mpu_certificate(PrimeGenerationContext *ctx, mp_int *p) +{ + return NULL; +} + +const PrimeGenerationPolicy primegen_probabilistic = { + probprime_add_progress_phase, + probprime_new_context, + probprime_free_context, + probprime_generate, + null_mpu_certificate, +}; + +/* ---------------------------------------------------------------------- + * Alternative provable-prime algorithm, based on the following paper: + * + * [MAURER] Maurer, U.M. Fast generation of prime numbers and secure + * public-key cryptographic parameters. J. Cryptology 8, 123–155 + * (1995). https://doi.org/10.1007/BF00202269 + */ + +typedef enum SubprimePolicy { + SPP_FAST, + SPP_MAURER_SIMPLE, + SPP_MAURER_COMPLEX, +} SubprimePolicy; + +typedef struct ProvablePrimePolicyExtra { + SubprimePolicy spp; +} ProvablePrimePolicyExtra; + +typedef struct ProvablePrimeContext ProvablePrimeContext; +struct ProvablePrimeContext { + Pockle *pockle; + PrimeGenerationContext pgc; + const ProvablePrimePolicyExtra *extra; +}; + +static PrimeGenerationContext *provableprime_new_context( + const PrimeGenerationPolicy *policy) +{ + ProvablePrimeContext *ppc = snew(ProvablePrimeContext); + ppc->pgc.vt = policy; + ppc->pockle = pockle_new(); + ppc->extra = policy->extra; + return &ppc->pgc; +} + +static void provableprime_free_context(PrimeGenerationContext *ctx) +{ + ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc); + pockle_free(ppc->pockle); + sfree(ppc); +} + +static ProgressPhase provableprime_add_progress_phase( + const PrimeGenerationPolicy *policy, + ProgressReceiver *prog, unsigned bits) +{ + /* + * Estimating the cost of making a _provable_ prime is difficult + * because of all the recursions to smaller sizes. + * + * Once you have enough factors of p-1 to certify primality of p, + * the remaining work in provable prime generation is not very + * different from probabilistic: you generate a random candidate, + * test its primality probabilistically, and use the witness value + * generated as a byproduct of that test for the full Pocklington + * verification. The expensive part, as usual, is made of modpows. + * + * The Pocklington test needs at least two modpows (one for the + * Fermat check, and one per known factor of p-1). + * + * The prior M-R step needs an unknown number, because we iterate + * until we find a value whose order is divisible by the largest + * power of 2 that divides p-1, say 2^j. That excludes half the + * possible witness values (specifically, the quadratic residues), + * so we expect to need on average two M-R operations to find one. + * But that's only if the number _is_ prime - as usual, it's also + * possible that we hit a non-prime and have to try again. + * + * So, if we were only estimating the cost of that final step, it + * would look a lot like the probabilistic version: we'd have to + * estimate the expected total number of modexps by knowing + * something about the density of primes among our candidate + * integers, and then multiply that by estimate_modexp_cost(bits). + * But the problem is that we also have to _find_ a smaller prime, + * so we have to recurse. + * + * In the MAURER_SIMPLE version of the algorithm, you recurse to + * any one of a range of possible smaller sizes i, each with + * probability proportional to 1/i. So your expected time to + * generate an n-bit prime is given by a horrible recurrence of + * the form E_n = S_n + (sum E_i/i) / (sum 1/i), in which S_n is + * the expected cost of the final step once you have your smaller + * primes, and both sums are over ceil(n/2) <= i <= n-20. + * + * At this point I ran out of effort to actually do the maths + * rigorously, so instead I did the empirical experiment of + * generating that sequence in Python and plotting it on a graph. + * My Python code is here, in case I need it again: + +from math import log + +alpha = log(3)/log(2) + 1 # exponent for modexp using Karatsuba mult + +E = [1] * 16 # assume generating tiny primes is trivial + +for n in range(len(E), 4096): + + # Expected time for sub-generations, as a weighted mean of prior + # values of the same sequence. + lo = (n+1)//2 + hi = n-20 + if lo <= hi: + subrange = range(lo, hi+1) + num = sum(E[i]/i for i in subrange) + den = sum(1/i for i in subrange) + else: + num, den = 0, 1 + + # Constant term (cost of final step). + # Similar to probprime_add_progress_phase. + winnow_factor = 1 if n < 32 else 19.76 + prob = winnow_factor / (n * log(2)) + cost = 4 * n**alpha / prob + + E.append(cost + num / den) + +for i, p in enumerate(E): + try: + print(log(i), log(p)) + except ValueError: + continue + + * The output loop prints the logs of both i and E_i, so that when + * I plot the resulting data file in gnuplot I get a log-log + * diagram. That showed me some early noise and then a very + * straight-looking line; feeding the straight part of the graph + * to linear-regression analysis reported that it fits the line + * + * log E_n = -1.7901825337965498 + 3.6199197179662517 * log(n) + * => E_n = 0.16692969657466802 * n^3.6199197179662517 + * + * So my somewhat empirical estimate is that Maurer prime + * generation costs about 0.167 * bits^3.62, in the same arbitrary + * time units used by estimate_modexp_cost. + */ + + return progress_add_linear(prog, 0.167 * pow(bits, 3.62)); +} + +static mp_int *primegen_small(Pockle *pockle, PrimeCandidateSource *pcs) +{ + assert(pcs_get_bits(pcs) <= 32); + + pcs_ready(pcs); + + while (true) { + mp_int *p = pcs_generate(pcs); + if (!p) { + pcs_free(pcs); + return NULL; + } + if (pockle_add_small_prime(pockle, p) == POCKLE_OK) { + pcs_free(pcs); + return p; + } + mp_free(p); + } +} + +#ifdef DEBUG_PRIMEGEN +static void timestamp(FILE *fp) +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + fprintf(fp, "%lu.%09lu: ", (unsigned long)ts.tv_sec, + (unsigned long)ts.tv_nsec); +} +static PRINTF_LIKE(1, 2) void debug_f(const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + timestamp(stderr); + vfprintf(stderr, fmt, ap); + fputc('\n', stderr); + va_end(ap); +} +static void debug_f_mp(const char *fmt, mp_int *x, ...) +{ + va_list ap; + va_start(ap, x); + timestamp(stderr); + vfprintf(stderr, fmt, ap); + mp_dump(stderr, "", x, "\n"); + va_end(ap); +} +#else +#define debug_f(...) ((void)0) +#define debug_f_mp(...) ((void)0) +#endif + +static double uniform_random_double(void) +{ + unsigned char randbuf[8]; + random_read(randbuf, 8); + return GET_64BIT_MSB_FIRST(randbuf) * 0x1.0p-64; +} + +static mp_int *mp_ceil_div(mp_int *n, mp_int *d) +{ + mp_int *nplus = mp_add(n, d); + mp_sub_integer_into(nplus, nplus, 1); + mp_int *toret = mp_div(nplus, d); + mp_free(nplus); + return toret; +} + +static mp_int *provableprime_generate_inner( + ProvablePrimeContext *ppc, PrimeCandidateSource *pcs, + ProgressReceiver *prog, double progress_origin, double progress_scale) +{ + unsigned bits = pcs_get_bits(pcs); + assert(bits > 1); + + if (bits <= 32) { + debug_f("ppgi(%u) -> small", bits); + return primegen_small(ppc->pockle, pcs); + } + + unsigned min_bits_needed, max_bits_needed; + { + /* + * Find the product of all the prime factors we already know + * about. + */ + mp_int *size_got = mp_from_integer(1); + size_t nfactors; + mp_int **factors = pcs_get_known_prime_factors(pcs, &nfactors); + for (size_t i = 0; i < nfactors; i++) { + mp_int *to_free = size_got; + size_got = mp_unsafe_shrink(mp_mul(size_got, factors[i])); + mp_free(to_free); + } + + /* + * Find the largest cofactor we might be able to use, and the + * smallest one we can get away with. + */ + mp_int *upperbound = pcs_get_upper_bound(pcs); + mp_int *size_needed = mp_nthroot(upperbound, 3, NULL); + debug_f_mp("upperbound = ", upperbound); + { + mp_int *to_free = upperbound; + upperbound = mp_unsafe_shrink(mp_div(upperbound, size_got)); + mp_free(to_free); + } + debug_f_mp("size_needed = ", size_needed); + { + mp_int *to_free = size_needed; + size_needed = mp_unsafe_shrink(mp_ceil_div(size_needed, size_got)); + mp_free(to_free); + } + + max_bits_needed = pcs_get_bits_remaining(pcs); + + /* + * We need a prime that is greater than or equal to + * 'size_needed' in order for the product of all our known + * factors of p-1 to exceed the cube root of the largest value + * p might take. + * + * Since pcs_new wants a size specified in bits, we must count + * the bits in size_needed and then add 1. Otherwise we might + * get a value with the same bit count as size_needed but + * slightly smaller than it. + * + * An exception is if size_needed = 1. In that case the + * product of existing known factors is _already_ enough, so + * we don't need to generate an extra factor at all. + */ + if (mp_hs_integer(size_needed, 2)) { + min_bits_needed = mp_get_nbits(size_needed) + 1; + } else { + min_bits_needed = 0; + } + + mp_free(upperbound); + mp_free(size_needed); + mp_free(size_got); + } + + double progress = 0.0; + + if (min_bits_needed) { + debug_f("ppgi(%u) recursing, need [%u,%u] more bits", + bits, min_bits_needed, max_bits_needed); + + unsigned *sizes = NULL; + size_t nsizes = 0, sizesize = 0; + + unsigned real_min = max_bits_needed / 2; + unsigned real_max = (max_bits_needed >= 20 ? + max_bits_needed - 20 : 0); + if (real_min < min_bits_needed) + real_min = min_bits_needed; + if (real_max < real_min) + real_max = real_min; + debug_f("ppgi(%u) revised bits interval = [%u,%u]", + bits, real_min, real_max); + + switch (ppc->extra->spp) { + case SPP_FAST: + /* + * Always pick the smallest subsidiary prime we can get + * away with: just over n/3 bits. + * + * This is not a good mode for cryptographic prime + * generation, because it skews the distribution of primes + * greatly, and worse, it skews them in a direction that + * heads away from the properties crypto algorithms tend + * to like. + * + * (For both discrete-log systems and RSA, people have + * tended to recommend in the past that p-1 should have a + * _large_ factor if possible. There's some disagreement + * on which algorithms this is really necessary for, but + * certainly I've never seen anyone recommend arranging a + * _small_ factor on purpose.) + * + * I originally implemented this mode because it was + * convenient for debugging - it wastes as little time as + * possible on finding a sub-prime and lets you get to the + * interesting part! And I leave it in the code because it + * might still be useful for _something_. Because it's + * cryptographically questionable, it's not selectable in + * the UI of either version of PuTTYgen proper; but it can + * be accessed through testcrypt, and if for some reason a + * definite prime is needed for non-crypto purposes, it + * may still be the fastest way to put your hands on one. + */ + debug_f("ppgi(%u) fast mode, just ask for %u bits", + bits, min_bits_needed); + sgrowarray(sizes, sizesize, nsizes); + sizes[nsizes++] = min_bits_needed; + break; + case SPP_MAURER_SIMPLE: { + /* + * Select the size of the subsidiary prime at random from + * sqrt(outputprime) up to outputprime/2^20, in such a way + * that the probability distribution matches that of the + * largest prime factor of a random n-bit number. + * + * Per [MAURER] section 3.4, the cumulative distribution + * function of this relative size is 1+log2(x), for x in + * [1/2,1]. You can generate a value from the distribution + * given by a cdf by applying the inverse cdf to a uniform + * value in [0,1]. Simplifying that in this case, what we + * have to do is raise 2 to the power of a random real + * number between -1 and 0. (And that gives you the number + * of _bits_ in the sub-prime, as a factor of the desired + * output number of bits.) + * + * We also require that the subsidiary prime q is at least + * 20 bits smaller than the output one, to give us a + * fighting chance of there being _any_ prime we can find + * such that q | p-1. + * + * (But these rules have to be applied in an order that + * still leaves us _some_ interval of possible sizes we + * can pick!) + */ + maurer_simple: + debug_f("ppgi(%u) Maurer simple mode", bits); + + unsigned sub_bits; + do { + double uniform = uniform_random_double(); + sub_bits = real_max * pow(2.0, uniform - 1) + 0.5; + debug_f(" ... %.6f -> %u?", uniform, sub_bits); + } while (!(real_min <= sub_bits && sub_bits <= real_max)); + + debug_f("ppgi(%u) asking for %u bits", bits, sub_bits); + sgrowarray(sizes, sizesize, nsizes); + sizes[nsizes++] = sub_bits; + + break; + } + case SPP_MAURER_COMPLEX: { + /* + * In this mode, we may generate multiple factors of p-1 + * which between them add up to at least n/2 bits, in such + * a way that those are guaranteed to be the largest + * factors of p-1 and that they have the same probability + * distribution as the largest k factors would have in a + * random integer. The idea is that this more elaborate + * procedure gets as close as possible to the same + * probability distribution you'd get by selecting a + * completely random prime (if you feasibly could). + * + * Algorithm from Appendix 1 of [MAURER]: we generate + * random real numbers that sum to at most 1, by choosing + * each one uniformly from the range [0, 1 - sum of all + * the previous ones]. We maintain them in a list in + * decreasing order, and we stop as soon as we find an + * initial subsequence of the list s_1,...,s_r such that + * s_1 + ... + s_{r-1} + 2 s_r > 1. In particular, this + * guarantees that the sum of that initial subsequence is + * at least 1/2, so we end up with enough factors to + * satisfy Pocklington. + */ + + if (max_bits_needed / 2 + 1 > real_max) { + /* Early exit path in the case where this algorithm + * can't possibly generate a value in the range we + * need. In that situation, fall back to Maurer + * simple. */ + debug_f("ppgi(%u) skipping GenerateSizeList, " + "real_max too small", bits); + goto maurer_simple; /* sorry! */ + } + + double *s = NULL; + size_t ns, ssize = 0; + + while (true) { + debug_f("ppgi(%u) starting GenerateSizeList", bits); + ns = 0; + double range = 1.0; + while (true) { + /* Generate the next number */ + double u = uniform_random_double() * range; + range -= u; + debug_f(" u_%"SIZEu" = %g", ns, u); + + /* Insert it in the list */ + sgrowarray(s, ssize, ns); + size_t i; + for (i = ns; i > 0 && s[i-1] < u; i--) + s[i] = s[i-1]; + s[i] = u; + ns++; + debug_f(" inserting as s[%"SIZEu"]", i); + + /* Look for a suitable initial subsequence */ + double sum = 0; + for (i = 0; i < ns; i++) { + sum += s[i]; + if (sum + s[i] > 1.0) { + debug_f(" s[0..%"SIZEu"] works!", i); + + /* Truncate the sequence here, and stop + * generating random real numbers. */ + ns = i+1; + goto got_list; + } + } + } + + got_list:; + /* + * Now translate those real numbers into actual bit + * counts, and do a last-minute check to make sure + * their product is going to be in range. + * + * We have to check both the min and max sizes of the + * total. A b-bit number is in [2^{b-1},2^b). So the + * product of numbers of sizes b_1,...,b_k is at least + * 2^{\sum (b_i-1)}, and less than 2^{\sum b_i}. + */ + nsizes = 0; + + unsigned min_total = 0, max_total = 0; + + for (size_t i = 0; i < ns; i++) { + /* These sizes are measured in actual entropy, so + * add 1 bit each time to account for the + * zero-information leading 1 */ + unsigned this_size = max_bits_needed * s[i] + 1; + debug_f(" bits[%"SIZEu"] = %u", i, this_size); + sgrowarray(sizes, sizesize, nsizes); + sizes[nsizes++] = this_size; + + min_total += this_size - 1; + max_total += this_size; + } + + debug_f(" total bits = [%u,%u)", min_total, max_total); + if (min_total < real_min || max_total > real_max+1) { + debug_f(" total out of range, try again"); + } else { + debug_f(" success! %"SIZEu" sub-primes totalling [%u,%u) " + "bits", nsizes, min_total, max_total); + break; + } + } + + smemclr(s, ssize * sizeof(*s)); + sfree(s); + break; + } + default: + unreachable("bad subprime policy"); + } + + for (size_t i = 0; i < nsizes; i++) { + unsigned sub_bits = sizes[i]; + double progress_in_this_prime = (double)sub_bits / bits; + mp_int *q = provableprime_generate_inner( + ppc, pcs_new(sub_bits), + prog, progress_origin + progress_scale * progress, + progress_scale * progress_in_this_prime); + progress += progress_in_this_prime; + assert(q); + debug_f_mp("ppgi(%u) got factor ", q, bits); + pcs_require_residue_1_mod_prime(pcs, q); + mp_free(q); + } + + smemclr(sizes, sizesize * sizeof(*sizes)); + sfree(sizes); + } else { + debug_f("ppgi(%u) no need to recurse", bits); + } + + debug_f("ppgi(%u) ready, %u bits remaining", + bits, pcs_get_bits_remaining(pcs)); + pcs_ready(pcs); + + while (true) { + mp_int *p = pcs_generate(pcs); + if (!p) { + pcs_free(pcs); + return NULL; + } + + debug_f_mp("provable_step p=", p); + + MillerRabin *mr = miller_rabin_new(p); + debug_f("provable_step mr setup done"); + mp_int *witness = miller_rabin_find_potential_primitive_root(mr); + miller_rabin_free(mr); + + if (!witness) { + debug_f("provable_step mr failed"); + mp_free(p); + continue; + } + + size_t nfactors; + mp_int **factors = pcs_get_known_prime_factors(pcs, &nfactors); + PockleStatus st = pockle_add_prime( + ppc->pockle, p, factors, nfactors, witness); + + if (st != POCKLE_OK) { + debug_f("provable_step proof failed %d", (int)st); + + /* + * Check by assertion that the error status is not one of + * the ones we ought to have ruled out already by + * construction. If there's a bug in this code that means + * we can _never_ pass this test (e.g. picking products of + * factors that never quite reach cbrt(n)), we'd rather + * fail an assertion than loop forever. + */ + assert(st == POCKLE_DISCRIMINANT_IS_SQUARE || + st == POCKLE_WITNESS_POWER_IS_1 || + st == POCKLE_WITNESS_POWER_NOT_COPRIME); + + mp_free(p); + if (witness) + mp_free(witness); + continue; + } + + mp_free(witness); + pcs_free(pcs); + debug_f_mp("ppgi(%u) done, got ", p, bits); + progress_report(prog, progress_origin + progress_scale); + return p; + } +} + +static mp_int *provableprime_generate( + PrimeGenerationContext *ctx, + PrimeCandidateSource *pcs, ProgressReceiver *prog) +{ + ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc); + mp_int *p = provableprime_generate_inner(ppc, pcs, prog, 0.0, 1.0); + + return p; +} + +static inline strbuf *provableprime_mpu_certificate( + PrimeGenerationContext *ctx, mp_int *p) +{ + ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc); + return pockle_mpu(ppc->pockle, p); +} + +#define DECLARE_POLICY(name, policy) \ + static const struct ProvablePrimePolicyExtra \ + pppextra_##name = {policy}; \ + const PrimeGenerationPolicy name = { \ + provableprime_add_progress_phase, \ + provableprime_new_context, \ + provableprime_free_context, \ + provableprime_generate, \ + provableprime_mpu_certificate, \ + &pppextra_##name, \ + } + +DECLARE_POLICY(primegen_provable_fast, SPP_FAST); +DECLARE_POLICY(primegen_provable_maurer_simple, SPP_MAURER_SIMPLE); +DECLARE_POLICY(primegen_provable_maurer_complex, SPP_MAURER_COMPLEX); + +/* ---------------------------------------------------------------------- + * Reusable null implementation of the progress-reporting API. + */ + +static inline ProgressPhase null_progress_add(void) { + ProgressPhase ph = { .n = 0 }; + return ph; +} +ProgressPhase null_progress_add_linear( + ProgressReceiver *prog, double c) { return null_progress_add(); } +ProgressPhase null_progress_add_probabilistic( + ProgressReceiver *prog, double c, double p) { return null_progress_add(); } +void null_progress_ready(ProgressReceiver *prog) {} +void null_progress_start_phase(ProgressReceiver *prog, ProgressPhase phase) {} +void null_progress_report(ProgressReceiver *prog, double progress) {} +void null_progress_report_attempt(ProgressReceiver *prog) {} +void null_progress_report_phase_complete(ProgressReceiver *prog) {} +const ProgressReceiverVtable null_progress_vt = { + .add_linear = null_progress_add_linear, + .add_probabilistic = null_progress_add_probabilistic, + .ready = null_progress_ready, + .start_phase = null_progress_start_phase, + .report = null_progress_report, + .report_attempt = null_progress_report_attempt, + .report_phase_complete = null_progress_report_phase_complete, +}; + +/* ---------------------------------------------------------------------- + * Helper function for progress estimation. + */ + +double estimate_modexp_cost(unsigned bits) +{ + /* + * A modexp of n bits goes roughly like O(n^2.58), on the grounds + * that our modmul is O(n^1.58) (Karatsuba) and you need O(n) of + * them in a modexp. + */ + return pow(bits, 2.58); +} |