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

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