9a8b0c8d |
1 | /* -*-c-*- |
2 | * |
34e4f738 |
3 | * $Id: pfilt.c,v 1.5 2004/04/01 12:50:09 mdw Exp $ |
9a8b0c8d |
4 | * |
5 | * Finding and testing prime numbers |
6 | * |
7 | * (c) 1999 Straylight/Edgeware |
8 | */ |
9 | |
10 | /*----- Licensing notice --------------------------------------------------* |
11 | * |
12 | * This file is part of Catacomb. |
13 | * |
14 | * Catacomb is free software; you can redistribute it and/or modify |
15 | * it under the terms of the GNU Library General Public License as |
16 | * published by the Free Software Foundation; either version 2 of the |
17 | * License, or (at your option) any later version. |
18 | * |
19 | * Catacomb is distributed in the hope that it will be useful, |
20 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
21 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
22 | * GNU Library General Public License for more details. |
23 | * |
24 | * You should have received a copy of the GNU Library General Public |
25 | * License along with Catacomb; if not, write to the Free |
26 | * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
27 | * MA 02111-1307, USA. |
28 | */ |
29 | |
30 | /*----- Revision history --------------------------------------------------* |
31 | * |
32 | * $Log: pfilt.c,v $ |
34e4f738 |
33 | * Revision 1.5 2004/04/01 12:50:09 mdw |
34 | * Add cyclic group abstraction, with test code. Separate off exponentation |
35 | * functions for better static linking. Fix a buttload of bugs on the way. |
36 | * Generally ensure that negative exponents do inversion correctly. Add |
37 | * table of standard prime-field subgroups. (Binary field subgroups are |
38 | * currently unimplemented but easy to add if anyone ever finds a good one.) |
39 | * |
9312c71f |
40 | * Revision 1.4 2000/10/08 12:14:57 mdw |
41 | * Remove vestiges of @primorial@. |
42 | * |
b8e79eeb |
43 | * Revision 1.3 2000/08/15 21:44:27 mdw |
44 | * (pfilt_smallfactor): New function for doing trial division the hard |
45 | * way. |
46 | * |
47 | * (pfilt_create): Use @mpx_udivn@ for computing residues, for improved |
48 | * performance. |
49 | * |
50 | * Pull the `small prime' test into a separate function, and do it |
51 | * properly. |
52 | * |
87b63d41 |
53 | * Revision 1.2 2000/06/17 11:54:27 mdw |
54 | * Use new MP memory management functions. |
55 | * |
9a8b0c8d |
56 | * Revision 1.1 1999/12/22 15:49:39 mdw |
57 | * Renamed from `pgen'. Reworking for new prime-search system. |
58 | * |
59 | * Revision 1.3 1999/12/10 23:28:35 mdw |
60 | * Track suggested destination changes. |
61 | * |
62 | * Revision 1.2 1999/11/20 22:23:05 mdw |
63 | * Add multiply-and-add function for Diffie-Hellman safe prime generation. |
64 | * |
65 | * Revision 1.1 1999/11/19 13:17:57 mdw |
66 | * Prime number generator and tester. |
67 | * |
68 | */ |
69 | |
70 | /*----- Header files ------------------------------------------------------*/ |
71 | |
72 | #include "mp.h" |
b8e79eeb |
73 | #include "mpint.h" |
9a8b0c8d |
74 | #include "pfilt.h" |
75 | #include "pgen.h" |
76 | #include "primetab.h" |
77 | |
78 | /*----- Main code ---------------------------------------------------------*/ |
79 | |
b8e79eeb |
80 | /* --- @smallenough@ --- * |
81 | * |
82 | * Arguments: @mp *m@ = integer to test |
83 | * |
84 | * Returns: One of the @PGEN@ result codes. |
85 | * |
86 | * Use: Assuming that @m@ has been tested by trial division on every |
87 | * prime in the small-primes array, this function will return |
88 | * @PGEN_DONE@ if the number is less than the square of the |
89 | * largest small prime. |
90 | */ |
91 | |
92 | static int smallenough(mp *m) |
93 | { |
94 | static mp *max = 0; |
95 | int rc = PGEN_TRY; |
96 | |
97 | if (!max) { |
98 | max = mp_fromuint(MP_NEW, MAXPRIME); |
99 | max = mp_sqr(max, max); |
100 | max->a->n--; /* Permanent allocation */ |
101 | } |
102 | if (MP_CMP(m, <, max)) |
103 | rc = PGEN_DONE; |
104 | return (rc); |
105 | } |
106 | |
107 | /* --- @pfilt_smallfactor@ --- * |
108 | * |
109 | * Arguments: @mp *m@ = integer to test |
110 | * |
111 | * Returns: One of the @PGEN@ result codes. |
112 | * |
113 | * Use: Tests a number by dividing by a number of small primes. This |
114 | * is a useful first step if you're testing random primes; for |
115 | * sequential searches, @pfilt_create@ works better. |
116 | */ |
117 | |
118 | int pfilt_smallfactor(mp *m) |
119 | { |
120 | int rc = PGEN_TRY; |
121 | int i; |
122 | size_t sz = MP_LEN(m); |
34e4f738 |
123 | mparena *a = m->a ? m->a : MPARENA_GLOBAL; |
124 | mpw *v = mpalloc(a, sz); |
b8e79eeb |
125 | |
126 | /* --- Fill in the residues --- */ |
127 | |
128 | for (i = 0; i < NPRIME; i++) { |
129 | if (!mpx_udivn(v, v + sz, m->v, m->vl, primetab[i])) { |
130 | if (MP_LEN(m) == 1 && m->v[0] == primetab[i]) |
131 | rc = PGEN_DONE; |
132 | else |
133 | rc = PGEN_FAIL; |
134 | } |
135 | } |
136 | |
137 | /* --- Check for small primes --- */ |
138 | |
139 | if (rc == PGEN_TRY) |
140 | rc = smallenough(m); |
141 | |
142 | /* --- Done --- */ |
143 | |
34e4f738 |
144 | mpfree(a, v); |
b8e79eeb |
145 | return (rc); |
146 | } |
147 | |
9a8b0c8d |
148 | /* --- @pfilt_create@ --- * |
149 | * |
150 | * Arguments: @pfilt *p@ = pointer to prime filtering context |
151 | * @mp *m@ = pointer to initial number to test |
152 | * |
153 | * Returns: One of the @PGEN@ result codes. |
154 | * |
155 | * Use: Tests an initial number for primality by computing its |
156 | * residue modulo various small prime numbers. This is fairly |
157 | * quick, but not particularly certain. If a @PGEN_TRY@ |
158 | * result is returned, perform Rabin-Miller tests to confirm. |
159 | */ |
160 | |
161 | int pfilt_create(pfilt *p, mp *m) |
162 | { |
163 | int rc = PGEN_TRY; |
164 | int i; |
b8e79eeb |
165 | size_t sz = MP_LEN(m); |
34e4f738 |
166 | mparena *a = m->a ? m->a : MPARENA_GLOBAL; |
167 | mpw *v = mpalloc(a, sz); |
9a8b0c8d |
168 | |
169 | /* --- Take a copy of the number --- */ |
170 | |
171 | mp_shrink(m); |
172 | p->m = MP_COPY(m); |
173 | |
174 | /* --- Fill in the residues --- */ |
175 | |
9a8b0c8d |
176 | for (i = 0; i < NPRIME; i++) { |
b8e79eeb |
177 | p->r[i] = mpx_udivn(v, v + sz, m->v, m->vl, primetab[i]); |
9a8b0c8d |
178 | if (!p->r[i] && rc == PGEN_TRY) { |
179 | if (MP_LEN(m) == 1 && m->v[0] == primetab[i]) |
180 | rc = PGEN_DONE; |
181 | else |
182 | rc = PGEN_FAIL; |
183 | } |
184 | } |
185 | |
b8e79eeb |
186 | /* --- Check for small primes --- */ |
187 | |
188 | if (rc == PGEN_TRY) |
189 | rc = smallenough(m); |
190 | |
9a8b0c8d |
191 | /* --- Done --- */ |
192 | |
34e4f738 |
193 | mpfree(a, v); |
9a8b0c8d |
194 | return (rc); |
195 | } |
196 | |
197 | /* --- @pfilt_destroy@ --- * |
198 | * |
199 | * Arguments: @pfilt *p@ = pointer to prime filtering context |
200 | * |
201 | * Returns: --- |
202 | * |
203 | * Use: Discards a context and all the resources it holds. |
204 | */ |
205 | |
206 | void pfilt_destroy(pfilt *p) |
207 | { |
208 | mp_drop(p->m); |
209 | } |
210 | |
211 | /* --- @pfilt_step@ --- * |
212 | * |
213 | * Arguments: @pfilt *p@ = pointer to prime filtering context |
214 | * @mpw step@ = how much to step the number |
215 | * |
216 | * Returns: One of the @PGEN@ result codes. |
217 | * |
218 | * Use: Steps a number by a small amount. Stepping is much faster |
219 | * than initializing with a new number. The test performed is |
220 | * the same simple one used by @primetab_create@, so @PGEN_TRY@ |
221 | * results should be followed up by a Rabin-Miller test. |
222 | */ |
223 | |
224 | int pfilt_step(pfilt *p, mpw step) |
225 | { |
226 | int rc = PGEN_TRY; |
227 | int i; |
228 | |
229 | /* --- Add the step on to the number --- */ |
230 | |
231 | p->m = mp_split(p->m); |
232 | mp_ensure(p->m, MP_LEN(p->m) + 1); |
233 | mpx_uaddn(p->m->v, p->m->vl, step); |
234 | mp_shrink(p->m); |
235 | |
236 | /* --- Update the residue table --- */ |
237 | |
238 | for (i = 0; i < NPRIME; i++) { |
239 | p->r[i] = (p->r[i] + step) % primetab[i]; |
240 | if (!p->r[i] && rc == PGEN_TRY) { |
241 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) |
242 | rc = PGEN_DONE; |
243 | else |
244 | rc = PGEN_FAIL; |
245 | } |
246 | } |
247 | |
b8e79eeb |
248 | /* --- Check for small primes --- */ |
9a8b0c8d |
249 | |
b8e79eeb |
250 | if (rc == PGEN_TRY) |
251 | rc = smallenough(p->m); |
9a8b0c8d |
252 | |
253 | /* --- Done --- */ |
254 | |
255 | return (rc); |
256 | } |
257 | |
258 | /* --- @pfilt_muladd@ --- * |
259 | * |
260 | * Arguments: @pfilt *p@ = destination prime filtering context |
261 | * @const pfilt *q@ = source prime filtering context |
262 | * @mpw m@ = number to multiply by |
263 | * @mpw a@ = number to add |
264 | * |
265 | * Returns: One of the @PGEN@ result codes. |
266 | * |
267 | * Use: Multiplies the number in a prime filtering context by a |
268 | * small value and then adds a small value. The destination |
269 | * should either be uninitialized or the same as the source. |
270 | * |
271 | * Common things to do include multiplying by 2 and adding 0 to |
272 | * turn a prime into a jump for finding other primes with @q@ as |
273 | * a factor of @p - 1@, or multiplying by 2 and adding 1. |
274 | */ |
275 | |
276 | int pfilt_muladd(pfilt *p, const pfilt *q, mpw m, mpw a) |
277 | { |
278 | int rc = PGEN_TRY; |
279 | int i; |
280 | |
281 | /* --- Multiply the big number --- */ |
282 | |
283 | { |
87b63d41 |
284 | mp *d = mp_new(MP_LEN(q->m) + 2, q->m->f); |
9a8b0c8d |
285 | mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m); |
286 | mpx_uaddn(d->v, d->vl, a); |
9a8b0c8d |
287 | if (p == q) |
288 | mp_drop(p->m); |
289 | mp_shrink(d); |
290 | p->m = d; |
291 | } |
292 | |
293 | /* --- Gallivant through the residue table --- */ |
294 | |
295 | for (i = 0; i < NPRIME; i++) { |
296 | p->r[i] = (q->r[i] * m + a) % primetab[i]; |
297 | if (!p->r[i] && rc == PGEN_TRY) { |
298 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) |
299 | rc = PGEN_DONE; |
300 | else |
301 | rc = PGEN_FAIL; |
302 | } |
303 | } |
304 | |
b8e79eeb |
305 | /* --- Check for small primes --- */ |
9a8b0c8d |
306 | |
b8e79eeb |
307 | if (rc == PGEN_TRY) |
308 | rc = smallenough(p->m); |
9a8b0c8d |
309 | |
310 | /* --- Finished --- */ |
311 | |
312 | return (rc); |
313 | } |
314 | |
315 | /* --- @pfilt_jump@ --- * |
316 | * |
317 | * Arguments: @pfilt *p@ = pointer to prime filtering context |
318 | * @const pfilt *j@ = pointer to another filtering context |
319 | * |
320 | * Returns: One of the @PGEN@ result codes. |
321 | * |
322 | * Use: Steps a number by a large amount. Even so, jumping is much |
323 | * faster than initializing a new number. The test peformed is |
324 | * the same simple one used by @primetab_create@, so @PGEN_TRY@ |
325 | * results should be followed up by a Rabin-Miller test. |
326 | * |
327 | * Note that the number stored in the @j@ context is probably |
328 | * better off being even than prime. The important thing is |
329 | * that all of the residues for the number have already been |
330 | * computed. |
331 | */ |
332 | |
333 | int pfilt_jump(pfilt *p, const pfilt *j) |
334 | { |
335 | int rc = PGEN_TRY; |
336 | int i; |
337 | |
338 | /* --- Add the step on --- */ |
339 | |
340 | p->m = mp_add(p->m, p->m, j->m); |
341 | |
342 | /* --- Update the residue table --- */ |
343 | |
344 | for (i = 0; i < NPRIME; i++) { |
345 | p->r[i] = p->r[i] + j->r[i]; |
346 | if (p->r[i] > primetab[i]) |
347 | p->r[i] -= primetab[i]; |
348 | if (!p->r[i] && rc == PGEN_TRY) { |
349 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) |
350 | rc = PGEN_DONE; |
351 | else |
352 | rc = PGEN_FAIL; |
353 | } |
354 | } |
355 | |
b8e79eeb |
356 | /* --- Check for small primes --- */ |
9a8b0c8d |
357 | |
b8e79eeb |
358 | if (rc == PGEN_TRY) |
359 | rc = smallenough(p->m); |
9a8b0c8d |
360 | |
361 | /* --- Done --- */ |
362 | |
363 | return (rc); |
364 | } |
365 | |
366 | /*----- That's all, folks -------------------------------------------------*/ |