chiark / gitweb /
math/group-parse.c (group-parse): Parse binary-group descriptions.
[catacomb] / math / pgen.c
CommitLineData
0f5ec153 1/* -*-c-*-
0f5ec153 2 *
581c854e 3 * Prime generation glue
0f5ec153 4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
0f5ec153 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 *
0f5ec153 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 *
0f5ec153 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
0f5ec153 28/*----- Header files ------------------------------------------------------*/
29
581c854e 30#include <assert.h>
31#include <stdio.h>
32#include <stdlib.h>
33#include <string.h>
34
35#include "fibrand.h"
36#include "grand.h"
0f5ec153 37#include "mp.h"
581c854e 38#include "mprand.h"
0f5ec153 39#include "pgen.h"
581c854e 40#include "pfilt.h"
0f5ec153 41#include "rabin.h"
42
581c854e 43/*----- Standard prime filter ---------------------------------------------*/
0f5ec153 44
581c854e 45/* --- @pgen_filter@ --- */
0f5ec153 46
581c854e 47int pgen_filter(int rq, pgen_event *ev, void *p)
0f5ec153 48{
581c854e 49 pgen_filterctx *f = p;
57541e97 50 int rc = PGEN_FAIL;
581c854e 51
52 switch (rq) {
53 case PGEN_BEGIN:
54 rc = pfilt_create(&f->f, ev->m);
55 mp_drop(ev->m);
56 break;
57 case PGEN_TRY:
58 mp_drop(ev->m);
581c854e 59 break;
60 case PGEN_DONE:
61 pfilt_destroy(&f->f);
62 return (PGEN_DONE);
57541e97 63 default:
64 rc = PGEN_ABORT;
65 break;
0f5ec153 66 }
45c0fd36 67
57541e97 68 if (rc == PGEN_FAIL && !((f->step | f->f.m->v[0]) & 1))
69 rc = pfilt_step(&f->f, 1);
581c854e 70 while (rc == PGEN_FAIL)
71 rc = pfilt_step(&f->f, f->step);
72 ev->m = MP_COPY(f->f.m);
0f5ec153 73 return (rc);
74}
75
581c854e 76/* --- @pgen_jump@ --- *
0f5ec153 77 *
581c854e 78 * Similar to the standard @pgen_filter@, but jumps in large steps rather
79 * than small ones.
0f5ec153 80 */
81
581c854e 82int pgen_jump(int rq, pgen_event *ev, void *p)
0f5ec153 83{
581c854e 84 pgen_jumpctx *f = p;
85 int rc = PGEN_ABORT;
0f5ec153 86
581c854e 87 switch (rq) {
8b021c3f 88 case PGEN_BEGIN: {
89 mp *g = MP_NEW;
90 mp_gcd(&g, 0, 0, ev->m, f->j->m);
91 if (MP_CMP(g, >, MP_ONE)) {
92 mp_drop(g);
93 return (PGEN_ABORT);
94 }
95 mp_drop(g);
581c854e 96 rc = pfilt_create(&f->f, ev->m);
97 mp_drop(ev->m);
8b021c3f 98 } break;
581c854e 99 case PGEN_TRY:
100 mp_drop(ev->m);
101 rc = pfilt_jump(&f->f, f->j);
102 break;
103 case PGEN_DONE:
104 pfilt_destroy(&f->f);
105 return (PGEN_DONE);
106 }
45c0fd36 107
581c854e 108 while (rc == PGEN_FAIL)
109 rc = pfilt_jump(&f->f, f->j);
110 ev->m = MP_COPY(f->f.m);
111 return (rc);
112}
0f5ec153 113
581c854e 114/*----- Standard prime test -----------------------------------------------*/
0f5ec153 115
581c854e 116/* --- @pgen_test@ --- */
0f5ec153 117
581c854e 118int pgen_test(int rq, pgen_event *ev, void *p)
119{
120 rabin *r = p;
121 int rc = PGEN_ABORT;
0f5ec153 122
581c854e 123 switch (rq) {
124 case PGEN_BEGIN:
125 rabin_create(r, ev->m);
126 rc = PGEN_TRY;
127 break;
283b9af0 128 case PGEN_TRY:
129 if (!ev->tests)
130 rc = rabin_rtest(r, MP_TWO);
131 else {
132 mp *a = mprand_range(MP_NEW, ev->m, ev->r, 0);
133 rc = rabin_rtest(r, a);
134 mp_drop(a);
135 }
136 break;
581c854e 137 case PGEN_DONE:
138 rabin_destroy(r);
139 rc = PGEN_DONE;
140 break;
0f5ec153 141 }
142
0f5ec153 143 return (rc);
144}
145
581c854e 146/*----- The main driver ---------------------------------------------------*/
147
148/* --- @pgen@ --- *
bd98b2df 149 *
581c854e 150 * Arguments: @const char *name@ = name of the value being searched for
151 * @mp *d@ = destination for the result integer
152 * @mp *m@ = start value to pass to stepper
153 * @pgen_proc *event@ = event handler function
154 * @void *ectx@ = context argument for event andler
155 * @unsigned steps@ = number of steps to take in search
156 * @pgen_proc *step@ = stepper function to use
157 * @void *sctx@ = context argument for stepper
158 * @unsigned tests@ = number of tests to make
159 * @pgen_proc *test@ = tester function to use
160 * @void *tctx@ = context argument for tester
bd98b2df 161 *
581c854e 162 * Returns: Pointer to final result, or null.
bd98b2df 163 *
581c854e 164 * Use: A generalized prime-number search skeleton. Yes, that's a
165 * scary number of arguments.
bd98b2df 166 */
167
581c854e 168mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
169 unsigned steps, pgen_proc *step, void *sctx,
170 unsigned tests, pgen_proc *test, void *tctx)
bd98b2df 171{
581c854e 172 pgen_event ev;
173 int rq, rc;
174 pgen_proc *proc;
175 void *ctx;
4108c8d2 176 int p;
177
178 enum { P_STEP, P_TEST };
bd98b2df 179
581c854e 180 /* --- Set up the initial event block --- */
bd98b2df 181
581c854e 182 ev.name = name;
183 if (m)
184 ev.m = MP_COPY(m);
185 else
186 ev.m = 0;
283b9af0 187 ev.steps = 0;
188 ev.tests = 0;
581c854e 189 ev.r = fibrand_create(0);
bd98b2df 190
581c854e 191 /* --- Tell the event handler we're under way --- */
bd98b2df 192
4108c8d2 193 if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
194 ev.r->ops->destroy(ev.r);
581c854e 195 return (0);
4108c8d2 196 }
bd98b2df 197
581c854e 198 /* --- Set up for the initial call --- */
bd98b2df 199
4108c8d2 200 proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
0f5ec153 201
581c854e 202 /* --- Enter the great maelstrom of state transitions --- */
0f5ec153 203
581c854e 204 for (;;) {
205 unsigned act = 0;
206
4108c8d2 207#define A_STEP 1u
208#define A_TEST 2u
209#define A_EVENT 4u
210#define A_ENDTEST 8u
211#define A_ENDSTEP 16u
212#define A_DONE 32u
581c854e 213
214 /* --- Call the procedure and decide what to do next --- */
215
216 rc = proc(rq, &ev, ctx);
217 switch (rc) {
218 case PGEN_TRY:
4108c8d2 219 if (p == P_TEST)
581c854e 220 rq = PGEN_TRY;
221 else {
222 act |= A_EVENT;
4108c8d2 223 proc = test; ctx = tctx; p = P_TEST;
581c854e 224 rq = PGEN_BEGIN;
225 }
226 break;
227 case PGEN_PASS:
228 act |= A_TEST | A_EVENT;
4108c8d2 229 if (p == P_TEST)
581c854e 230 rq = PGEN_TRY;
231 else {
4108c8d2 232 proc = test; ctx = tctx; p = P_TEST;
581c854e 233 rq = PGEN_BEGIN;
234 }
235 break;
236 case PGEN_FAIL:
237 act |= A_STEP;
4108c8d2 238 if (p == P_TEST) {
581c854e 239 act |= A_ENDTEST | A_EVENT;
4108c8d2 240 proc = step; ctx = sctx; p = P_STEP;
581c854e 241 }
242 rq = PGEN_TRY;
243 break;
244 case PGEN_DONE:
245 act |= A_EVENT | A_DONE | A_ENDSTEP;
4108c8d2 246 if (p == P_TEST)
581c854e 247 act |= A_ENDTEST;
248 break;
249 case PGEN_ABORT:
250 act |= A_EVENT | A_DONE;
4108c8d2 251 if (p == P_TEST || rq == PGEN_TRY)
581c854e 252 act |= A_ENDSTEP;
4108c8d2 253 if (p == P_TEST && rq != PGEN_BEGIN)
581c854e 254 act |= A_ENDTEST;
255 break;
256 default:
257 assert(((void)"Invalid response from function", 0));
258 break;
259 }
0f5ec153 260
581c854e 261 /* --- If decrementing counters is requested, do that --- */
0f5ec153 262
581c854e 263 if ((act & A_STEP) && steps) {
283b9af0 264 ev.steps++;
265 if (ev.steps == steps) {
581c854e 266 act |= A_EVENT | A_ENDSTEP | A_DONE;
267 rc = PGEN_ABORT;
268 }
283b9af0 269 ev.tests = 0;
581c854e 270 }
0f5ec153 271
581c854e 272 if ((act & A_TEST) && tests) {
283b9af0 273 ev.tests++;
274 if (ev.tests == tests) {
581c854e 275 act |= A_ENDTEST | A_ENDSTEP | A_DONE;
276 rc = PGEN_DONE;
277 }
0f5ec153 278 }
0f5ec153 279
581c854e 280 /* --- Report an event if so directed --- */
bd98b2df 281
581c854e 282 if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
283 rc = PGEN_ABORT;
284 if (!(act & A_DONE)) {
285 act |= A_ENDSTEP | A_DONE;
8501f5f0 286 if (p == P_TEST && rq != PGEN_BEGIN)
581c854e 287 act |= A_ENDTEST;
288 }
45c0fd36 289 }
bd98b2df 290
581c854e 291 /* --- Close down tester and stepper functions --- */
0f5ec153 292
581c854e 293 if (act & A_ENDTEST)
294 test(PGEN_DONE, &ev, tctx);
295 if (act & A_ENDSTEP)
296 step(PGEN_DONE, &ev, sctx);
297
298 /* --- Stop the entire test if necessary --- */
299
300 if (act & A_DONE)
301 break;
302 }
303
304 /* --- Tidy up and return --- */
305
306 if (rc == PGEN_ABORT) {
307 mp_drop(ev.m);
308 ev.m = 0;
309 }
310 ev.r->ops->destroy(ev.r);
eab06f16 311 mp_drop(d);
581c854e 312
313 return (ev.m);
0f5ec153 314}
315
34e4f738 316/* --- @pgen_primep@ --- *
317 *
318 * Arguments: @mp *p@ = a number to check
319 * @grand *gr@ = a random number source
320 *
321 * Returns: Nonzero if @p@ is really prime.
322 */
323
324int pgen_primep(mp *p, grand *gr)
325{
4108c8d2 326 int i;
34e4f738 327 rabin r;
328 mp *x = MP_NEW;
329
a69a3efd 330 if (MP_NEGP(p)) return (0);
34e4f738 331 switch (pfilt_smallfactor(p)) {
332 case PGEN_DONE: return (1);
333 case PGEN_FAIL: return (0);
334 }
335 rabin_create(&r, p);
4108c8d2 336 for (i = 32; i; i--) {
34e4f738 337 x = mprand_range(x, p, gr, 0);
338 if (rabin_rtest(&r, x) == PGEN_FAIL)
339 break;
34e4f738 340 }
341 MP_DROP(x);
342 rabin_destroy(&r);
343 return (!i);
344}
345
581c854e 346/*----- Test rig ----------------------------------------------------------*/
0f5ec153 347
348#ifdef TEST_RIG
349
350#include <mLib/testrig.h>
351
d31f4a79
MW
352static int t_primep(dstr *v)
353{
354 mp *m = *(mp **)v[0].buf;
355 int e = *(int *)v[1].buf;
356 int r;
357 grand *rng;
358 int ok = 1;
359
360 rng = fibrand_create(0);
361 r = pgen_primep(m, rng);
362 GR_DESTROY(rng);
363 if (e != r) {
364 fputs("\n*** primep failed", stderr);
365 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
366 fprintf(stderr, "\nexpected %d", e);
367 fprintf(stderr, "\nreported %d", r);
368 fputc('\n', stderr);
369 ok = 0;
370 }
371
372 mp_drop(m);
373 assert(mparena_count(MPARENA_GLOBAL) == 0);
45c0fd36 374 return (ok);
d31f4a79
MW
375}
376
0f5ec153 377static int verify(dstr *v)
378{
379 mp *m = *(mp **)v[0].buf;
581c854e 380 mp *q = *(mp **)v[1].buf;
381 mp *p;
0f5ec153 382 int ok = 1;
383
581c854e 384 pgen_filterctx pf;
385 rabin r;
0f5ec153 386
581c854e 387 pf.step = 2;
388 p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
389 rabin_iters(mp_bits(m)), pgen_test, &r);
22bab86c 390 if (!p || !MP_EQ(p, q)) {
581c854e 391 fputs("\n*** pgen failed", stderr);
0f5ec153 392 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
581c854e 393 fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
394 fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
0f5ec153 395 fputc('\n', stderr);
396 ok = 0;
397 }
398
399 mp_drop(m);
581c854e 400 mp_drop(q);
eab06f16 401 mp_drop(p);
d94f85ac 402 assert(mparena_count(MPARENA_GLOBAL) == 0);
0f5ec153 403 return (ok);
404}
405
406static test_chunk tests[] = {
407 { "pgen", verify, { &type_mp, &type_mp, 0 } },
d31f4a79 408 { "primep", t_primep, { &type_mp, &type_int, 0 } },
0f5ec153 409 { 0, 0, { 0 } }
410};
411
412int main(int argc, char *argv[])
413{
414 sub_init();
0f00dc4c 415 test_run(argc, argv, tests, SRCDIR "/t/pgen");
0f5ec153 416 return (0);
417}
0f5ec153 418#endif
419
420/*----- That's all, folks -------------------------------------------------*/