3 * Generate `strong' prime numbers
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
30 #include <mLib/dstr.h>
31 #include <mLib/macros.h>
41 /*----- Main code ---------------------------------------------------------*/
43 /* Oh, just shut up. */
44 CLANG_WARNING("-Wempty-body")
46 /* --- @strongprime_setup@ --- *
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
57 * Returns: A starting point for a `strong' prime search, or zero.
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
67 * %$9 \cdot 2^{2N-4} > 2^{2N-1}$%
69 * i.e., it will be (at least) a %$2 N$%-bit value.
72 mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
73 grand *r, unsigned n, pgen_proc *event, void *ectx)
77 unsigned slop, nb, u, i;
84 /* --- Figure out how large the smaller primes should be --- *
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
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)$%.
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.
104 for (i = 1; i && nbits >> i; i <<= 1); assert(i);
105 for (slop = 6, nb = nbits; nb > 1; i >>= 1) {
107 if (u) { slop += i; nb = u; }
109 if (nbits/2 <= slop) return (0);
111 /* --- Choose two primes %$s$% and %$t$% of half the required size --- */
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)
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)
128 /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- *
130 * Then %$r \equiv 1 \pmod{t}$%, i.e., %$r - 1$% is a multiple of %$t$%.
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);
139 q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j,
140 rabin_iters(nb + slop), pgen_test, &rb);
145 /* --- Select a suitable congruence class for %$p$% --- *
147 * This computes %$p_0 = 2 s (s^{-1} \bmod r) - 1$%. Then %$p_0 + 1$% is
148 * clearly a multiple of %$s$%, and
150 * %$p_0 - 1 \equiv 2 s s^{-1} - 2 \equiv 0 \pmod{r}$%
152 * is a multiple of %$r$%.
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);
160 /* --- Pick a starting point for the search --- *
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}$.
168 x = mp_mul(MP_NEW, q, s);
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);
177 /* --- Return the result --- */
185 /* --- Tidy up if something failed --- */
197 /* --- @strongprime@ --- *
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
207 * Returns: A `strong' prime, or zero.
209 * Use: Finds `strong' primes. A strong prime %$p$% is such that
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$%.
216 mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
217 unsigned n, pgen_proc *event, void *ectx)
225 p = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
226 if (!p) { mp_drop(d); return (0); }
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); }
236 /*----- That's all, folks -------------------------------------------------*/