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