chiark - git - mdw - catacomb/blob - pgen.c

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