1*b1066cf3Schristos /* $NetBSD: moduli.c,v 1.17 2023/07/26 17:58:15 christos Exp $ */
2*b1066cf3Schristos /* $OpenBSD: moduli.c,v 1.39 2023/03/02 06:41:56 dtucker Exp $ */
3ca32bd8dSchristos /*
4ca32bd8dSchristos * Copyright 1994 Phil Karn <karn@qualcomm.com>
5ca32bd8dSchristos * Copyright 1996-1998, 2003 William Allen Simpson <wsimpson@greendragon.com>
6ca32bd8dSchristos * Copyright 2000 Niels Provos <provos@citi.umich.edu>
7ca32bd8dSchristos * All rights reserved.
8ca32bd8dSchristos *
9ca32bd8dSchristos * Redistribution and use in source and binary forms, with or without
10ca32bd8dSchristos * modification, are permitted provided that the following conditions
11ca32bd8dSchristos * are met:
12ca32bd8dSchristos * 1. Redistributions of source code must retain the above copyright
13ca32bd8dSchristos * notice, this list of conditions and the following disclaimer.
14ca32bd8dSchristos * 2. Redistributions in binary form must reproduce the above copyright
15ca32bd8dSchristos * notice, this list of conditions and the following disclaimer in the
16ca32bd8dSchristos * documentation and/or other materials provided with the distribution.
17ca32bd8dSchristos *
18ca32bd8dSchristos * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
19ca32bd8dSchristos * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20ca32bd8dSchristos * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21ca32bd8dSchristos * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
22ca32bd8dSchristos * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23ca32bd8dSchristos * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24ca32bd8dSchristos * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25ca32bd8dSchristos * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26ca32bd8dSchristos * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27ca32bd8dSchristos * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28ca32bd8dSchristos */
29ca32bd8dSchristos
30ca32bd8dSchristos /*
31ca32bd8dSchristos * Two-step process to generate safe primes for DHGEX
32ca32bd8dSchristos *
33ca32bd8dSchristos * Sieve candidates for "safe" primes,
34ca32bd8dSchristos * suitable for use as Diffie-Hellman moduli;
35ca32bd8dSchristos * that is, where q = (p-1)/2 is also prime.
36ca32bd8dSchristos *
37ca32bd8dSchristos * First step: generate candidate primes (memory intensive)
38ca32bd8dSchristos * Second step: test primes' safety (processor intensive)
39ca32bd8dSchristos */
40313c6c94Schristos #include "includes.h"
41*b1066cf3Schristos __RCSID("$NetBSD: moduli.c,v 1.17 2023/07/26 17:58:15 christos Exp $");
42ca32bd8dSchristos
43ca32bd8dSchristos #include <sys/types.h>
44ca32bd8dSchristos
45ca32bd8dSchristos #include <openssl/bn.h>
46ca32bd8dSchristos #include <openssl/dh.h>
47ca32bd8dSchristos
48091c4109Schristos #include <errno.h>
49ca32bd8dSchristos #include <stdio.h>
50ca32bd8dSchristos #include <stdlib.h>
51ca32bd8dSchristos #include <string.h>
52ca32bd8dSchristos #include <stdarg.h>
53ca32bd8dSchristos #include <time.h>
54091c4109Schristos #include <unistd.h>
55e4d43b82Schristos #include <limits.h>
56ca32bd8dSchristos
57ca32bd8dSchristos #include "xmalloc.h"
58ca32bd8dSchristos #include "dh.h"
59ca32bd8dSchristos #include "log.h"
608a4530f9Schristos #include "misc.h"
61ca32bd8dSchristos
62ca32bd8dSchristos /*
63ca32bd8dSchristos * File output defines
64ca32bd8dSchristos */
65ca32bd8dSchristos
66ca32bd8dSchristos /* need line long enough for largest moduli plus headers */
67ca32bd8dSchristos #define QLINESIZE (100+8192)
68ca32bd8dSchristos
69ca32bd8dSchristos /*
70ca32bd8dSchristos * Size: decimal.
71ca32bd8dSchristos * Specifies the number of the most significant bit (0 to M).
72ca32bd8dSchristos * WARNING: internally, usually 1 to N.
73ca32bd8dSchristos */
74ca32bd8dSchristos #define QSIZE_MINIMUM (511)
75ca32bd8dSchristos
76ca32bd8dSchristos /*
77ca32bd8dSchristos * Prime sieving defines
78ca32bd8dSchristos */
79ca32bd8dSchristos
80ca32bd8dSchristos /* Constant: assuming 8 bit bytes and 32 bit words */
81ca32bd8dSchristos #define SHIFT_BIT (3)
82ca32bd8dSchristos #define SHIFT_BYTE (2)
83ca32bd8dSchristos #define SHIFT_WORD (SHIFT_BIT+SHIFT_BYTE)
84ca32bd8dSchristos #define SHIFT_MEGABYTE (20)
85ca32bd8dSchristos #define SHIFT_MEGAWORD (SHIFT_MEGABYTE-SHIFT_BYTE)
86ca32bd8dSchristos
87ca32bd8dSchristos /*
88ca32bd8dSchristos * Using virtual memory can cause thrashing. This should be the largest
89ca32bd8dSchristos * number that is supported without a large amount of disk activity --
90ca32bd8dSchristos * that would increase the run time from hours to days or weeks!
91ca32bd8dSchristos */
92ca32bd8dSchristos #define LARGE_MINIMUM (8UL) /* megabytes */
93ca32bd8dSchristos
94ca32bd8dSchristos /*
95ca32bd8dSchristos * Do not increase this number beyond the unsigned integer bit size.
96ca32bd8dSchristos * Due to a multiple of 4, it must be LESS than 128 (yielding 2**30 bits).
97ca32bd8dSchristos */
98ca32bd8dSchristos #define LARGE_MAXIMUM (127UL) /* megabytes */
99ca32bd8dSchristos
100ca32bd8dSchristos /*
101ca32bd8dSchristos * Constant: when used with 32-bit integers, the largest sieve prime
102ca32bd8dSchristos * has to be less than 2**32.
103ca32bd8dSchristos */
104ca32bd8dSchristos #define SMALL_MAXIMUM (0xffffffffUL)
105ca32bd8dSchristos
106ca32bd8dSchristos /* Constant: can sieve all primes less than 2**32, as 65537**2 > 2**32-1. */
107ca32bd8dSchristos #define TINY_NUMBER (1UL<<16)
108ca32bd8dSchristos
109ca32bd8dSchristos /* Ensure enough bit space for testing 2*q. */
110ca32bd8dSchristos #define TEST_MAXIMUM (1UL<<16)
111ca32bd8dSchristos #define TEST_MINIMUM (QSIZE_MINIMUM + 1)
112ca32bd8dSchristos /* real TEST_MINIMUM (1UL << (SHIFT_WORD - TEST_POWER)) */
113ca32bd8dSchristos #define TEST_POWER (3) /* 2**n, n < SHIFT_WORD */
114ca32bd8dSchristos
115ca32bd8dSchristos /* bit operations on 32-bit words */
116ca32bd8dSchristos #define BIT_CLEAR(a,n) ((a)[(n)>>SHIFT_WORD] &= ~(1L << ((n) & 31)))
117ca32bd8dSchristos #define BIT_SET(a,n) ((a)[(n)>>SHIFT_WORD] |= (1L << ((n) & 31)))
118ca32bd8dSchristos #define BIT_TEST(a,n) ((a)[(n)>>SHIFT_WORD] & (1L << ((n) & 31)))
119ca32bd8dSchristos
120ca32bd8dSchristos /*
121ca32bd8dSchristos * Prime testing defines
122ca32bd8dSchristos */
123ca32bd8dSchristos
124ca32bd8dSchristos /* Minimum number of primality tests to perform */
125ca32bd8dSchristos #define TRIAL_MINIMUM (4)
126ca32bd8dSchristos
127ca32bd8dSchristos /*
128ca32bd8dSchristos * Sieving data (XXX - move to struct)
129ca32bd8dSchristos */
130ca32bd8dSchristos
131ca32bd8dSchristos /* sieve 2**16 */
132ca32bd8dSchristos static u_int32_t *TinySieve, tinybits;
133ca32bd8dSchristos
134ca32bd8dSchristos /* sieve 2**30 in 2**16 parts */
135ca32bd8dSchristos static u_int32_t *SmallSieve, smallbits, smallbase;
136ca32bd8dSchristos
137ca32bd8dSchristos /* sieve relative to the initial value */
138ca32bd8dSchristos static u_int32_t *LargeSieve, largewords, largetries, largenumbers;
139ca32bd8dSchristos static u_int32_t largebits, largememory; /* megabytes */
140ca32bd8dSchristos static BIGNUM *largebase;
141ca32bd8dSchristos
142ca32bd8dSchristos int gen_candidates(FILE *, u_int32_t, u_int32_t, BIGNUM *);
1432649c700Schristos int prime_test(FILE *, FILE *, u_int32_t, u_int32_t, char *, unsigned long,
1442649c700Schristos unsigned long);
145ca32bd8dSchristos
146ca32bd8dSchristos /*
147ca32bd8dSchristos * print moduli out in consistent form,
148ca32bd8dSchristos */
149ca32bd8dSchristos static int
qfileout(FILE * ofile,u_int32_t otype,u_int32_t otests,u_int32_t otries,u_int32_t osize,u_int32_t ogenerator,BIGNUM * omodulus)150ca32bd8dSchristos qfileout(FILE * ofile, u_int32_t otype, u_int32_t otests, u_int32_t otries,
151ca32bd8dSchristos u_int32_t osize, u_int32_t ogenerator, BIGNUM * omodulus)
152ca32bd8dSchristos {
153ca32bd8dSchristos struct tm *gtm;
154ca32bd8dSchristos time_t time_now;
155ca32bd8dSchristos int res;
156ca32bd8dSchristos
157ca32bd8dSchristos time(&time_now);
158ca32bd8dSchristos gtm = gmtime(&time_now);
159cd4ada6aSchristos if (gtm == NULL)
160cd4ada6aSchristos return -1;
161ca32bd8dSchristos
162ca32bd8dSchristos res = fprintf(ofile, "%04d%02d%02d%02d%02d%02d %u %u %u %u %x ",
163ca32bd8dSchristos gtm->tm_year + 1900, gtm->tm_mon + 1, gtm->tm_mday,
164ca32bd8dSchristos gtm->tm_hour, gtm->tm_min, gtm->tm_sec,
165ca32bd8dSchristos otype, otests, otries, osize, ogenerator);
166ca32bd8dSchristos
167ca32bd8dSchristos if (res < 0)
168ca32bd8dSchristos return (-1);
169ca32bd8dSchristos
170ca32bd8dSchristos if (BN_print_fp(ofile, omodulus) < 1)
171ca32bd8dSchristos return (-1);
172ca32bd8dSchristos
173ca32bd8dSchristos res = fprintf(ofile, "\n");
174ca32bd8dSchristos fflush(ofile);
175ca32bd8dSchristos
176ca32bd8dSchristos return (res > 0 ? 0 : -1);
177ca32bd8dSchristos }
178ca32bd8dSchristos
179ca32bd8dSchristos
180ca32bd8dSchristos /*
181ca32bd8dSchristos ** Sieve p's and q's with small factors
182ca32bd8dSchristos */
183ca32bd8dSchristos static void
sieve_large(u_int32_t s32)184e160b4e8Schristos sieve_large(u_int32_t s32)
185ca32bd8dSchristos {
186e160b4e8Schristos u_int64_t r, u, s = s32;
187ca32bd8dSchristos
188e160b4e8Schristos debug3("sieve_large %u", s32);
189ca32bd8dSchristos largetries++;
190ca32bd8dSchristos /* r = largebase mod s */
191e160b4e8Schristos r = BN_mod_word(largebase, s32);
192ca32bd8dSchristos if (r == 0)
193ca32bd8dSchristos u = 0; /* s divides into largebase exactly */
194ca32bd8dSchristos else
195ca32bd8dSchristos u = s - r; /* largebase+u is first entry divisible by s */
196ca32bd8dSchristos
197e160b4e8Schristos if (u < largebits * 2ULL) {
198ca32bd8dSchristos /*
199ca32bd8dSchristos * The sieve omits p's and q's divisible by 2, so ensure that
200ca32bd8dSchristos * largebase+u is odd. Then, step through the sieve in
201ca32bd8dSchristos * increments of 2*s
202ca32bd8dSchristos */
203ca32bd8dSchristos if (u & 0x1)
204ca32bd8dSchristos u += s; /* Make largebase+u odd, and u even */
205ca32bd8dSchristos
206ca32bd8dSchristos /* Mark all multiples of 2*s */
207ca32bd8dSchristos for (u /= 2; u < largebits; u += s)
208ca32bd8dSchristos BIT_SET(LargeSieve, u);
209ca32bd8dSchristos }
210ca32bd8dSchristos
211ca32bd8dSchristos /* r = p mod s */
212ca32bd8dSchristos r = (2 * r + 1) % s;
213ca32bd8dSchristos if (r == 0)
214ca32bd8dSchristos u = 0; /* s divides p exactly */
215ca32bd8dSchristos else
216ca32bd8dSchristos u = s - r; /* p+u is first entry divisible by s */
217ca32bd8dSchristos
218e160b4e8Schristos if (u < largebits * 4ULL) {
219ca32bd8dSchristos /*
220ca32bd8dSchristos * The sieve omits p's divisible by 4, so ensure that
221ca32bd8dSchristos * largebase+u is not. Then, step through the sieve in
222ca32bd8dSchristos * increments of 4*s
223ca32bd8dSchristos */
224ca32bd8dSchristos while (u & 0x3) {
225ca32bd8dSchristos if (SMALL_MAXIMUM - u < s)
226ca32bd8dSchristos return;
227ca32bd8dSchristos u += s;
228ca32bd8dSchristos }
229ca32bd8dSchristos
230ca32bd8dSchristos /* Mark all multiples of 4*s */
231ca32bd8dSchristos for (u /= 4; u < largebits; u += s)
232ca32bd8dSchristos BIT_SET(LargeSieve, u);
233ca32bd8dSchristos }
234ca32bd8dSchristos }
235ca32bd8dSchristos
236ca32bd8dSchristos /*
237ca32bd8dSchristos * list candidates for Sophie-Germain primes (where q = (p-1)/2)
238ca32bd8dSchristos * to standard output.
239ca32bd8dSchristos * The list is checked against small known primes (less than 2**30).
240ca32bd8dSchristos */
241ca32bd8dSchristos int
gen_candidates(FILE * out,u_int32_t memory,u_int32_t power,BIGNUM * start)242ca32bd8dSchristos gen_candidates(FILE *out, u_int32_t memory, u_int32_t power, BIGNUM *start)
243ca32bd8dSchristos {
244ca32bd8dSchristos BIGNUM *q;
245ca32bd8dSchristos u_int32_t j, r, s, t;
246ca32bd8dSchristos u_int32_t smallwords = TINY_NUMBER >> 6;
247ca32bd8dSchristos u_int32_t tinywords = TINY_NUMBER >> 6;
248ca32bd8dSchristos time_t time_start, time_stop;
249ca32bd8dSchristos u_int32_t i;
250ca32bd8dSchristos int ret = 0;
251ca32bd8dSchristos
252ca32bd8dSchristos largememory = memory;
253ca32bd8dSchristos
254ca32bd8dSchristos if (memory != 0 &&
255ca32bd8dSchristos (memory < LARGE_MINIMUM || memory > LARGE_MAXIMUM)) {
256ca32bd8dSchristos error("Invalid memory amount (min %ld, max %ld)",
257ca32bd8dSchristos LARGE_MINIMUM, LARGE_MAXIMUM);
258ca32bd8dSchristos return (-1);
259ca32bd8dSchristos }
260ca32bd8dSchristos
261ca32bd8dSchristos /*
262ca32bd8dSchristos * Set power to the length in bits of the prime to be generated.
263ca32bd8dSchristos * This is changed to 1 less than the desired safe prime moduli p.
264ca32bd8dSchristos */
265ca32bd8dSchristos if (power > TEST_MAXIMUM) {
266ca32bd8dSchristos error("Too many bits: %u > %lu", power, TEST_MAXIMUM);
267ca32bd8dSchristos return (-1);
268ca32bd8dSchristos } else if (power < TEST_MINIMUM) {
269ca32bd8dSchristos error("Too few bits: %u < %u", power, TEST_MINIMUM);
270ca32bd8dSchristos return (-1);
271ca32bd8dSchristos }
272ca32bd8dSchristos power--; /* decrement before squaring */
273ca32bd8dSchristos
274ca32bd8dSchristos /*
275ca32bd8dSchristos * The density of ordinary primes is on the order of 1/bits, so the
276ca32bd8dSchristos * density of safe primes should be about (1/bits)**2. Set test range
277ca32bd8dSchristos * to something well above bits**2 to be reasonably sure (but not
278ca32bd8dSchristos * guaranteed) of catching at least one safe prime.
279ca32bd8dSchristos */
280ca32bd8dSchristos largewords = ((power * power) >> (SHIFT_WORD - TEST_POWER));
281ca32bd8dSchristos
282ca32bd8dSchristos /*
283ca32bd8dSchristos * Need idea of how much memory is available. We don't have to use all
284ca32bd8dSchristos * of it.
285ca32bd8dSchristos */
286ca32bd8dSchristos if (largememory > LARGE_MAXIMUM) {
287ca32bd8dSchristos logit("Limited memory: %u MB; limit %lu MB",
288ca32bd8dSchristos largememory, LARGE_MAXIMUM);
289ca32bd8dSchristos largememory = LARGE_MAXIMUM;
290ca32bd8dSchristos }
291ca32bd8dSchristos
292ca32bd8dSchristos if (largewords <= (largememory << SHIFT_MEGAWORD)) {
293ca32bd8dSchristos logit("Increased memory: %u MB; need %u bytes",
294ca32bd8dSchristos largememory, (largewords << SHIFT_BYTE));
295ca32bd8dSchristos largewords = (largememory << SHIFT_MEGAWORD);
296ca32bd8dSchristos } else if (largememory > 0) {
297ca32bd8dSchristos logit("Decreased memory: %u MB; want %u bytes",
298ca32bd8dSchristos largememory, (largewords << SHIFT_BYTE));
299ca32bd8dSchristos largewords = (largememory << SHIFT_MEGAWORD);
300ca32bd8dSchristos }
301ca32bd8dSchristos
302ca32bd8dSchristos TinySieve = xcalloc(tinywords, sizeof(u_int32_t));
303ca32bd8dSchristos tinybits = tinywords << SHIFT_WORD;
304ca32bd8dSchristos
305ca32bd8dSchristos SmallSieve = xcalloc(smallwords, sizeof(u_int32_t));
306ca32bd8dSchristos smallbits = smallwords << SHIFT_WORD;
307ca32bd8dSchristos
308ca32bd8dSchristos /*
309ca32bd8dSchristos * dynamically determine available memory
310ca32bd8dSchristos */
311ca32bd8dSchristos while ((LargeSieve = calloc(largewords, sizeof(u_int32_t))) == NULL)
312ca32bd8dSchristos largewords -= (1L << (SHIFT_MEGAWORD - 2)); /* 1/4 MB chunks */
313ca32bd8dSchristos
314ca32bd8dSchristos largebits = largewords << SHIFT_WORD;
315ca32bd8dSchristos largenumbers = largebits * 2; /* even numbers excluded */
316ca32bd8dSchristos
317ca32bd8dSchristos /* validation check: count the number of primes tried */
318ca32bd8dSchristos largetries = 0;
319ca32bd8dSchristos if ((q = BN_new()) == NULL)
320ca32bd8dSchristos fatal("BN_new failed");
321ca32bd8dSchristos
322ca32bd8dSchristos /*
323ca32bd8dSchristos * Generate random starting point for subprime search, or use
324ca32bd8dSchristos * specified parameter.
325ca32bd8dSchristos */
326ca32bd8dSchristos if ((largebase = BN_new()) == NULL)
327ca32bd8dSchristos fatal("BN_new failed");
328ca32bd8dSchristos if (start == NULL) {
329ca32bd8dSchristos if (BN_rand(largebase, power, 1, 1) == 0)
330ca32bd8dSchristos fatal("BN_rand failed");
331ca32bd8dSchristos } else {
332ca32bd8dSchristos if (BN_copy(largebase, start) == NULL)
333ca32bd8dSchristos fatal("BN_copy: failed");
334ca32bd8dSchristos }
335ca32bd8dSchristos
336ca32bd8dSchristos /* ensure odd */
337ca32bd8dSchristos if (BN_set_bit(largebase, 0) == 0)
338ca32bd8dSchristos fatal("BN_set_bit: failed");
339ca32bd8dSchristos
340ca32bd8dSchristos time(&time_start);
341ca32bd8dSchristos
342ca32bd8dSchristos logit("%.24s Sieve next %u plus %u-bit", ctime(&time_start),
343ca32bd8dSchristos largenumbers, power);
344ca32bd8dSchristos debug2("start point: 0x%s", BN_bn2hex(largebase));
345ca32bd8dSchristos
346ca32bd8dSchristos /*
347ca32bd8dSchristos * TinySieve
348ca32bd8dSchristos */
349ca32bd8dSchristos for (i = 0; i < tinybits; i++) {
350ca32bd8dSchristos if (BIT_TEST(TinySieve, i))
351ca32bd8dSchristos continue; /* 2*i+3 is composite */
352ca32bd8dSchristos
353ca32bd8dSchristos /* The next tiny prime */
354ca32bd8dSchristos t = 2 * i + 3;
355ca32bd8dSchristos
356ca32bd8dSchristos /* Mark all multiples of t */
357ca32bd8dSchristos for (j = i + t; j < tinybits; j += t)
358ca32bd8dSchristos BIT_SET(TinySieve, j);
359ca32bd8dSchristos
360ca32bd8dSchristos sieve_large(t);
361ca32bd8dSchristos }
362ca32bd8dSchristos
363ca32bd8dSchristos /*
364ca32bd8dSchristos * Start the small block search at the next possible prime. To avoid
365ca32bd8dSchristos * fencepost errors, the last pass is skipped.
366ca32bd8dSchristos */
367ca32bd8dSchristos for (smallbase = TINY_NUMBER + 3;
368ca32bd8dSchristos smallbase < (SMALL_MAXIMUM - TINY_NUMBER);
369ca32bd8dSchristos smallbase += TINY_NUMBER) {
370ca32bd8dSchristos for (i = 0; i < tinybits; i++) {
371ca32bd8dSchristos if (BIT_TEST(TinySieve, i))
372ca32bd8dSchristos continue; /* 2*i+3 is composite */
373ca32bd8dSchristos
374ca32bd8dSchristos /* The next tiny prime */
375ca32bd8dSchristos t = 2 * i + 3;
376ca32bd8dSchristos r = smallbase % t;
377ca32bd8dSchristos
378ca32bd8dSchristos if (r == 0) {
379ca32bd8dSchristos s = 0; /* t divides into smallbase exactly */
380ca32bd8dSchristos } else {
381ca32bd8dSchristos /* smallbase+s is first entry divisible by t */
382ca32bd8dSchristos s = t - r;
383ca32bd8dSchristos }
384ca32bd8dSchristos
385ca32bd8dSchristos /*
386ca32bd8dSchristos * The sieve omits even numbers, so ensure that
387ca32bd8dSchristos * smallbase+s is odd. Then, step through the sieve
388ca32bd8dSchristos * in increments of 2*t
389ca32bd8dSchristos */
390ca32bd8dSchristos if (s & 1)
391ca32bd8dSchristos s += t; /* Make smallbase+s odd, and s even */
392ca32bd8dSchristos
393ca32bd8dSchristos /* Mark all multiples of 2*t */
394ca32bd8dSchristos for (s /= 2; s < smallbits; s += t)
395ca32bd8dSchristos BIT_SET(SmallSieve, s);
396ca32bd8dSchristos }
397ca32bd8dSchristos
398ca32bd8dSchristos /*
399ca32bd8dSchristos * SmallSieve
400ca32bd8dSchristos */
401ca32bd8dSchristos for (i = 0; i < smallbits; i++) {
402ca32bd8dSchristos if (BIT_TEST(SmallSieve, i))
403ca32bd8dSchristos continue; /* 2*i+smallbase is composite */
404ca32bd8dSchristos
405ca32bd8dSchristos /* The next small prime */
406ca32bd8dSchristos sieve_large((2 * i) + smallbase);
407ca32bd8dSchristos }
408ca32bd8dSchristos
409ca32bd8dSchristos memset(SmallSieve, 0, smallwords << SHIFT_BYTE);
410ca32bd8dSchristos }
411ca32bd8dSchristos
412ca32bd8dSchristos time(&time_stop);
413ca32bd8dSchristos
414ffae97bbSchristos logit("%.24s Sieved with %u small primes in %lld seconds",
415ffae97bbSchristos ctime(&time_stop), largetries, (long long)(time_stop - time_start));
416ca32bd8dSchristos
417ca32bd8dSchristos for (j = r = 0; j < largebits; j++) {
418ca32bd8dSchristos if (BIT_TEST(LargeSieve, j))
419ca32bd8dSchristos continue; /* Definitely composite, skip */
420ca32bd8dSchristos
421ca32bd8dSchristos debug2("test q = largebase+%u", 2 * j);
422ca32bd8dSchristos if (BN_set_word(q, 2 * j) == 0)
423ca32bd8dSchristos fatal("BN_set_word failed");
424ca32bd8dSchristos if (BN_add(q, q, largebase) == 0)
425ca32bd8dSchristos fatal("BN_add failed");
426ca32bd8dSchristos if (qfileout(out, MODULI_TYPE_SOPHIE_GERMAIN,
427ca32bd8dSchristos MODULI_TESTS_SIEVE, largetries,
428ca32bd8dSchristos (power - 1) /* MSB */, (0), q) == -1) {
429ca32bd8dSchristos ret = -1;
430ca32bd8dSchristos break;
431ca32bd8dSchristos }
432ca32bd8dSchristos
433ca32bd8dSchristos r++; /* count q */
434ca32bd8dSchristos }
435ca32bd8dSchristos
436ca32bd8dSchristos time(&time_stop);
437ca32bd8dSchristos
43800a838c4Schristos free(LargeSieve);
43900a838c4Schristos free(SmallSieve);
44000a838c4Schristos free(TinySieve);
441ca32bd8dSchristos
442ca32bd8dSchristos logit("%.24s Found %u candidates", ctime(&time_stop), r);
443ca32bd8dSchristos
444ca32bd8dSchristos return (ret);
445ca32bd8dSchristos }
446ca32bd8dSchristos
447091c4109Schristos static void
write_checkpoint(char * cpfile,u_int32_t lineno)448091c4109Schristos write_checkpoint(char *cpfile, u_int32_t lineno)
449091c4109Schristos {
450091c4109Schristos FILE *fp;
451e4d43b82Schristos char tmp[PATH_MAX];
452*b1066cf3Schristos int r, writeok, closeok;
453091c4109Schristos
454091c4109Schristos r = snprintf(tmp, sizeof(tmp), "%s.XXXXXXXXXX", cpfile);
455cd4ada6aSchristos if (r < 0 || r >= PATH_MAX) {
456091c4109Schristos logit("write_checkpoint: temp pathname too long");
457091c4109Schristos return;
458091c4109Schristos }
459091c4109Schristos if ((r = mkstemp(tmp)) == -1) {
460091c4109Schristos logit("mkstemp(%s): %s", tmp, strerror(errno));
461091c4109Schristos return;
462091c4109Schristos }
463091c4109Schristos if ((fp = fdopen(r, "w")) == NULL) {
464091c4109Schristos logit("write_checkpoint: fdopen: %s", strerror(errno));
465e4d43b82Schristos unlink(tmp);
466091c4109Schristos close(r);
467091c4109Schristos return;
468091c4109Schristos }
469*b1066cf3Schristos writeok = (fprintf(fp, "%lu\n", (unsigned long)lineno) > 0);
470*b1066cf3Schristos closeok = (fclose(fp) == 0);
471*b1066cf3Schristos if (writeok && closeok && rename(tmp, cpfile) == 0) {
472091c4109Schristos debug3("wrote checkpoint line %lu to '%s'",
473091c4109Schristos (unsigned long)lineno, cpfile);
474*b1066cf3Schristos } else {
475091c4109Schristos logit("failed to write to checkpoint file '%s': %s", cpfile,
476091c4109Schristos strerror(errno));
477*b1066cf3Schristos (void)unlink(tmp);
478*b1066cf3Schristos }
479091c4109Schristos }
480091c4109Schristos
481091c4109Schristos static unsigned long
read_checkpoint(char * cpfile)482091c4109Schristos read_checkpoint(char *cpfile)
483091c4109Schristos {
484091c4109Schristos FILE *fp;
485091c4109Schristos unsigned long lineno = 0;
486091c4109Schristos
487091c4109Schristos if ((fp = fopen(cpfile, "r")) == NULL)
488091c4109Schristos return 0;
489091c4109Schristos if (fscanf(fp, "%lu\n", &lineno) < 1)
490091c4109Schristos logit("Failed to load checkpoint from '%s'", cpfile);
491091c4109Schristos else
492091c4109Schristos logit("Loaded checkpoint from '%s' line %lu", cpfile, lineno);
493091c4109Schristos fclose(fp);
494091c4109Schristos return lineno;
495091c4109Schristos }
496091c4109Schristos
4978a4530f9Schristos static unsigned long
count_lines(FILE * f)4988a4530f9Schristos count_lines(FILE *f)
4998a4530f9Schristos {
5008a4530f9Schristos unsigned long count = 0;
5018a4530f9Schristos char lp[QLINESIZE + 1];
5028a4530f9Schristos
5038a4530f9Schristos if (fseek(f, 0, SEEK_SET) != 0) {
5048a4530f9Schristos debug("input file is not seekable");
5058a4530f9Schristos return ULONG_MAX;
5068a4530f9Schristos }
5078a4530f9Schristos while (fgets(lp, QLINESIZE + 1, f) != NULL)
5088a4530f9Schristos count++;
5098a4530f9Schristos rewind(f);
5108a4530f9Schristos debug("input file has %lu lines", count);
5118a4530f9Schristos return count;
5128a4530f9Schristos }
5138a4530f9Schristos
5148a4530f9Schristos static char *
fmt_time(time_t seconds)5158a4530f9Schristos fmt_time(time_t seconds)
5168a4530f9Schristos {
5178a4530f9Schristos int day, hr, min;
5188a4530f9Schristos static char buf[128];
5198a4530f9Schristos
5208a4530f9Schristos min = (seconds / 60) % 60;
5218a4530f9Schristos hr = (seconds / 60 / 60) % 24;
5228a4530f9Schristos day = seconds / 60 / 60 / 24;
5238a4530f9Schristos if (day > 0)
5248a4530f9Schristos snprintf(buf, sizeof buf, "%dd %d:%02d", day, hr, min);
5258a4530f9Schristos else
5268a4530f9Schristos snprintf(buf, sizeof buf, "%d:%02d", hr, min);
5278a4530f9Schristos return buf;
5288a4530f9Schristos }
5298a4530f9Schristos
5308a4530f9Schristos static void
print_progress(unsigned long start_lineno,unsigned long current_lineno,unsigned long end_lineno)5318a4530f9Schristos print_progress(unsigned long start_lineno, unsigned long current_lineno,
5328a4530f9Schristos unsigned long end_lineno)
5338a4530f9Schristos {
5348a4530f9Schristos static time_t time_start, time_prev;
5358a4530f9Schristos time_t time_now, elapsed;
5368a4530f9Schristos unsigned long num_to_process, processed, remaining, percent, eta;
5378a4530f9Schristos double time_per_line;
5388a4530f9Schristos char *eta_str;
5398a4530f9Schristos
5408a4530f9Schristos time_now = monotime();
5418a4530f9Schristos if (time_start == 0) {
5428a4530f9Schristos time_start = time_prev = time_now;
5438a4530f9Schristos return;
5448a4530f9Schristos }
5458a4530f9Schristos /* print progress after 1m then once per 5m */
5468a4530f9Schristos if (time_now - time_prev < 5 * 60)
5478a4530f9Schristos return;
5488a4530f9Schristos time_prev = time_now;
5498a4530f9Schristos elapsed = time_now - time_start;
5508a4530f9Schristos processed = current_lineno - start_lineno;
5518a4530f9Schristos remaining = end_lineno - current_lineno;
5528a4530f9Schristos num_to_process = end_lineno - start_lineno;
5538a4530f9Schristos time_per_line = (double)elapsed / processed;
5548a4530f9Schristos /* if we don't know how many we're processing just report count+time */
5558a4530f9Schristos time(&time_now);
5568a4530f9Schristos if (end_lineno == ULONG_MAX) {
5578a4530f9Schristos logit("%.24s processed %lu in %s", ctime(&time_now),
5588a4530f9Schristos processed, fmt_time(elapsed));
5598a4530f9Schristos return;
5608a4530f9Schristos }
5618a4530f9Schristos percent = 100 * processed / num_to_process;
5628a4530f9Schristos eta = time_per_line * remaining;
5638a4530f9Schristos eta_str = xstrdup(fmt_time(eta));
5648a4530f9Schristos logit("%.24s processed %lu of %lu (%lu%%) in %s, ETA %s",
5658a4530f9Schristos ctime(&time_now), processed, num_to_process, percent,
5668a4530f9Schristos fmt_time(elapsed), eta_str);
5678a4530f9Schristos free(eta_str);
5688a4530f9Schristos }
5698a4530f9Schristos
570ca32bd8dSchristos /*
571ca32bd8dSchristos * perform a Miller-Rabin primality test
572ca32bd8dSchristos * on the list of candidates
573ca32bd8dSchristos * (checking both q and p)
574ca32bd8dSchristos * The result is a list of so-call "safe" primes
575ca32bd8dSchristos */
576ca32bd8dSchristos int
prime_test(FILE * in,FILE * out,u_int32_t trials,u_int32_t generator_wanted,char * checkpoint_file,unsigned long start_lineno,unsigned long num_lines)577091c4109Schristos prime_test(FILE *in, FILE *out, u_int32_t trials, u_int32_t generator_wanted,
5782649c700Schristos char *checkpoint_file, unsigned long start_lineno, unsigned long num_lines)
579ca32bd8dSchristos {
580ca32bd8dSchristos BIGNUM *q, *p, *a;
581ca32bd8dSchristos char *cp, *lp;
582ca32bd8dSchristos u_int32_t count_in = 0, count_out = 0, count_possible = 0;
583ca32bd8dSchristos u_int32_t generator_known, in_tests, in_tries, in_type, in_size;
5842649c700Schristos unsigned long last_processed = 0, end_lineno;
585ca32bd8dSchristos time_t time_start, time_stop;
586aa36fcacSchristos int res, is_prime;
587ca32bd8dSchristos
588ca32bd8dSchristos if (trials < TRIAL_MINIMUM) {
589ca32bd8dSchristos error("Minimum primality trials is %d", TRIAL_MINIMUM);
590ca32bd8dSchristos return (-1);
591ca32bd8dSchristos }
592ca32bd8dSchristos
5938a4530f9Schristos if (num_lines == 0)
5948a4530f9Schristos end_lineno = count_lines(in);
5958a4530f9Schristos else
5968a4530f9Schristos end_lineno = start_lineno + num_lines;
5978a4530f9Schristos
598ca32bd8dSchristos time(&time_start);
599ca32bd8dSchristos
600ca32bd8dSchristos if ((p = BN_new()) == NULL)
601ca32bd8dSchristos fatal("BN_new failed");
602ca32bd8dSchristos if ((q = BN_new()) == NULL)
603ca32bd8dSchristos fatal("BN_new failed");
604ca32bd8dSchristos
605ca32bd8dSchristos debug2("%.24s Final %u Miller-Rabin trials (%x generator)",
606ca32bd8dSchristos ctime(&time_start), trials, generator_wanted);
607ca32bd8dSchristos
608091c4109Schristos if (checkpoint_file != NULL)
609091c4109Schristos last_processed = read_checkpoint(checkpoint_file);
610ee85abc4Schristos last_processed = start_lineno = MAXIMUM(last_processed, start_lineno);
6118a4530f9Schristos if (end_lineno == ULONG_MAX)
6128a4530f9Schristos debug("process from line %lu from pipe", last_processed);
6132649c700Schristos else
6148a4530f9Schristos debug("process from line %lu to line %lu", last_processed,
6158a4530f9Schristos end_lineno);
616091c4109Schristos
617ca32bd8dSchristos res = 0;
618ca32bd8dSchristos lp = xmalloc(QLINESIZE + 1);
6192649c700Schristos while (fgets(lp, QLINESIZE + 1, in) != NULL && count_in < end_lineno) {
620ca32bd8dSchristos count_in++;
621091c4109Schristos if (count_in <= last_processed) {
6228a4530f9Schristos debug3("skipping line %u, before checkpoint or "
6238a4530f9Schristos "specified start line", count_in);
624091c4109Schristos continue;
625091c4109Schristos }
6268a4530f9Schristos if (checkpoint_file != NULL)
627091c4109Schristos write_checkpoint(checkpoint_file, count_in);
6288a4530f9Schristos print_progress(start_lineno, count_in, end_lineno);
629ca32bd8dSchristos if (strlen(lp) < 14 || *lp == '!' || *lp == '#') {
630ca32bd8dSchristos debug2("%10u: comment or short line", count_in);
631ca32bd8dSchristos continue;
632ca32bd8dSchristos }
633ca32bd8dSchristos
634ca32bd8dSchristos /* XXX - fragile parser */
635ca32bd8dSchristos /* time */
636ca32bd8dSchristos cp = &lp[14]; /* (skip) */
637ca32bd8dSchristos
638ca32bd8dSchristos /* type */
639ca32bd8dSchristos in_type = strtoul(cp, &cp, 10);
640ca32bd8dSchristos
641ca32bd8dSchristos /* tests */
642ca32bd8dSchristos in_tests = strtoul(cp, &cp, 10);
643ca32bd8dSchristos
644ca32bd8dSchristos if (in_tests & MODULI_TESTS_COMPOSITE) {
645ca32bd8dSchristos debug2("%10u: known composite", count_in);
646ca32bd8dSchristos continue;
647ca32bd8dSchristos }
648ca32bd8dSchristos
649ca32bd8dSchristos /* tries */
650ca32bd8dSchristos in_tries = strtoul(cp, &cp, 10);
651ca32bd8dSchristos
652ca32bd8dSchristos /* size (most significant bit) */
653ca32bd8dSchristos in_size = strtoul(cp, &cp, 10);
654ca32bd8dSchristos
655ca32bd8dSchristos /* generator (hex) */
656ca32bd8dSchristos generator_known = strtoul(cp, &cp, 16);
657ca32bd8dSchristos
658ca32bd8dSchristos /* Skip white space */
659ca32bd8dSchristos cp += strspn(cp, " ");
660ca32bd8dSchristos
661ca32bd8dSchristos /* modulus (hex) */
662ca32bd8dSchristos switch (in_type) {
663ca32bd8dSchristos case MODULI_TYPE_SOPHIE_GERMAIN:
664ca32bd8dSchristos debug2("%10u: (%u) Sophie-Germain", count_in, in_type);
665ca32bd8dSchristos a = q;
666ca32bd8dSchristos if (BN_hex2bn(&a, cp) == 0)
667ca32bd8dSchristos fatal("BN_hex2bn failed");
668ca32bd8dSchristos /* p = 2*q + 1 */
669ca32bd8dSchristos if (BN_lshift(p, q, 1) == 0)
670ca32bd8dSchristos fatal("BN_lshift failed");
671ca32bd8dSchristos if (BN_add_word(p, 1) == 0)
672ca32bd8dSchristos fatal("BN_add_word failed");
673ca32bd8dSchristos in_size += 1;
674ca32bd8dSchristos generator_known = 0;
675ca32bd8dSchristos break;
676ca32bd8dSchristos case MODULI_TYPE_UNSTRUCTURED:
677ca32bd8dSchristos case MODULI_TYPE_SAFE:
678ca32bd8dSchristos case MODULI_TYPE_SCHNORR:
679ca32bd8dSchristos case MODULI_TYPE_STRONG:
680ca32bd8dSchristos case MODULI_TYPE_UNKNOWN:
681ca32bd8dSchristos debug2("%10u: (%u)", count_in, in_type);
682ca32bd8dSchristos a = p;
683ca32bd8dSchristos if (BN_hex2bn(&a, cp) == 0)
684ca32bd8dSchristos fatal("BN_hex2bn failed");
685ca32bd8dSchristos /* q = (p-1) / 2 */
686ca32bd8dSchristos if (BN_rshift(q, p, 1) == 0)
687ca32bd8dSchristos fatal("BN_rshift failed");
688ca32bd8dSchristos break;
689ca32bd8dSchristos default:
690ca32bd8dSchristos debug2("Unknown prime type");
691ca32bd8dSchristos break;
692ca32bd8dSchristos }
693ca32bd8dSchristos
694ca32bd8dSchristos /*
695ca32bd8dSchristos * due to earlier inconsistencies in interpretation, check
696ca32bd8dSchristos * the proposed bit size.
697ca32bd8dSchristos */
698ca32bd8dSchristos if ((u_int32_t)BN_num_bits(p) != (in_size + 1)) {
699ca32bd8dSchristos debug2("%10u: bit size %u mismatch", count_in, in_size);
700ca32bd8dSchristos continue;
701ca32bd8dSchristos }
702ca32bd8dSchristos if (in_size < QSIZE_MINIMUM) {
703ca32bd8dSchristos debug2("%10u: bit size %u too short", count_in, in_size);
704ca32bd8dSchristos continue;
705ca32bd8dSchristos }
706ca32bd8dSchristos
707ca32bd8dSchristos if (in_tests & MODULI_TESTS_MILLER_RABIN)
708ca32bd8dSchristos in_tries += trials;
709ca32bd8dSchristos else
710ca32bd8dSchristos in_tries = trials;
711ca32bd8dSchristos
712ca32bd8dSchristos /*
713ca32bd8dSchristos * guess unknown generator
714ca32bd8dSchristos */
715ca32bd8dSchristos if (generator_known == 0) {
716ca32bd8dSchristos if (BN_mod_word(p, 24) == 11)
717ca32bd8dSchristos generator_known = 2;
718ca32bd8dSchristos else {
719ca32bd8dSchristos u_int32_t r = BN_mod_word(p, 10);
720ca32bd8dSchristos
721ca32bd8dSchristos if (r == 3 || r == 7)
722ca32bd8dSchristos generator_known = 5;
723ca32bd8dSchristos }
724ca32bd8dSchristos }
725ca32bd8dSchristos /*
726ca32bd8dSchristos * skip tests when desired generator doesn't match
727ca32bd8dSchristos */
728ca32bd8dSchristos if (generator_wanted > 0 &&
729ca32bd8dSchristos generator_wanted != generator_known) {
730ca32bd8dSchristos debug2("%10u: generator %d != %d",
731ca32bd8dSchristos count_in, generator_known, generator_wanted);
732ca32bd8dSchristos continue;
733ca32bd8dSchristos }
734ca32bd8dSchristos
735ca32bd8dSchristos /*
736ca32bd8dSchristos * Primes with no known generator are useless for DH, so
737ca32bd8dSchristos * skip those.
738ca32bd8dSchristos */
739ca32bd8dSchristos if (generator_known == 0) {
740ca32bd8dSchristos debug2("%10u: no known generator", count_in);
741ca32bd8dSchristos continue;
742ca32bd8dSchristos }
743ca32bd8dSchristos
744ca32bd8dSchristos count_possible++;
745ca32bd8dSchristos
746ca32bd8dSchristos /*
747ca32bd8dSchristos * The (1/4)^N performance bound on Miller-Rabin is
748ca32bd8dSchristos * extremely pessimistic, so don't spend a lot of time
749ca32bd8dSchristos * really verifying that q is prime until after we know
750ca32bd8dSchristos * that p is also prime. A single pass will weed out the
751ca32bd8dSchristos * vast majority of composite q's.
752ca32bd8dSchristos */
753ed75d7a8Schristos is_prime = BN_is_prime_ex(q, 1, NULL, NULL);
754aa36fcacSchristos if (is_prime < 0)
755aa36fcacSchristos fatal("BN_is_prime_ex failed");
756aa36fcacSchristos if (is_prime == 0) {
757ca32bd8dSchristos debug("%10u: q failed first possible prime test",
758ca32bd8dSchristos count_in);
759ca32bd8dSchristos continue;
760ca32bd8dSchristos }
761ca32bd8dSchristos
762ca32bd8dSchristos /*
763ca32bd8dSchristos * q is possibly prime, so go ahead and really make sure
764ca32bd8dSchristos * that p is prime. If it is, then we can go back and do
765ca32bd8dSchristos * the same for q. If p is composite, chances are that
766ca32bd8dSchristos * will show up on the first Rabin-Miller iteration so it
767ca32bd8dSchristos * doesn't hurt to specify a high iteration count.
768ca32bd8dSchristos */
769ed75d7a8Schristos is_prime = BN_is_prime_ex(p, trials, NULL, NULL);
770aa36fcacSchristos if (is_prime < 0)
771aa36fcacSchristos fatal("BN_is_prime_ex failed");
772aa36fcacSchristos if (is_prime == 0) {
773ca32bd8dSchristos debug("%10u: p is not prime", count_in);
774ca32bd8dSchristos continue;
775ca32bd8dSchristos }
776ca32bd8dSchristos debug("%10u: p is almost certainly prime", count_in);
777ca32bd8dSchristos
778ca32bd8dSchristos /* recheck q more rigorously */
779ed75d7a8Schristos is_prime = BN_is_prime_ex(q, trials - 1, NULL, NULL);
780aa36fcacSchristos if (is_prime < 0)
781aa36fcacSchristos fatal("BN_is_prime_ex failed");
782aa36fcacSchristos if (is_prime == 0) {
783ca32bd8dSchristos debug("%10u: q is not prime", count_in);
784ca32bd8dSchristos continue;
785ca32bd8dSchristos }
786ca32bd8dSchristos debug("%10u: q is almost certainly prime", count_in);
787ca32bd8dSchristos
788ca32bd8dSchristos if (qfileout(out, MODULI_TYPE_SAFE,
789ca32bd8dSchristos in_tests | MODULI_TESTS_MILLER_RABIN,
790ca32bd8dSchristos in_tries, in_size, generator_known, p)) {
791ca32bd8dSchristos res = -1;
792ca32bd8dSchristos break;
793ca32bd8dSchristos }
794ca32bd8dSchristos
795ca32bd8dSchristos count_out++;
796ca32bd8dSchristos }
797ca32bd8dSchristos
798ca32bd8dSchristos time(&time_stop);
79900a838c4Schristos free(lp);
800ca32bd8dSchristos BN_free(p);
801ca32bd8dSchristos BN_free(q);
802ca32bd8dSchristos
803091c4109Schristos if (checkpoint_file != NULL)
804091c4109Schristos unlink(checkpoint_file);
805091c4109Schristos
806ca32bd8dSchristos logit("%.24s Found %u safe primes of %u candidates in %ld seconds",
807ca32bd8dSchristos ctime(&time_stop), count_out, count_possible,
808ca32bd8dSchristos (long) (time_stop - time_start));
809ca32bd8dSchristos
810ca32bd8dSchristos return (res);
811ca32bd8dSchristos }
812