Commit | Line | Data |
---|---|---|
9a8b0c8d | 1 | /* -*-c-*- |
2 | * | |
b817bfc6 | 3 | * $Id: pfilt.c,v 1.6 2004/04/08 01:36:15 mdw Exp $ |
9a8b0c8d | 4 | * |
5 | * Finding and testing prime numbers | |
6 | * | |
7 | * (c) 1999 Straylight/Edgeware | |
8 | */ | |
9 | ||
45c0fd36 | 10 | /*----- Licensing notice --------------------------------------------------* |
9a8b0c8d | 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. | |
45c0fd36 | 18 | * |
9a8b0c8d | 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. | |
45c0fd36 | 23 | * |
9a8b0c8d | 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 | ||
9a8b0c8d | 30 | /*----- Header files ------------------------------------------------------*/ |
31 | ||
32 | #include "mp.h" | |
b8e79eeb | 33 | #include "mpint.h" |
9a8b0c8d | 34 | #include "pfilt.h" |
35 | #include "pgen.h" | |
36 | #include "primetab.h" | |
37 | ||
38 | /*----- Main code ---------------------------------------------------------*/ | |
39 | ||
b8e79eeb | 40 | /* --- @smallenough@ --- * |
41 | * | |
42 | * Arguments: @mp *m@ = integer to test | |
43 | * | |
44 | * Returns: One of the @PGEN@ result codes. | |
45 | * | |
46 | * Use: Assuming that @m@ has been tested by trial division on every | |
47 | * prime in the small-primes array, this function will return | |
48 | * @PGEN_DONE@ if the number is less than the square of the | |
49 | * largest small prime. | |
50 | */ | |
51 | ||
52 | static int smallenough(mp *m) | |
53 | { | |
54 | static mp *max = 0; | |
55 | int rc = PGEN_TRY; | |
56 | ||
57 | if (!max) { | |
58 | max = mp_fromuint(MP_NEW, MAXPRIME); | |
59 | max = mp_sqr(max, max); | |
60 | max->a->n--; /* Permanent allocation */ | |
61 | } | |
d31f4a79 MW |
62 | if (MP_CMP(m, <=, MP_ONE)) |
63 | rc = PGEN_FAIL; | |
64 | else if (MP_CMP(m, <, max)) | |
b8e79eeb | 65 | rc = PGEN_DONE; |
66 | return (rc); | |
67 | } | |
68 | ||
69 | /* --- @pfilt_smallfactor@ --- * | |
70 | * | |
71 | * Arguments: @mp *m@ = integer to test | |
72 | * | |
73 | * Returns: One of the @PGEN@ result codes. | |
74 | * | |
75 | * Use: Tests a number by dividing by a number of small primes. This | |
76 | * is a useful first step if you're testing random primes; for | |
77 | * sequential searches, @pfilt_create@ works better. | |
78 | */ | |
79 | ||
80 | int pfilt_smallfactor(mp *m) | |
81 | { | |
82 | int rc = PGEN_TRY; | |
83 | int i; | |
84 | size_t sz = MP_LEN(m); | |
34e4f738 | 85 | mparena *a = m->a ? m->a : MPARENA_GLOBAL; |
86 | mpw *v = mpalloc(a, sz); | |
b8e79eeb | 87 | |
88 | /* --- Fill in the residues --- */ | |
89 | ||
90 | for (i = 0; i < NPRIME; i++) { | |
91 | if (!mpx_udivn(v, v + sz, m->v, m->vl, primetab[i])) { | |
92 | if (MP_LEN(m) == 1 && m->v[0] == primetab[i]) | |
93 | rc = PGEN_DONE; | |
94 | else | |
95 | rc = PGEN_FAIL; | |
d31f4a79 | 96 | break; |
b8e79eeb | 97 | } |
98 | } | |
99 | ||
100 | /* --- Check for small primes --- */ | |
101 | ||
102 | if (rc == PGEN_TRY) | |
103 | rc = smallenough(m); | |
104 | ||
105 | /* --- Done --- */ | |
106 | ||
34e4f738 | 107 | mpfree(a, v); |
b8e79eeb | 108 | return (rc); |
109 | } | |
110 | ||
9a8b0c8d | 111 | /* --- @pfilt_create@ --- * |
112 | * | |
113 | * Arguments: @pfilt *p@ = pointer to prime filtering context | |
114 | * @mp *m@ = pointer to initial number to test | |
115 | * | |
116 | * Returns: One of the @PGEN@ result codes. | |
117 | * | |
118 | * Use: Tests an initial number for primality by computing its | |
119 | * residue modulo various small prime numbers. This is fairly | |
120 | * quick, but not particularly certain. If a @PGEN_TRY@ | |
121 | * result is returned, perform Rabin-Miller tests to confirm. | |
122 | */ | |
123 | ||
124 | int pfilt_create(pfilt *p, mp *m) | |
125 | { | |
126 | int rc = PGEN_TRY; | |
127 | int i; | |
b8e79eeb | 128 | size_t sz = MP_LEN(m); |
34e4f738 | 129 | mparena *a = m->a ? m->a : MPARENA_GLOBAL; |
130 | mpw *v = mpalloc(a, sz); | |
9a8b0c8d | 131 | |
132 | /* --- Take a copy of the number --- */ | |
133 | ||
134 | mp_shrink(m); | |
135 | p->m = MP_COPY(m); | |
136 | ||
137 | /* --- Fill in the residues --- */ | |
138 | ||
9a8b0c8d | 139 | for (i = 0; i < NPRIME; i++) { |
b8e79eeb | 140 | p->r[i] = mpx_udivn(v, v + sz, m->v, m->vl, primetab[i]); |
9a8b0c8d | 141 | if (!p->r[i] && rc == PGEN_TRY) { |
142 | if (MP_LEN(m) == 1 && m->v[0] == primetab[i]) | |
143 | rc = PGEN_DONE; | |
144 | else | |
145 | rc = PGEN_FAIL; | |
146 | } | |
147 | } | |
148 | ||
b8e79eeb | 149 | /* --- Check for small primes --- */ |
150 | ||
151 | if (rc == PGEN_TRY) | |
152 | rc = smallenough(m); | |
153 | ||
9a8b0c8d | 154 | /* --- Done --- */ |
155 | ||
34e4f738 | 156 | mpfree(a, v); |
9a8b0c8d | 157 | return (rc); |
158 | } | |
159 | ||
160 | /* --- @pfilt_destroy@ --- * | |
161 | * | |
162 | * Arguments: @pfilt *p@ = pointer to prime filtering context | |
163 | * | |
164 | * Returns: --- | |
165 | * | |
166 | * Use: Discards a context and all the resources it holds. | |
167 | */ | |
168 | ||
169 | void pfilt_destroy(pfilt *p) | |
170 | { | |
171 | mp_drop(p->m); | |
172 | } | |
173 | ||
174 | /* --- @pfilt_step@ --- * | |
175 | * | |
176 | * Arguments: @pfilt *p@ = pointer to prime filtering context | |
177 | * @mpw step@ = how much to step the number | |
178 | * | |
179 | * Returns: One of the @PGEN@ result codes. | |
180 | * | |
181 | * Use: Steps a number by a small amount. Stepping is much faster | |
182 | * than initializing with a new number. The test performed is | |
183 | * the same simple one used by @primetab_create@, so @PGEN_TRY@ | |
184 | * results should be followed up by a Rabin-Miller test. | |
185 | */ | |
186 | ||
187 | int pfilt_step(pfilt *p, mpw step) | |
188 | { | |
189 | int rc = PGEN_TRY; | |
190 | int i; | |
191 | ||
192 | /* --- Add the step on to the number --- */ | |
193 | ||
194 | p->m = mp_split(p->m); | |
195 | mp_ensure(p->m, MP_LEN(p->m) + 1); | |
196 | mpx_uaddn(p->m->v, p->m->vl, step); | |
197 | mp_shrink(p->m); | |
198 | ||
199 | /* --- Update the residue table --- */ | |
200 | ||
201 | for (i = 0; i < NPRIME; i++) { | |
202 | p->r[i] = (p->r[i] + step) % primetab[i]; | |
203 | if (!p->r[i] && rc == PGEN_TRY) { | |
204 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) | |
205 | rc = PGEN_DONE; | |
206 | else | |
207 | rc = PGEN_FAIL; | |
208 | } | |
209 | } | |
210 | ||
b8e79eeb | 211 | /* --- Check for small primes --- */ |
9a8b0c8d | 212 | |
b8e79eeb | 213 | if (rc == PGEN_TRY) |
214 | rc = smallenough(p->m); | |
9a8b0c8d | 215 | |
216 | /* --- Done --- */ | |
217 | ||
218 | return (rc); | |
219 | } | |
220 | ||
221 | /* --- @pfilt_muladd@ --- * | |
222 | * | |
223 | * Arguments: @pfilt *p@ = destination prime filtering context | |
224 | * @const pfilt *q@ = source prime filtering context | |
225 | * @mpw m@ = number to multiply by | |
226 | * @mpw a@ = number to add | |
227 | * | |
228 | * Returns: One of the @PGEN@ result codes. | |
229 | * | |
230 | * Use: Multiplies the number in a prime filtering context by a | |
231 | * small value and then adds a small value. The destination | |
232 | * should either be uninitialized or the same as the source. | |
233 | * | |
234 | * Common things to do include multiplying by 2 and adding 0 to | |
235 | * turn a prime into a jump for finding other primes with @q@ as | |
236 | * a factor of @p - 1@, or multiplying by 2 and adding 1. | |
237 | */ | |
238 | ||
239 | int pfilt_muladd(pfilt *p, const pfilt *q, mpw m, mpw a) | |
240 | { | |
241 | int rc = PGEN_TRY; | |
242 | int i; | |
243 | ||
244 | /* --- Multiply the big number --- */ | |
245 | ||
246 | { | |
87b63d41 | 247 | mp *d = mp_new(MP_LEN(q->m) + 2, q->m->f); |
9a8b0c8d | 248 | mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m); |
249 | mpx_uaddn(d->v, d->vl, a); | |
9a8b0c8d | 250 | if (p == q) |
251 | mp_drop(p->m); | |
252 | mp_shrink(d); | |
253 | p->m = d; | |
254 | } | |
255 | ||
256 | /* --- Gallivant through the residue table --- */ | |
45c0fd36 | 257 | |
9a8b0c8d | 258 | for (i = 0; i < NPRIME; i++) { |
259 | p->r[i] = (q->r[i] * m + a) % primetab[i]; | |
260 | if (!p->r[i] && rc == PGEN_TRY) { | |
261 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) | |
262 | rc = PGEN_DONE; | |
263 | else | |
264 | rc = PGEN_FAIL; | |
265 | } | |
266 | } | |
267 | ||
b8e79eeb | 268 | /* --- Check for small primes --- */ |
9a8b0c8d | 269 | |
b8e79eeb | 270 | if (rc == PGEN_TRY) |
271 | rc = smallenough(p->m); | |
9a8b0c8d | 272 | |
273 | /* --- Finished --- */ | |
274 | ||
275 | return (rc); | |
276 | } | |
277 | ||
278 | /* --- @pfilt_jump@ --- * | |
279 | * | |
280 | * Arguments: @pfilt *p@ = pointer to prime filtering context | |
281 | * @const pfilt *j@ = pointer to another filtering context | |
282 | * | |
283 | * Returns: One of the @PGEN@ result codes. | |
284 | * | |
285 | * Use: Steps a number by a large amount. Even so, jumping is much | |
286 | * faster than initializing a new number. The test peformed is | |
287 | * the same simple one used by @primetab_create@, so @PGEN_TRY@ | |
288 | * results should be followed up by a Rabin-Miller test. | |
289 | * | |
290 | * Note that the number stored in the @j@ context is probably | |
291 | * better off being even than prime. The important thing is | |
292 | * that all of the residues for the number have already been | |
293 | * computed. | |
294 | */ | |
295 | ||
296 | int pfilt_jump(pfilt *p, const pfilt *j) | |
297 | { | |
298 | int rc = PGEN_TRY; | |
299 | int i; | |
300 | ||
301 | /* --- Add the step on --- */ | |
302 | ||
303 | p->m = mp_add(p->m, p->m, j->m); | |
304 | ||
305 | /* --- Update the residue table --- */ | |
306 | ||
307 | for (i = 0; i < NPRIME; i++) { | |
308 | p->r[i] = p->r[i] + j->r[i]; | |
309 | if (p->r[i] > primetab[i]) | |
310 | p->r[i] -= primetab[i]; | |
311 | if (!p->r[i] && rc == PGEN_TRY) { | |
312 | if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i]) | |
313 | rc = PGEN_DONE; | |
314 | else | |
315 | rc = PGEN_FAIL; | |
316 | } | |
317 | } | |
318 | ||
b8e79eeb | 319 | /* --- Check for small primes --- */ |
9a8b0c8d | 320 | |
b8e79eeb | 321 | if (rc == PGEN_TRY) |
322 | rc = smallenough(p->m); | |
9a8b0c8d | 323 | |
324 | /* --- Done --- */ | |
325 | ||
326 | return (rc); | |
327 | } | |
328 | ||
329 | /*----- That's all, folks -------------------------------------------------*/ |