git.distorted.org.uk Git - u/mdw/catacomb/blob - pgen.c

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