1 /*
2 * Prime generation.
3 */
4
5 #include <assert.h>
6 #include <math.h>
7
8 #include "ssh.h"
9 #include "mpint.h"
10 #include "mpunsafe.h"
11 #include "sshkeygen.h"
12
13 /* ----------------------------------------------------------------------
14 * Standard probabilistic prime-generation algorithm:
15 *
16 * - get a number from our PrimeCandidateSource which will at least
17 * avoid being divisible by any prime under 2^16
18 *
19 * - perform the Miller-Rabin primality test enough times to
20 * ensure the probability of it being composite is 2^-80 or
21 * less
22 *
23 * - go back to square one if any M-R test fails.
24 */
25
probprime_new_context(const PrimeGenerationPolicy * policy)26 static PrimeGenerationContext *probprime_new_context(
27 const PrimeGenerationPolicy *policy)
28 {
29 PrimeGenerationContext *ctx = snew(PrimeGenerationContext);
30 ctx->vt = policy;
31 return ctx;
32 }
33
probprime_free_context(PrimeGenerationContext * ctx)34 static void probprime_free_context(PrimeGenerationContext *ctx)
35 {
36 sfree(ctx);
37 }
38
probprime_add_progress_phase(const PrimeGenerationPolicy * policy,ProgressReceiver * prog,unsigned bits)39 static ProgressPhase probprime_add_progress_phase(
40 const PrimeGenerationPolicy *policy,
41 ProgressReceiver *prog, unsigned bits)
42 {
43 /*
44 * The density of primes near x is 1/(log x). When x is about 2^b,
45 * that's 1/(b log 2).
46 *
47 * But we're only doing the expensive part of the process (the M-R
48 * checks) for a number that passes the initial winnowing test of
49 * having no factor less than 2^16 (at least, unless the prime is
50 * so small that PrimeCandidateSource gives up on that winnowing).
51 * The density of _those_ numbers is about 1/19.76. So the odds of
52 * hitting a prime per expensive attempt are boosted by a factor
53 * of 19.76.
54 */
55 const double log_2 = 0.693147180559945309417232121458;
56 double winnow_factor = (bits < 32 ? 1.0 : 19.76);
57 double prob = winnow_factor / (bits * log_2);
58
59 /*
60 * Estimate the cost of prime generation as the cost of the M-R
61 * modexps.
62 */
63 double cost = (miller_rabin_checks_needed(bits) *
64 estimate_modexp_cost(bits));
65 return progress_add_probabilistic(prog, cost, prob);
66 }
67
probprime_generate(PrimeGenerationContext * ctx,PrimeCandidateSource * pcs,ProgressReceiver * prog)68 static mp_int *probprime_generate(
69 PrimeGenerationContext *ctx,
70 PrimeCandidateSource *pcs, ProgressReceiver *prog)
71 {
72 pcs_ready(pcs);
73
74 while (true) {
75 progress_report_attempt(prog);
76
77 mp_int *p = pcs_generate(pcs);
78 if (!p) {
79 pcs_free(pcs);
80 return NULL;
81 }
82
83 MillerRabin *mr = miller_rabin_new(p);
84 bool known_bad = false;
85 unsigned nchecks = miller_rabin_checks_needed(mp_get_nbits(p));
86 for (unsigned check = 0; check < nchecks; check++) {
87 if (!miller_rabin_test_random(mr)) {
88 known_bad = true;
89 break;
90 }
91 }
92 miller_rabin_free(mr);
93
94 if (!known_bad) {
95 /*
96 * We have a prime!
97 */
98 pcs_free(pcs);
99 return p;
100 }
101
102 mp_free(p);
103 }
104 }
105
null_mpu_certificate(PrimeGenerationContext * ctx,mp_int * p)106 static strbuf *null_mpu_certificate(PrimeGenerationContext *ctx, mp_int *p)
107 {
108 return NULL;
109 }
110
111 const PrimeGenerationPolicy primegen_probabilistic = {
112 probprime_add_progress_phase,
113 probprime_new_context,
114 probprime_free_context,
115 probprime_generate,
116 null_mpu_certificate,
117 };
118
119 /* ----------------------------------------------------------------------
120 * Alternative provable-prime algorithm, based on the following paper:
121 *
122 * [MAURER] Maurer, U.M. Fast generation of prime numbers and secure
123 * public-key cryptographic parameters. J. Cryptology 8, 123–155
124 * (1995). https://doi.org/10.1007/BF00202269
125 */
126
127 typedef enum SubprimePolicy {
128 SPP_FAST,
129 SPP_MAURER_SIMPLE,
130 SPP_MAURER_COMPLEX,
131 } SubprimePolicy;
132
133 typedef struct ProvablePrimePolicyExtra {
134 SubprimePolicy spp;
135 } ProvablePrimePolicyExtra;
136
137 typedef struct ProvablePrimeContext ProvablePrimeContext;
138 struct ProvablePrimeContext {
139 Pockle *pockle;
140 PrimeGenerationContext pgc;
141 const ProvablePrimePolicyExtra *extra;
142 };
143
provableprime_new_context(const PrimeGenerationPolicy * policy)144 static PrimeGenerationContext *provableprime_new_context(
145 const PrimeGenerationPolicy *policy)
146 {
147 ProvablePrimeContext *ppc = snew(ProvablePrimeContext);
148 ppc->pgc.vt = policy;
149 ppc->pockle = pockle_new();
150 ppc->extra = policy->extra;
151 return &ppc->pgc;
152 }
153
provableprime_free_context(PrimeGenerationContext * ctx)154 static void provableprime_free_context(PrimeGenerationContext *ctx)
155 {
156 ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc);
157 pockle_free(ppc->pockle);
158 sfree(ppc);
159 }
160
provableprime_add_progress_phase(const PrimeGenerationPolicy * policy,ProgressReceiver * prog,unsigned bits)161 static ProgressPhase provableprime_add_progress_phase(
162 const PrimeGenerationPolicy *policy,
163 ProgressReceiver *prog, unsigned bits)
164 {
165 /*
166 * Estimating the cost of making a _provable_ prime is difficult
167 * because of all the recursions to smaller sizes.
168 *
169 * Once you have enough factors of p-1 to certify primality of p,
170 * the remaining work in provable prime generation is not very
171 * different from probabilistic: you generate a random candidate,
172 * test its primality probabilistically, and use the witness value
173 * generated as a byproduct of that test for the full Pocklington
174 * verification. The expensive part, as usual, is made of modpows.
175 *
176 * The Pocklington test needs at least two modpows (one for the
177 * Fermat check, and one per known factor of p-1).
178 *
179 * The prior M-R step needs an unknown number, because we iterate
180 * until we find a value whose order is divisible by the largest
181 * power of 2 that divides p-1, say 2^j. That excludes half the
182 * possible witness values (specifically, the quadratic residues),
183 * so we expect to need on average two M-R operations to find one.
184 * But that's only if the number _is_ prime - as usual, it's also
185 * possible that we hit a non-prime and have to try again.
186 *
187 * So, if we were only estimating the cost of that final step, it
188 * would look a lot like the probabilistic version: we'd have to
189 * estimate the expected total number of modexps by knowing
190 * something about the density of primes among our candidate
191 * integers, and then multiply that by estimate_modexp_cost(bits).
192 * But the problem is that we also have to _find_ a smaller prime,
193 * so we have to recurse.
194 *
195 * In the MAURER_SIMPLE version of the algorithm, you recurse to
196 * any one of a range of possible smaller sizes i, each with
197 * probability proportional to 1/i. So your expected time to
198 * generate an n-bit prime is given by a horrible recurrence of
199 * the form E_n = S_n + (sum E_i/i) / (sum 1/i), in which S_n is
200 * the expected cost of the final step once you have your smaller
201 * primes, and both sums are over ceil(n/2) <= i <= n-20.
202 *
203 * At this point I ran out of effort to actually do the maths
204 * rigorously, so instead I did the empirical experiment of
205 * generating that sequence in Python and plotting it on a graph.
206 * My Python code is here, in case I need it again:
207
208 from math import log
209
210 alpha = log(3)/log(2) + 1 # exponent for modexp using Karatsuba mult
211
212 E = [1] * 16 # assume generating tiny primes is trivial
213
214 for n in range(len(E), 4096):
215
216 # Expected time for sub-generations, as a weighted mean of prior
217 # values of the same sequence.
218 lo = (n+1)//2
219 hi = n-20
220 if lo <= hi:
221 subrange = range(lo, hi+1)
222 num = sum(E[i]/i for i in subrange)
223 den = sum(1/i for i in subrange)
224 else:
225 num, den = 0, 1
226
227 # Constant term (cost of final step).
228 # Similar to probprime_add_progress_phase.
229 winnow_factor = 1 if n < 32 else 19.76
230 prob = winnow_factor / (n * log(2))
231 cost = 4 * n**alpha / prob
232
233 E.append(cost + num / den)
234
235 for i, p in enumerate(E):
236 try:
237 print(log(i), log(p))
238 except ValueError:
239 continue
240
241 * The output loop prints the logs of both i and E_i, so that when
242 * I plot the resulting data file in gnuplot I get a log-log
243 * diagram. That showed me some early noise and then a very
244 * straight-looking line; feeding the straight part of the graph
245 * to linear-regression analysis reported that it fits the line
246 *
247 * log E_n = -1.7901825337965498 + 3.6199197179662517 * log(n)
248 * => E_n = 0.16692969657466802 * n^3.6199197179662517
249 *
250 * So my somewhat empirical estimate is that Maurer prime
251 * generation costs about 0.167 * bits^3.62, in the same arbitrary
252 * time units used by estimate_modexp_cost.
253 */
254
255 return progress_add_linear(prog, 0.167 * pow(bits, 3.62));
256 }
257
primegen_small(Pockle * pockle,PrimeCandidateSource * pcs)258 static mp_int *primegen_small(Pockle *pockle, PrimeCandidateSource *pcs)
259 {
260 assert(pcs_get_bits(pcs) <= 32);
261
262 pcs_ready(pcs);
263
264 while (true) {
265 mp_int *p = pcs_generate(pcs);
266 if (!p) {
267 pcs_free(pcs);
268 return NULL;
269 }
270 if (pockle_add_small_prime(pockle, p) == POCKLE_OK) {
271 pcs_free(pcs);
272 return p;
273 }
274 mp_free(p);
275 }
276 }
277
278 #ifdef DEBUG_PRIMEGEN
timestamp(FILE * fp)279 static void timestamp(FILE *fp)
280 {
281 struct timespec ts;
282 clock_gettime(CLOCK_MONOTONIC, &ts);
283 fprintf(fp, "%lu.%09lu: ", (unsigned long)ts.tv_sec,
284 (unsigned long)ts.tv_nsec);
285 }
debug_f(const char * fmt,...)286 static PRINTF_LIKE(1, 2) void debug_f(const char *fmt, ...)
287 {
288 va_list ap;
289 va_start(ap, fmt);
290 timestamp(stderr);
291 vfprintf(stderr, fmt, ap);
292 fputc('\n', stderr);
293 va_end(ap);
294 }
debug_f_mp(const char * fmt,mp_int * x,...)295 static void debug_f_mp(const char *fmt, mp_int *x, ...)
296 {
297 va_list ap;
298 va_start(ap, x);
299 timestamp(stderr);
300 vfprintf(stderr, fmt, ap);
301 mp_dump(stderr, "", x, "\n");
302 va_end(ap);
303 }
304 #else
305 #define debug_f(...) ((void)0)
306 #define debug_f_mp(...) ((void)0)
307 #endif
308
uniform_random_double(void)309 static double uniform_random_double(void)
310 {
311 unsigned char randbuf[8];
312 random_read(randbuf, 8);
313 return GET_64BIT_MSB_FIRST(randbuf) * 0x1.0p-64;
314 }
315
mp_ceil_div(mp_int * n,mp_int * d)316 static mp_int *mp_ceil_div(mp_int *n, mp_int *d)
317 {
318 mp_int *nplus = mp_add(n, d);
319 mp_sub_integer_into(nplus, nplus, 1);
320 mp_int *toret = mp_div(nplus, d);
321 mp_free(nplus);
322 return toret;
323 }
324
provableprime_generate_inner(ProvablePrimeContext * ppc,PrimeCandidateSource * pcs,ProgressReceiver * prog,double progress_origin,double progress_scale)325 static mp_int *provableprime_generate_inner(
326 ProvablePrimeContext *ppc, PrimeCandidateSource *pcs,
327 ProgressReceiver *prog, double progress_origin, double progress_scale)
328 {
329 unsigned bits = pcs_get_bits(pcs);
330 assert(bits > 1);
331
332 if (bits <= 32) {
333 debug_f("ppgi(%u) -> small", bits);
334 return primegen_small(ppc->pockle, pcs);
335 }
336
337 unsigned min_bits_needed, max_bits_needed;
338 {
339 /*
340 * Find the product of all the prime factors we already know
341 * about.
342 */
343 mp_int *size_got = mp_from_integer(1);
344 size_t nfactors;
345 mp_int **factors = pcs_get_known_prime_factors(pcs, &nfactors);
346 for (size_t i = 0; i < nfactors; i++) {
347 mp_int *to_free = size_got;
348 size_got = mp_unsafe_shrink(mp_mul(size_got, factors[i]));
349 mp_free(to_free);
350 }
351
352 /*
353 * Find the largest cofactor we might be able to use, and the
354 * smallest one we can get away with.
355 */
356 mp_int *upperbound = pcs_get_upper_bound(pcs);
357 mp_int *size_needed = mp_nthroot(upperbound, 3, NULL);
358 debug_f_mp("upperbound = ", upperbound);
359 {
360 mp_int *to_free = upperbound;
361 upperbound = mp_unsafe_shrink(mp_div(upperbound, size_got));
362 mp_free(to_free);
363 }
364 debug_f_mp("size_needed = ", size_needed);
365 {
366 mp_int *to_free = size_needed;
367 size_needed = mp_unsafe_shrink(mp_ceil_div(size_needed, size_got));
368 mp_free(to_free);
369 }
370
371 max_bits_needed = pcs_get_bits_remaining(pcs);
372
373 /*
374 * We need a prime that is greater than or equal to
375 * 'size_needed' in order for the product of all our known
376 * factors of p-1 to exceed the cube root of the largest value
377 * p might take.
378 *
379 * Since pcs_new wants a size specified in bits, we must count
380 * the bits in size_needed and then add 1. Otherwise we might
381 * get a value with the same bit count as size_needed but
382 * slightly smaller than it.
383 *
384 * An exception is if size_needed = 1. In that case the
385 * product of existing known factors is _already_ enough, so
386 * we don't need to generate an extra factor at all.
387 */
388 if (mp_hs_integer(size_needed, 2)) {
389 min_bits_needed = mp_get_nbits(size_needed) + 1;
390 } else {
391 min_bits_needed = 0;
392 }
393
394 mp_free(upperbound);
395 mp_free(size_needed);
396 mp_free(size_got);
397 }
398
399 double progress = 0.0;
400
401 if (min_bits_needed) {
402 debug_f("ppgi(%u) recursing, need [%u,%u] more bits",
403 bits, min_bits_needed, max_bits_needed);
404
405 unsigned *sizes = NULL;
406 size_t nsizes = 0, sizesize = 0;
407
408 unsigned real_min = max_bits_needed / 2;
409 unsigned real_max = (max_bits_needed >= 20 ?
410 max_bits_needed - 20 : 0);
411 if (real_min < min_bits_needed)
412 real_min = min_bits_needed;
413 if (real_max < real_min)
414 real_max = real_min;
415 debug_f("ppgi(%u) revised bits interval = [%u,%u]",
416 bits, real_min, real_max);
417
418 switch (ppc->extra->spp) {
419 case SPP_FAST:
420 /*
421 * Always pick the smallest subsidiary prime we can get
422 * away with: just over n/3 bits.
423 *
424 * This is not a good mode for cryptographic prime
425 * generation, because it skews the distribution of primes
426 * greatly, and worse, it skews them in a direction that
427 * heads away from the properties crypto algorithms tend
428 * to like.
429 *
430 * (For both discrete-log systems and RSA, people have
431 * tended to recommend in the past that p-1 should have a
432 * _large_ factor if possible. There's some disagreement
433 * on which algorithms this is really necessary for, but
434 * certainly I've never seen anyone recommend arranging a
435 * _small_ factor on purpose.)
436 *
437 * I originally implemented this mode because it was
438 * convenient for debugging - it wastes as little time as
439 * possible on finding a sub-prime and lets you get to the
440 * interesting part! And I leave it in the code because it
441 * might still be useful for _something_. Because it's
442 * cryptographically questionable, it's not selectable in
443 * the UI of either version of PuTTYgen proper; but it can
444 * be accessed through testcrypt, and if for some reason a
445 * definite prime is needed for non-crypto purposes, it
446 * may still be the fastest way to put your hands on one.
447 */
448 debug_f("ppgi(%u) fast mode, just ask for %u bits",
449 bits, min_bits_needed);
450 sgrowarray(sizes, sizesize, nsizes);
451 sizes[nsizes++] = min_bits_needed;
452 break;
453 case SPP_MAURER_SIMPLE: {
454 /*
455 * Select the size of the subsidiary prime at random from
456 * sqrt(outputprime) up to outputprime/2^20, in such a way
457 * that the probability distribution matches that of the
458 * largest prime factor of a random n-bit number.
459 *
460 * Per [MAURER] section 3.4, the cumulative distribution
461 * function of this relative size is 1+log2(x), for x in
462 * [1/2,1]. You can generate a value from the distribution
463 * given by a cdf by applying the inverse cdf to a uniform
464 * value in [0,1]. Simplifying that in this case, what we
465 * have to do is raise 2 to the power of a random real
466 * number between -1 and 0. (And that gives you the number
467 * of _bits_ in the sub-prime, as a factor of the desired
468 * output number of bits.)
469 *
470 * We also require that the subsidiary prime q is at least
471 * 20 bits smaller than the output one, to give us a
472 * fighting chance of there being _any_ prime we can find
473 * such that q | p-1.
474 *
475 * (But these rules have to be applied in an order that
476 * still leaves us _some_ interval of possible sizes we
477 * can pick!)
478 */
479 maurer_simple:
480 debug_f("ppgi(%u) Maurer simple mode", bits);
481
482 unsigned sub_bits;
483 do {
484 double uniform = uniform_random_double();
485 sub_bits = real_max * pow(2.0, uniform - 1) + 0.5;
486 debug_f(" ... %.6f -> %u?", uniform, sub_bits);
487 } while (!(real_min <= sub_bits && sub_bits <= real_max));
488
489 debug_f("ppgi(%u) asking for %u bits", bits, sub_bits);
490 sgrowarray(sizes, sizesize, nsizes);
491 sizes[nsizes++] = sub_bits;
492
493 break;
494 }
495 case SPP_MAURER_COMPLEX: {
496 /*
497 * In this mode, we may generate multiple factors of p-1
498 * which between them add up to at least n/2 bits, in such
499 * a way that those are guaranteed to be the largest
500 * factors of p-1 and that they have the same probability
501 * distribution as the largest k factors would have in a
502 * random integer. The idea is that this more elaborate
503 * procedure gets as close as possible to the same
504 * probability distribution you'd get by selecting a
505 * completely random prime (if you feasibly could).
506 *
507 * Algorithm from Appendix 1 of [MAURER]: we generate
508 * random real numbers that sum to at most 1, by choosing
509 * each one uniformly from the range [0, 1 - sum of all
510 * the previous ones]. We maintain them in a list in
511 * decreasing order, and we stop as soon as we find an
512 * initial subsequence of the list s_1,...,s_r such that
513 * s_1 + ... + s_{r-1} + 2 s_r > 1. In particular, this
514 * guarantees that the sum of that initial subsequence is
515 * at least 1/2, so we end up with enough factors to
516 * satisfy Pocklington.
517 */
518
519 if (max_bits_needed / 2 + 1 > real_max) {
520 /* Early exit path in the case where this algorithm
521 * can't possibly generate a value in the range we
522 * need. In that situation, fall back to Maurer
523 * simple. */
524 debug_f("ppgi(%u) skipping GenerateSizeList, "
525 "real_max too small", bits);
526 goto maurer_simple; /* sorry! */
527 }
528
529 double *s = NULL;
530 size_t ns, ssize = 0;
531
532 while (true) {
533 debug_f("ppgi(%u) starting GenerateSizeList", bits);
534 ns = 0;
535 double range = 1.0;
536 while (true) {
537 /* Generate the next number */
538 double u = uniform_random_double() * range;
539 range -= u;
540 debug_f(" u_%"SIZEu" = %g", ns, u);
541
542 /* Insert it in the list */
543 sgrowarray(s, ssize, ns);
544 size_t i;
545 for (i = ns; i > 0 && s[i-1] < u; i--)
546 s[i] = s[i-1];
547 s[i] = u;
548 ns++;
549 debug_f(" inserting as s[%"SIZEu"]", i);
550
551 /* Look for a suitable initial subsequence */
552 double sum = 0;
553 for (i = 0; i < ns; i++) {
554 sum += s[i];
555 if (sum + s[i] > 1.0) {
556 debug_f(" s[0..%"SIZEu"] works!", i);
557
558 /* Truncate the sequence here, and stop
559 * generating random real numbers. */
560 ns = i+1;
561 goto got_list;
562 }
563 }
564 }
565
566 got_list:;
567 /*
568 * Now translate those real numbers into actual bit
569 * counts, and do a last-minute check to make sure
570 * their product is going to be in range.
571 *
572 * We have to check both the min and max sizes of the
573 * total. A b-bit number is in [2^{b-1},2^b). So the
574 * product of numbers of sizes b_1,...,b_k is at least
575 * 2^{\sum (b_i-1)}, and less than 2^{\sum b_i}.
576 */
577 nsizes = 0;
578
579 unsigned min_total = 0, max_total = 0;
580
581 for (size_t i = 0; i < ns; i++) {
582 /* These sizes are measured in actual entropy, so
583 * add 1 bit each time to account for the
584 * zero-information leading 1 */
585 unsigned this_size = max_bits_needed * s[i] + 1;
586 debug_f(" bits[%"SIZEu"] = %u", i, this_size);
587 sgrowarray(sizes, sizesize, nsizes);
588 sizes[nsizes++] = this_size;
589
590 min_total += this_size - 1;
591 max_total += this_size;
592 }
593
594 debug_f(" total bits = [%u,%u)", min_total, max_total);
595 if (min_total < real_min || max_total > real_max+1) {
596 debug_f(" total out of range, try again");
597 } else {
598 debug_f(" success! %"SIZEu" sub-primes totalling [%u,%u) "
599 "bits", nsizes, min_total, max_total);
600 break;
601 }
602 }
603
604 smemclr(s, ssize * sizeof(*s));
605 sfree(s);
606 break;
607 }
608 default:
609 unreachable("bad subprime policy");
610 }
611
612 for (size_t i = 0; i < nsizes; i++) {
613 unsigned sub_bits = sizes[i];
614 double progress_in_this_prime = (double)sub_bits / bits;
615 mp_int *q = provableprime_generate_inner(
616 ppc, pcs_new(sub_bits),
617 prog, progress_origin + progress_scale * progress,
618 progress_scale * progress_in_this_prime);
619 progress += progress_in_this_prime;
620 assert(q);
621 debug_f_mp("ppgi(%u) got factor ", q, bits);
622 pcs_require_residue_1_mod_prime(pcs, q);
623 mp_free(q);
624 }
625
626 smemclr(sizes, sizesize * sizeof(*sizes));
627 sfree(sizes);
628 } else {
629 debug_f("ppgi(%u) no need to recurse", bits);
630 }
631
632 debug_f("ppgi(%u) ready, %u bits remaining",
633 bits, pcs_get_bits_remaining(pcs));
634 pcs_ready(pcs);
635
636 while (true) {
637 mp_int *p = pcs_generate(pcs);
638 if (!p) {
639 pcs_free(pcs);
640 return NULL;
641 }
642
643 debug_f_mp("provable_step p=", p);
644
645 MillerRabin *mr = miller_rabin_new(p);
646 debug_f("provable_step mr setup done");
647 mp_int *witness = miller_rabin_find_potential_primitive_root(mr);
648 miller_rabin_free(mr);
649
650 if (!witness) {
651 debug_f("provable_step mr failed");
652 mp_free(p);
653 continue;
654 }
655
656 size_t nfactors;
657 mp_int **factors = pcs_get_known_prime_factors(pcs, &nfactors);
658 PockleStatus st = pockle_add_prime(
659 ppc->pockle, p, factors, nfactors, witness);
660
661 if (st != POCKLE_OK) {
662 debug_f("provable_step proof failed %d", (int)st);
663
664 /*
665 * Check by assertion that the error status is not one of
666 * the ones we ought to have ruled out already by
667 * construction. If there's a bug in this code that means
668 * we can _never_ pass this test (e.g. picking products of
669 * factors that never quite reach cbrt(n)), we'd rather
670 * fail an assertion than loop forever.
671 */
672 assert(st == POCKLE_DISCRIMINANT_IS_SQUARE ||
673 st == POCKLE_WITNESS_POWER_IS_1 ||
674 st == POCKLE_WITNESS_POWER_NOT_COPRIME);
675
676 mp_free(p);
677 if (witness)
678 mp_free(witness);
679 continue;
680 }
681
682 mp_free(witness);
683 pcs_free(pcs);
684 debug_f_mp("ppgi(%u) done, got ", p, bits);
685 progress_report(prog, progress_origin + progress_scale);
686 return p;
687 }
688 }
689
provableprime_generate(PrimeGenerationContext * ctx,PrimeCandidateSource * pcs,ProgressReceiver * prog)690 static mp_int *provableprime_generate(
691 PrimeGenerationContext *ctx,
692 PrimeCandidateSource *pcs, ProgressReceiver *prog)
693 {
694 ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc);
695 mp_int *p = provableprime_generate_inner(ppc, pcs, prog, 0.0, 1.0);
696
697 return p;
698 }
699
provableprime_mpu_certificate(PrimeGenerationContext * ctx,mp_int * p)700 static inline strbuf *provableprime_mpu_certificate(
701 PrimeGenerationContext *ctx, mp_int *p)
702 {
703 ProvablePrimeContext *ppc = container_of(ctx, ProvablePrimeContext, pgc);
704 return pockle_mpu(ppc->pockle, p);
705 }
706
707 #define DECLARE_POLICY(name, policy) \
708 static const struct ProvablePrimePolicyExtra \
709 pppextra_##name = {policy}; \
710 const PrimeGenerationPolicy name = { \
711 provableprime_add_progress_phase, \
712 provableprime_new_context, \
713 provableprime_free_context, \
714 provableprime_generate, \
715 provableprime_mpu_certificate, \
716 &pppextra_##name, \
717 }
718
719 DECLARE_POLICY(primegen_provable_fast, SPP_FAST);
720 DECLARE_POLICY(primegen_provable_maurer_simple, SPP_MAURER_SIMPLE);
721 DECLARE_POLICY(primegen_provable_maurer_complex, SPP_MAURER_COMPLEX);
722
723 /* ----------------------------------------------------------------------
724 * Reusable null implementation of the progress-reporting API.
725 */
726
null_progress_add(void)727 static inline ProgressPhase null_progress_add(void) {
728 ProgressPhase ph = { .n = 0 };
729 return ph;
730 }
null_progress_add_linear(ProgressReceiver * prog,double c)731 ProgressPhase null_progress_add_linear(
732 ProgressReceiver *prog, double c) { return null_progress_add(); }
null_progress_add_probabilistic(ProgressReceiver * prog,double c,double p)733 ProgressPhase null_progress_add_probabilistic(
734 ProgressReceiver *prog, double c, double p) { return null_progress_add(); }
null_progress_ready(ProgressReceiver * prog)735 void null_progress_ready(ProgressReceiver *prog) {}
null_progress_start_phase(ProgressReceiver * prog,ProgressPhase phase)736 void null_progress_start_phase(ProgressReceiver *prog, ProgressPhase phase) {}
null_progress_report(ProgressReceiver * prog,double progress)737 void null_progress_report(ProgressReceiver *prog, double progress) {}
null_progress_report_attempt(ProgressReceiver * prog)738 void null_progress_report_attempt(ProgressReceiver *prog) {}
null_progress_report_phase_complete(ProgressReceiver * prog)739 void null_progress_report_phase_complete(ProgressReceiver *prog) {}
740 const ProgressReceiverVtable null_progress_vt = {
741 .add_linear = null_progress_add_linear,
742 .add_probabilistic = null_progress_add_probabilistic,
743 .ready = null_progress_ready,
744 .start_phase = null_progress_start_phase,
745 .report = null_progress_report,
746 .report_attempt = null_progress_report_attempt,
747 .report_phase_complete = null_progress_report_phase_complete,
748 };
749
750 /* ----------------------------------------------------------------------
751 * Helper function for progress estimation.
752 */
753
estimate_modexp_cost(unsigned bits)754 double estimate_modexp_cost(unsigned bits)
755 {
756 /*
757 * A modexp of n bits goes roughly like O(n^2.58), on the grounds
758 * that our modmul is O(n^1.58) (Karatsuba) and you need O(n) of
759 * them in a modexp.
760 */
761 return pow(bits, 2.58);
762 }
763