git.distorted.org.uk Git - u/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_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   int p;
 178
 179   enum { P_STEP, P_TEST };
 180
 181   /* --- Set up the initial event block --- */
 182
 183   ev.name = name;
 184   if (m)
 185     ev.m = MP_COPY(m);
 186   else
 187     ev.m = 0;
 188   ev.steps = 0;
 189   ev.tests = 0;
 190   ev.r = fibrand_create(0);
 191
 192   /* --- Tell the event handler we're under way --- */
 193
 194   if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
 195     ev.r->ops->destroy(ev.r);
 196     return (0);
 197   }
 198
 199   /* --- Set up for the initial call --- */
 200
 201   proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
 202
 203   /* --- Enter the great maelstrom of state transitions --- */
 204
 205   for (;;) {
 206     unsigned act = 0;
 207
 208 #define A_STEP 1u
 209 #define A_TEST 2u
 210 #define A_EVENT 4u
 211 #define A_ENDTEST 8u
 212 #define A_ENDSTEP 16u
 213 #define A_DONE 32u
 214
 215     /* --- Call the procedure and decide what to do next --- */
 216
 217     rc = proc(rq, &ev, ctx);
 218     switch (rc) {
 219       case PGEN_TRY:
 220         if (p == P_TEST)
 221           rq = PGEN_TRY;
 222         else {
 223           act |= A_EVENT;
 224           proc = test; ctx = tctx; p = P_TEST;
 225           rq = PGEN_BEGIN;
 226         }
 227         break;
 228       case PGEN_PASS:
 229         act |= A_TEST | A_EVENT;
 230         if (p == P_TEST)
 231           rq = PGEN_TRY;
 232         else {
 233           proc = test; ctx = tctx; p = P_TEST;
 234           rq = PGEN_BEGIN;
 235         }
 236         break;
 237       case PGEN_FAIL:
 238         act |= A_STEP;
 239         if (p == P_TEST) {
 240           act |= A_ENDTEST | A_EVENT;
 241           proc = step; ctx = sctx; p = P_STEP;
 242         }
 243         rq = PGEN_TRY;
 244         break;
 245       case PGEN_DONE:
 246         act |= A_EVENT | A_DONE | A_ENDSTEP;
 247         if (p == P_TEST)
 248           act |= A_ENDTEST;
 249         break;
 250       case PGEN_ABORT:
 251         act |= A_EVENT | A_DONE;
 252         if (p == P_TEST || rq == PGEN_TRY)
 253           act |= A_ENDSTEP;
 254         if (p == P_TEST && rq != PGEN_BEGIN)
 255           act |= A_ENDTEST;
 256         break;
 257       default:
 258         assert(((void)"Invalid response from function", 0));
 259         break;
 260     }
 261
 262     /* --- If decrementing counters is requested, do that --- */
 263
 264     if ((act & A_STEP) && steps) {
 265       ev.steps++;
 266       if (ev.steps == steps) {
 267         act |= A_EVENT | A_ENDSTEP | A_DONE;
 268         rc = PGEN_ABORT;
 269       }
 270       ev.tests = 0;
 271     }
 272
 273     if ((act & A_TEST) && tests) {
 274       ev.tests++;
 275       if (ev.tests == tests) {
 276         act |= A_ENDTEST | A_ENDSTEP | A_DONE;
 277         rc = PGEN_DONE;
 278       }
 279     }
 280
 281     /* --- Report an event if so directed --- */
 282
 283     if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
 284       rc = PGEN_ABORT;
 285       if (!(act & A_DONE)) {
 286         act |= A_ENDSTEP | A_DONE;
 287         if (p == P_TEST)
 288           act |= A_ENDTEST;
 289       }
 290     }
 291
 292     /* --- Close down tester and stepper functions --- */
 293
 294     if (act & A_ENDTEST)
 295       test(PGEN_DONE, &ev, tctx);
 296     if (act & A_ENDSTEP)
 297       step(PGEN_DONE, &ev, sctx);
 298
 299     /* --- Stop the entire test if necessary --- */
 300
 301     if (act & A_DONE)
 302       break;
 303   }
 304
 305   /* --- Tidy up and return --- */
 306
 307   if (rc == PGEN_ABORT) {
 308     mp_drop(ev.m);
 309     ev.m = 0;
 310   }
 311   ev.r->ops->destroy(ev.r);
 312   mp_drop(d);
 313
 314   return (ev.m);
 315 }
 316
 317 /* --- @pgen_primep@ --- *
 318  *
 319  * Arguments:   @mp *p@ = a number to check
 320  *              @grand *gr@ = a random number source
 321  *
 322  * Returns:     Nonzero if @p@ is really prime.
 323  */
 324
 325 int pgen_primep(mp *p, grand *gr)
 326 {
 327   int i;
 328   rabin r;
 329   mp *x = MP_NEW;
 330
 331   if (MP_NEGP(p)) return (0);
 332   switch (pfilt_smallfactor(p)) {
 333     case PGEN_DONE: return (1);
 334     case PGEN_FAIL: return (0);
 335   }
 336   rabin_create(&r, p);
 337   for (i = 32; i; i--) {
 338     x = mprand_range(x, p, gr, 0);
 339     if (rabin_rtest(&r, x) == PGEN_FAIL)
 340       break;
 341   }
 342   MP_DROP(x);
 343   rabin_destroy(&r);
 344   return (!i);
 345 }
 346
 347 /*----- Test rig ----------------------------------------------------------*/
 348
 349 #ifdef TEST_RIG
 350
 351 #include <mLib/testrig.h>
 352
 353 static int verify(dstr *v)
 354 {
 355   mp *m = *(mp **)v[0].buf;
 356   mp *q = *(mp **)v[1].buf;
 357   mp *p;
 358   int ok = 1;
 359
 360   pgen_filterctx pf;
 361   rabin r;
 362
 363   pf.step = 2;
 364   p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
 365            rabin_iters(mp_bits(m)), pgen_test, &r);
 366   if (!p || !MP_EQ(p, q)) {
 367     fputs("\n*** pgen failed", stderr);
 368     fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
 369     fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
 370     fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
 371     fputc('\n', stderr);
 372     ok = 0;
 373   }
 374
 375   mp_drop(m);
 376   mp_drop(q);
 377   mp_drop(p);
 378   assert(mparena_count(MPARENA_GLOBAL) == 0);
 379   return (ok);
 380 }
 381
 382 static test_chunk tests[] = {
 383   { "pgen", verify, { &type_mp, &type_mp, 0 } },
 384   { 0, 0, { 0 } }
 385 };
 386
 387 int main(int argc, char *argv[])
 388 {
 389   sub_init();
 390   test_run(argc, argv, tests, SRCDIR "/tests/pgen");
 391   return (0);
 392 }
 393 #endif
 394
 395 /*----- That's all, folks -------------------------------------------------*/