Add trunc and truncf.
[dragonfly.git] / crypto / openssh-4 / moduli.c
1 /* $OpenBSD: moduli.c,v 1.20 2007/02/24 03:30:11 ray Exp $ */
2 /*
3  * Copyright 1994 Phil Karn <karn@qualcomm.com>
4  * Copyright 1996-1998, 2003 William Allen Simpson <wsimpson@greendragon.com>
5  * Copyright 2000 Niels Provos <provos@citi.umich.edu>
6  * All rights reserved.
7  *
8  * Redistribution and use in source and binary forms, with or without
9  * modification, are permitted provided that the following conditions
10  * are met:
11  * 1. Redistributions of source code must retain the above copyright
12  *    notice, this list of conditions and the following disclaimer.
13  * 2. Redistributions in binary form must reproduce the above copyright
14  *    notice, this list of conditions and the following disclaimer in the
15  *    documentation and/or other materials provided with the distribution.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
28
29 /*
30  * Two-step process to generate safe primes for DHGEX
31  *
32  *  Sieve candidates for "safe" primes,
33  *  suitable for use as Diffie-Hellman moduli;
34  *  that is, where q = (p-1)/2 is also prime.
35  *
36  * First step: generate candidate primes (memory intensive)
37  * Second step: test primes' safety (processor intensive)
38  */
39
40 #include "includes.h"
41
42 #include <sys/types.h>
43
44 #include <openssl/bn.h>
45
46 #include <stdio.h>
47 #include <stdlib.h>
48 #include <string.h>
49 #include <stdarg.h>
50 #include <time.h>
51
52 #include "xmalloc.h"
53 #include "log.h"
54
55 /*
56  * File output defines
57  */
58
59 /* need line long enough for largest moduli plus headers */
60 #define QLINESIZE               (100+8192)
61
62 /* Type: decimal.
63  * Specifies the internal structure of the prime modulus.
64  */
65 #define QTYPE_UNKNOWN           (0)
66 #define QTYPE_UNSTRUCTURED      (1)
67 #define QTYPE_SAFE              (2)
68 #define QTYPE_SCHNORR           (3)
69 #define QTYPE_SOPHIE_GERMAIN    (4)
70 #define QTYPE_STRONG            (5)
71
72 /* Tests: decimal (bit field).
73  * Specifies the methods used in checking for primality.
74  * Usually, more than one test is used.
75  */
76 #define QTEST_UNTESTED          (0x00)
77 #define QTEST_COMPOSITE         (0x01)
78 #define QTEST_SIEVE             (0x02)
79 #define QTEST_MILLER_RABIN      (0x04)
80 #define QTEST_JACOBI            (0x08)
81 #define QTEST_ELLIPTIC          (0x10)
82
83 /*
84  * Size: decimal.
85  * Specifies the number of the most significant bit (0 to M).
86  * WARNING: internally, usually 1 to N.
87  */
88 #define QSIZE_MINIMUM           (511)
89
90 /*
91  * Prime sieving defines
92  */
93
94 /* Constant: assuming 8 bit bytes and 32 bit words */
95 #define SHIFT_BIT       (3)
96 #define SHIFT_BYTE      (2)
97 #define SHIFT_WORD      (SHIFT_BIT+SHIFT_BYTE)
98 #define SHIFT_MEGABYTE  (20)
99 #define SHIFT_MEGAWORD  (SHIFT_MEGABYTE-SHIFT_BYTE)
100
101 /*
102  * Using virtual memory can cause thrashing.  This should be the largest
103  * number that is supported without a large amount of disk activity --
104  * that would increase the run time from hours to days or weeks!
105  */
106 #define LARGE_MINIMUM   (8UL)   /* megabytes */
107
108 /*
109  * Do not increase this number beyond the unsigned integer bit size.
110  * Due to a multiple of 4, it must be LESS than 128 (yielding 2**30 bits).
111  */
112 #define LARGE_MAXIMUM   (127UL) /* megabytes */
113
114 /*
115  * Constant: when used with 32-bit integers, the largest sieve prime
116  * has to be less than 2**32.
117  */
118 #define SMALL_MAXIMUM   (0xffffffffUL)
119
120 /* Constant: can sieve all primes less than 2**32, as 65537**2 > 2**32-1. */
121 #define TINY_NUMBER     (1UL<<16)
122
123 /* Ensure enough bit space for testing 2*q. */
124 #define TEST_MAXIMUM    (1UL<<16)
125 #define TEST_MINIMUM    (QSIZE_MINIMUM + 1)
126 /* real TEST_MINIMUM    (1UL << (SHIFT_WORD - TEST_POWER)) */
127 #define TEST_POWER      (3)     /* 2**n, n < SHIFT_WORD */
128
129 /* bit operations on 32-bit words */
130 #define BIT_CLEAR(a,n)  ((a)[(n)>>SHIFT_WORD] &= ~(1L << ((n) & 31)))
131 #define BIT_SET(a,n)    ((a)[(n)>>SHIFT_WORD] |= (1L << ((n) & 31)))
132 #define BIT_TEST(a,n)   ((a)[(n)>>SHIFT_WORD] & (1L << ((n) & 31)))
133
134 /*
135  * Prime testing defines
136  */
137
138 /* Minimum number of primality tests to perform */
139 #define TRIAL_MINIMUM   (4)
140
141 /*
142  * Sieving data (XXX - move to struct)
143  */
144
145 /* sieve 2**16 */
146 static u_int32_t *TinySieve, tinybits;
147
148 /* sieve 2**30 in 2**16 parts */
149 static u_int32_t *SmallSieve, smallbits, smallbase;
150
151 /* sieve relative to the initial value */
152 static u_int32_t *LargeSieve, largewords, largetries, largenumbers;
153 static u_int32_t largebits, largememory;        /* megabytes */
154 static BIGNUM *largebase;
155
156 int gen_candidates(FILE *, u_int32_t, u_int32_t, BIGNUM *);
157 int prime_test(FILE *, FILE *, u_int32_t, u_int32_t);
158
159 /*
160  * print moduli out in consistent form,
161  */
162 static int
163 qfileout(FILE * ofile, u_int32_t otype, u_int32_t otests, u_int32_t otries,
164     u_int32_t osize, u_int32_t ogenerator, BIGNUM * omodulus)
165 {
166         struct tm *gtm;
167         time_t time_now;
168         int res;
169
170         time(&time_now);
171         gtm = gmtime(&time_now);
172
173         res = fprintf(ofile, "%04d%02d%02d%02d%02d%02d %u %u %u %u %x ",
174             gtm->tm_year + 1900, gtm->tm_mon + 1, gtm->tm_mday,
175             gtm->tm_hour, gtm->tm_min, gtm->tm_sec,
176             otype, otests, otries, osize, ogenerator);
177
178         if (res < 0)
179                 return (-1);
180
181         if (BN_print_fp(ofile, omodulus) < 1)
182                 return (-1);
183
184         res = fprintf(ofile, "\n");
185         fflush(ofile);
186
187         return (res > 0 ? 0 : -1);
188 }
189
190
191 /*
192  ** Sieve p's and q's with small factors
193  */
194 static void
195 sieve_large(u_int32_t s)
196 {
197         u_int32_t r, u;
198
199         debug3("sieve_large %u", s);
200         largetries++;
201         /* r = largebase mod s */
202         r = BN_mod_word(largebase, s);
203         if (r == 0)
204                 u = 0; /* s divides into largebase exactly */
205         else
206                 u = s - r; /* largebase+u is first entry divisible by s */
207
208         if (u < largebits * 2) {
209                 /*
210                  * The sieve omits p's and q's divisible by 2, so ensure that
211                  * largebase+u is odd. Then, step through the sieve in
212                  * increments of 2*s
213                  */
214                 if (u & 0x1)
215                         u += s; /* Make largebase+u odd, and u even */
216
217                 /* Mark all multiples of 2*s */
218                 for (u /= 2; u < largebits; u += s)
219                         BIT_SET(LargeSieve, u);
220         }
221
222         /* r = p mod s */
223         r = (2 * r + 1) % s;
224         if (r == 0)
225                 u = 0; /* s divides p exactly */
226         else
227                 u = s - r; /* p+u is first entry divisible by s */
228
229         if (u < largebits * 4) {
230                 /*
231                  * The sieve omits p's divisible by 4, so ensure that
232                  * largebase+u is not. Then, step through the sieve in
233                  * increments of 4*s
234                  */
235                 while (u & 0x3) {
236                         if (SMALL_MAXIMUM - u < s)
237                                 return;
238                         u += s;
239                 }
240
241                 /* Mark all multiples of 4*s */
242                 for (u /= 4; u < largebits; u += s)
243                         BIT_SET(LargeSieve, u);
244         }
245 }
246
247 /*
248  * list candidates for Sophie-Germain primes (where q = (p-1)/2)
249  * to standard output.
250  * The list is checked against small known primes (less than 2**30).
251  */
252 int
253 gen_candidates(FILE *out, u_int32_t memory, u_int32_t power, BIGNUM *start)
254 {
255         BIGNUM *q;
256         u_int32_t j, r, s, t;
257         u_int32_t smallwords = TINY_NUMBER >> 6;
258         u_int32_t tinywords = TINY_NUMBER >> 6;
259         time_t time_start, time_stop;
260         u_int32_t i;
261         int ret = 0;
262
263         largememory = memory;
264
265         if (memory != 0 &&
266             (memory < LARGE_MINIMUM || memory > LARGE_MAXIMUM)) {
267                 error("Invalid memory amount (min %ld, max %ld)",
268                     LARGE_MINIMUM, LARGE_MAXIMUM);
269                 return (-1);
270         }
271
272         /*
273          * Set power to the length in bits of the prime to be generated.
274          * This is changed to 1 less than the desired safe prime moduli p.
275          */
276         if (power > TEST_MAXIMUM) {
277                 error("Too many bits: %u > %lu", power, TEST_MAXIMUM);
278                 return (-1);
279         } else if (power < TEST_MINIMUM) {
280                 error("Too few bits: %u < %u", power, TEST_MINIMUM);
281                 return (-1);
282         }
283         power--; /* decrement before squaring */
284
285         /*
286          * The density of ordinary primes is on the order of 1/bits, so the
287          * density of safe primes should be about (1/bits)**2. Set test range
288          * to something well above bits**2 to be reasonably sure (but not
289          * guaranteed) of catching at least one safe prime.
290          */
291         largewords = ((power * power) >> (SHIFT_WORD - TEST_POWER));
292
293         /*
294          * Need idea of how much memory is available. We don't have to use all
295          * of it.
296          */
297         if (largememory > LARGE_MAXIMUM) {
298                 logit("Limited memory: %u MB; limit %lu MB",
299                     largememory, LARGE_MAXIMUM);
300                 largememory = LARGE_MAXIMUM;
301         }
302
303         if (largewords <= (largememory << SHIFT_MEGAWORD)) {
304                 logit("Increased memory: %u MB; need %u bytes",
305                     largememory, (largewords << SHIFT_BYTE));
306                 largewords = (largememory << SHIFT_MEGAWORD);
307         } else if (largememory > 0) {
308                 logit("Decreased memory: %u MB; want %u bytes",
309                     largememory, (largewords << SHIFT_BYTE));
310                 largewords = (largememory << SHIFT_MEGAWORD);
311         }
312
313         TinySieve = xcalloc(tinywords, sizeof(u_int32_t));
314         tinybits = tinywords << SHIFT_WORD;
315
316         SmallSieve = xcalloc(smallwords, sizeof(u_int32_t));
317         smallbits = smallwords << SHIFT_WORD;
318
319         /*
320          * dynamically determine available memory
321          */
322         while ((LargeSieve = calloc(largewords, sizeof(u_int32_t))) == NULL)
323                 largewords -= (1L << (SHIFT_MEGAWORD - 2)); /* 1/4 MB chunks */
324
325         largebits = largewords << SHIFT_WORD;
326         largenumbers = largebits * 2;   /* even numbers excluded */
327
328         /* validation check: count the number of primes tried */
329         largetries = 0;
330         if ((q = BN_new()) == NULL)
331                 fatal("BN_new failed");
332
333         /*
334          * Generate random starting point for subprime search, or use
335          * specified parameter.
336          */
337         if ((largebase = BN_new()) == NULL)
338                 fatal("BN_new failed");
339         if (start == NULL) {
340                 if (BN_rand(largebase, power, 1, 1) == 0)
341                         fatal("BN_rand failed");
342         } else {
343                 if (BN_copy(largebase, start) == NULL)
344                         fatal("BN_copy: failed");
345         }
346
347         /* ensure odd */
348         if (BN_set_bit(largebase, 0) == 0)
349                 fatal("BN_set_bit: failed");
350
351         time(&time_start);
352
353         logit("%.24s Sieve next %u plus %u-bit", ctime(&time_start),
354             largenumbers, power);
355         debug2("start point: 0x%s", BN_bn2hex(largebase));
356
357         /*
358          * TinySieve
359          */
360         for (i = 0; i < tinybits; i++) {
361                 if (BIT_TEST(TinySieve, i))
362                         continue; /* 2*i+3 is composite */
363
364                 /* The next tiny prime */
365                 t = 2 * i + 3;
366
367                 /* Mark all multiples of t */
368                 for (j = i + t; j < tinybits; j += t)
369                         BIT_SET(TinySieve, j);
370
371                 sieve_large(t);
372         }
373
374         /*
375          * Start the small block search at the next possible prime. To avoid
376          * fencepost errors, the last pass is skipped.
377          */
378         for (smallbase = TINY_NUMBER + 3;
379             smallbase < (SMALL_MAXIMUM - TINY_NUMBER);
380             smallbase += TINY_NUMBER) {
381                 for (i = 0; i < tinybits; i++) {
382                         if (BIT_TEST(TinySieve, i))
383                                 continue; /* 2*i+3 is composite */
384
385                         /* The next tiny prime */
386                         t = 2 * i + 3;
387                         r = smallbase % t;
388
389                         if (r == 0) {
390                                 s = 0; /* t divides into smallbase exactly */
391                         } else {
392                                 /* smallbase+s is first entry divisible by t */
393                                 s = t - r;
394                         }
395
396                         /*
397                          * The sieve omits even numbers, so ensure that
398                          * smallbase+s is odd. Then, step through the sieve
399                          * in increments of 2*t
400                          */
401                         if (s & 1)
402                                 s += t; /* Make smallbase+s odd, and s even */
403
404                         /* Mark all multiples of 2*t */
405                         for (s /= 2; s < smallbits; s += t)
406                                 BIT_SET(SmallSieve, s);
407                 }
408
409                 /*
410                  * SmallSieve
411                  */
412                 for (i = 0; i < smallbits; i++) {
413                         if (BIT_TEST(SmallSieve, i))
414                                 continue; /* 2*i+smallbase is composite */
415
416                         /* The next small prime */
417                         sieve_large((2 * i) + smallbase);
418                 }
419
420                 memset(SmallSieve, 0, smallwords << SHIFT_BYTE);
421         }
422
423         time(&time_stop);
424
425         logit("%.24s Sieved with %u small primes in %ld seconds",
426             ctime(&time_stop), largetries, (long) (time_stop - time_start));
427
428         for (j = r = 0; j < largebits; j++) {
429                 if (BIT_TEST(LargeSieve, j))
430                         continue; /* Definitely composite, skip */
431
432                 debug2("test q = largebase+%u", 2 * j);
433                 if (BN_set_word(q, 2 * j) == 0)
434                         fatal("BN_set_word failed");
435                 if (BN_add(q, q, largebase) == 0)
436                         fatal("BN_add failed");
437                 if (qfileout(out, QTYPE_SOPHIE_GERMAIN, QTEST_SIEVE,
438                     largetries, (power - 1) /* MSB */, (0), q) == -1) {
439                         ret = -1;
440                         break;
441                 }
442
443                 r++; /* count q */
444         }
445
446         time(&time_stop);
447
448         xfree(LargeSieve);
449         xfree(SmallSieve);
450         xfree(TinySieve);
451
452         logit("%.24s Found %u candidates", ctime(&time_stop), r);
453
454         return (ret);
455 }
456
457 /*
458  * perform a Miller-Rabin primality test
459  * on the list of candidates
460  * (checking both q and p)
461  * The result is a list of so-call "safe" primes
462  */
463 int
464 prime_test(FILE *in, FILE *out, u_int32_t trials, u_int32_t generator_wanted)
465 {
466         BIGNUM *q, *p, *a;
467         BN_CTX *ctx;
468         char *cp, *lp;
469         u_int32_t count_in = 0, count_out = 0, count_possible = 0;
470         u_int32_t generator_known, in_tests, in_tries, in_type, in_size;
471         time_t time_start, time_stop;
472         int res;
473
474         if (trials < TRIAL_MINIMUM) {
475                 error("Minimum primality trials is %d", TRIAL_MINIMUM);
476                 return (-1);
477         }
478
479         time(&time_start);
480
481         if ((p = BN_new()) == NULL)
482                 fatal("BN_new failed");
483         if ((q = BN_new()) == NULL)
484                 fatal("BN_new failed");
485         if ((ctx = BN_CTX_new()) == NULL)
486                 fatal("BN_CTX_new failed");
487
488         debug2("%.24s Final %u Miller-Rabin trials (%x generator)",
489             ctime(&time_start), trials, generator_wanted);
490
491         res = 0;
492         lp = xmalloc(QLINESIZE + 1);
493         while (fgets(lp, QLINESIZE + 1, in) != NULL) {
494                 count_in++;
495                 if (strlen(lp) < 14 || *lp == '!' || *lp == '#') {
496                         debug2("%10u: comment or short line", count_in);
497                         continue;
498                 }
499
500                 /* XXX - fragile parser */
501                 /* time */
502                 cp = &lp[14];   /* (skip) */
503
504                 /* type */
505                 in_type = strtoul(cp, &cp, 10);
506
507                 /* tests */
508                 in_tests = strtoul(cp, &cp, 10);
509
510                 if (in_tests & QTEST_COMPOSITE) {
511                         debug2("%10u: known composite", count_in);
512                         continue;
513                 }
514
515                 /* tries */
516                 in_tries = strtoul(cp, &cp, 10);
517
518                 /* size (most significant bit) */
519                 in_size = strtoul(cp, &cp, 10);
520
521                 /* generator (hex) */
522                 generator_known = strtoul(cp, &cp, 16);
523
524                 /* Skip white space */
525                 cp += strspn(cp, " ");
526
527                 /* modulus (hex) */
528                 switch (in_type) {
529                 case QTYPE_SOPHIE_GERMAIN:
530                         debug2("%10u: (%u) Sophie-Germain", count_in, in_type);
531                         a = q;
532                         if (BN_hex2bn(&a, cp) == 0)
533                                 fatal("BN_hex2bn failed");
534                         /* p = 2*q + 1 */
535                         if (BN_lshift(p, q, 1) == 0)
536                                 fatal("BN_lshift failed");
537                         if (BN_add_word(p, 1) == 0)
538                                 fatal("BN_add_word failed");
539                         in_size += 1;
540                         generator_known = 0;
541                         break;
542                 case QTYPE_UNSTRUCTURED:
543                 case QTYPE_SAFE:
544                 case QTYPE_SCHNORR:
545                 case QTYPE_STRONG:
546                 case QTYPE_UNKNOWN:
547                         debug2("%10u: (%u)", count_in, in_type);
548                         a = p;
549                         if (BN_hex2bn(&a, cp) == 0)
550                                 fatal("BN_hex2bn failed");
551                         /* q = (p-1) / 2 */
552                         if (BN_rshift(q, p, 1) == 0)
553                                 fatal("BN_rshift failed");
554                         break;
555                 default:
556                         debug2("Unknown prime type");
557                         break;
558                 }
559
560                 /*
561                  * due to earlier inconsistencies in interpretation, check
562                  * the proposed bit size.
563                  */
564                 if ((u_int32_t)BN_num_bits(p) != (in_size + 1)) {
565                         debug2("%10u: bit size %u mismatch", count_in, in_size);
566                         continue;
567                 }
568                 if (in_size < QSIZE_MINIMUM) {
569                         debug2("%10u: bit size %u too short", count_in, in_size);
570                         continue;
571                 }
572
573                 if (in_tests & QTEST_MILLER_RABIN)
574                         in_tries += trials;
575                 else
576                         in_tries = trials;
577
578                 /*
579                  * guess unknown generator
580                  */
581                 if (generator_known == 0) {
582                         if (BN_mod_word(p, 24) == 11)
583                                 generator_known = 2;
584                         else if (BN_mod_word(p, 12) == 5)
585                                 generator_known = 3;
586                         else {
587                                 u_int32_t r = BN_mod_word(p, 10);
588
589                                 if (r == 3 || r == 7)
590                                         generator_known = 5;
591                         }
592                 }
593                 /*
594                  * skip tests when desired generator doesn't match
595                  */
596                 if (generator_wanted > 0 &&
597                     generator_wanted != generator_known) {
598                         debug2("%10u: generator %d != %d",
599                             count_in, generator_known, generator_wanted);
600                         continue;
601                 }
602
603                 /*
604                  * Primes with no known generator are useless for DH, so
605                  * skip those.
606                  */
607                 if (generator_known == 0) {
608                         debug2("%10u: no known generator", count_in);
609                         continue;
610                 }
611
612                 count_possible++;
613
614                 /*
615                  * The (1/4)^N performance bound on Miller-Rabin is
616                  * extremely pessimistic, so don't spend a lot of time
617                  * really verifying that q is prime until after we know
618                  * that p is also prime. A single pass will weed out the
619                  * vast majority of composite q's.
620                  */
621                 if (BN_is_prime(q, 1, NULL, ctx, NULL) <= 0) {
622                         debug("%10u: q failed first possible prime test",
623                             count_in);
624                         continue;
625                 }
626
627                 /*
628                  * q is possibly prime, so go ahead and really make sure
629                  * that p is prime. If it is, then we can go back and do
630                  * the same for q. If p is composite, chances are that
631                  * will show up on the first Rabin-Miller iteration so it
632                  * doesn't hurt to specify a high iteration count.
633                  */
634                 if (!BN_is_prime(p, trials, NULL, ctx, NULL)) {
635                         debug("%10u: p is not prime", count_in);
636                         continue;
637                 }
638                 debug("%10u: p is almost certainly prime", count_in);
639
640                 /* recheck q more rigorously */
641                 if (!BN_is_prime(q, trials - 1, NULL, ctx, NULL)) {
642                         debug("%10u: q is not prime", count_in);
643                         continue;
644                 }
645                 debug("%10u: q is almost certainly prime", count_in);
646
647                 if (qfileout(out, QTYPE_SAFE, (in_tests | QTEST_MILLER_RABIN),
648                     in_tries, in_size, generator_known, p)) {
649                         res = -1;
650                         break;
651                 }
652
653                 count_out++;
654         }
655
656         time(&time_stop);
657         xfree(lp);
658         BN_free(p);
659         BN_free(q);
660         BN_CTX_free(ctx);
661
662         logit("%.24s Found %u safe primes of %u candidates in %ld seconds",
663             ctime(&time_stop), count_out, count_possible,
664             (long) (time_stop - time_start));
665
666         return (res);
667 }