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

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