chiark / gitweb /
math/strongprime.c: Muffle an irritating style warning from Clang.
[catacomb] / math / strongprime.c
1 /* -*-c-*-
2  *
3  * Generate `strong' prime numbers
4  *
5  * (c) 1999 Straylight/Edgeware
6  */
7
8 /*----- Licensing notice --------------------------------------------------*
9  *
10  * This file is part of Catacomb.
11  *
12  * Catacomb is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU Library General Public License as
14  * published by the Free Software Foundation; either version 2 of the
15  * License, or (at your option) any later version.
16  *
17  * Catacomb is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20  * GNU Library General Public License for more details.
21  *
22  * You should have received a copy of the GNU Library General Public
23  * License along with Catacomb; if not, write to the Free
24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25  * MA 02111-1307, USA.
26  */
27
28 /*----- Header files ------------------------------------------------------*/
29
30 #include <mLib/dstr.h>
31 #include <mLib/macros.h>
32
33 #include "grand.h"
34 #include "mp.h"
35 #include "mpmont.h"
36 #include "mprand.h"
37 #include "pgen.h"
38 #include "pfilt.h"
39 #include "rabin.h"
40
41 /*----- Main code ---------------------------------------------------------*/
42
43 /* Oh, just shut up. */
44 CLANG_WARNING("-Wempty-body")
45
46 /* --- @strongprime_setup@ --- *
47  *
48  * Arguments:   @const char *name@ = pointer to name root
49  *              @mp *d@ = destination for search start point
50  *              @pfilt *f@ = where to store filter jump context
51  *              @unsigned nbits@ = number of bits wanted
52  *              @grand *r@ = random number source
53  *              @unsigned n@ = number of attempts to make
54  *              @pgen_proc *event@ = event handler function
55  *              @void *ectx@ = argument for the event handler
56  *
57  * Returns:     A starting point for a `strong' prime search, or zero.
58  *
59  * Use:         Sets up for a strong prime search, so that primes with
60  *              particular properties can be found.  It's probably important
61  *              to note that the number left in the filter context @f@ is
62  *              congruent to 2 (mod 4); that the jump value is twice the
63  *              product of two large primes; and that the starting point is
64  *              at least %$3 \cdot 2^{N-2}$%.  (Hence, if you multiply two
65  *              such numbers, the product is at least
66  *
67  *                      %$9 \cdot 2^{2N-4} > 2^{2N-1}$%
68  *
69  *              i.e., it will be (at least) a %$2 N$%-bit value.
70  */
71
72 mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
73                       grand *r, unsigned n, pgen_proc *event, void *ectx)
74 {
75   mp *s, *t, *q;
76   dstr dn = DSTR_INIT;
77   unsigned slop, nb, u, i;
78
79   mp *rr = d;
80   pgen_filterctx c;
81   pgen_jumpctx j;
82   rabin rb;
83
84   /* --- Figure out how large the smaller primes should be --- *
85    *
86    * We want them to be `as large as possible', subject to the constraint
87    * that we produce a number of the requested size at the end.  This is
88    * tricky, because the final prime search is going to involve quite large
89    * jumps from its starting point; the size of the jumps are basically
90    * determined by our choice here, and if they're too big then we won't find
91    * a prime in time.
92    *
93    * Let's suppose we're trying to make an %$N$%-bit prime.  The expected
94    * number of steps tends to increase linearly with size, i.e., we need to
95    * take about %2^k N$% steps for some %$k$%.  If we're jumping by a
96    * %$J$%-bit quantity each time, from an %$N$%-bit starting point, then we
97    * will only be able to find a match if %$2^k N 2^{J-1} \le 2^{N-1}$%,
98    * i.e., if %$J \le N - (k + \log_2 N)$%.
99    *
100    * Experimentation shows that taking %$k + \log_2 N = 12$% works well for
101    * %$N = 1024$%, so %$k = 2$%.  Add a few extra bits for luck.
102    */
103
104   for (i = 1; i && nbits >> i; i <<= 1); assert(i);
105   for (slop = 6, nb = nbits; nb > 1; i >>= 1) {
106     u = nb >> i;
107     if (u) { slop += i; nb = u; }
108   }
109   if (nbits/2 <= slop) return (0);
110
111   /* --- Choose two primes %$s$% and %$t$% of half the required size --- */
112
113   nb = nbits/2 - slop;
114   c.step = 1;
115
116   rr = mprand(rr, nb, r, 1);
117   DRESET(&dn); dstr_putf(&dn, "%s [s]", name);
118   if ((s = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
119                 rabin_iters(nb), pgen_test, &rb)) == 0)
120     goto fail_s;
121
122   rr = mprand(rr, nb, r, 1);
123   DRESET(&dn); dstr_putf(&dn, "%s [t]", name);
124   if ((t = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
125                 rabin_iters(nb), pgen_test, &rb)) == 0)
126     goto fail_t;
127
128   /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- *
129    *
130    * Then %$r \equiv 1 \pmod{t}$%, i.e., %$r - 1$% is a multiple of %$t$%.
131    */
132
133   rr = mp_lsl(rr, t, 1);
134   pfilt_create(&c.f, rr);
135   rr = mp_lsl(rr, rr, slop - 1);
136   rr = mp_add(rr, rr, MP_ONE);
137   DRESET(&dn); dstr_putf(&dn, "%s [r]", name);
138   j.j = &c.f;
139   q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j,
140            rabin_iters(nb + slop), pgen_test, &rb);
141   pfilt_destroy(&c.f);
142   if (!q)
143     goto fail_r;
144
145   /* --- Select a suitable congruence class for %$p$% --- *
146    *
147    * This computes %$p_0 = 2 s (s^{-1} \bmod r) - 1$%.  Then %$p_0 + 1$% is
148    * clearly a multiple of %$s$%, and
149    *
150    *    %$p_0 - 1 \equiv 2 s s^{-1} - 2 \equiv 0 \pmod{r}$%
151    *
152    * is a multiple of %$r$%.
153    */
154
155   rr = mp_modinv(rr, s, q);
156   rr = mp_mul(rr, rr, s);
157   rr = mp_lsl(rr, rr, 1);
158   rr = mp_sub(rr, rr, MP_ONE);
159
160   /* --- Pick a starting point for the search --- *
161    *
162    * Select %$3 \cdot 2^{N-2} < p_1 < 2^N$% at random, only with
163    * %$p_1 \equiv p_0 \pmod{2 r s}$.
164    */
165
166   {
167     mp *x, *y;
168     x = mp_mul(MP_NEW, q, s);
169     x = mp_lsl(x, x, 1);
170     pfilt_create(f, x); /* %$2 r s$% */
171     y = mprand(MP_NEW, nbits, r, 0);
172     y = mp_setbit(y, y, nbits - 2);
173     rr = mp_leastcongruent(rr, y, rr, x);
174     mp_drop(x); mp_drop(y);
175   }
176
177   /* --- Return the result --- */
178
179   mp_drop(q);
180   mp_drop(t);
181   mp_drop(s);
182   dstr_destroy(&dn);
183   return (rr);
184
185   /* --- Tidy up if something failed --- */
186
187 fail_r:
188   mp_drop(t);
189 fail_t:
190   mp_drop(s);
191 fail_s:
192   mp_drop(rr);
193   dstr_destroy(&dn);
194   return (0);
195 }
196
197 /* --- @strongprime@ --- *
198  *
199  * Arguments:   @const char *name@ = pointer to name root
200  *              @mp *d@ = destination integer
201  *              @unsigned nbits@ = number of bits wanted
202  *              @grand *r@ = random number source
203  *              @unsigned n@ = number of attempts to make
204  *              @pgen_proc *event@ = event handler function
205  *              @void *ectx@ = argument for the event handler
206  *
207  * Returns:     A `strong' prime, or zero.
208  *
209  * Use:         Finds `strong' primes.  A strong prime %$p$% is such that
210  *
211  *                * %$p - 1$% has a large prime factor %$r$%,
212  *                * %$p + 1$% has a large prime factor %$s$%, and
213  *                * %$r - 1$% has a large prime factor %$t$%.
214  */
215
216 mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
217                 unsigned n, pgen_proc *event, void *ectx)
218 {
219   mp *p;
220   pfilt f;
221   pgen_jumpctx j;
222   rabin rb;
223
224   if (d) mp_copy(d);
225   p = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
226   if (!p) { mp_drop(d); return (0); }
227   j.j = &f;
228   p = pgen(name, p, p, event, ectx, n, pgen_jump, &j,
229            rabin_iters(nbits), pgen_test, &rb);
230   if (mp_bits(p) != nbits) { mp_drop(p); return (0); }
231   pfilt_destroy(&f);
232   mp_drop(d);
233   return (p);
234 }
235
236 /*----- That's all, folks -------------------------------------------------*/