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

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