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

   1 /* -*-c-*-
   2  *
   3  * Simultaneous prime search
   4  *
   5  * (c) 2006 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 <stdlib.h>
  31
  32 #include <mLib/macros.h>
  33
  34 #include "mprand.h"
  35 #include "pgen.h"
  36
  37 /*----- Main code ---------------------------------------------------------*/
  38
  39 /* --- @rcmerge@ --- *
  40  *
  41  * Arguments:   @int a, b@ = a pair of @PGEN@ return codes
  42  *
  43  * Returns:     The overall return code capturing both.
  44  */
  45
  46 static int rcmerge(int a, int b)
  47 {
  48 #define FROB(state)                                                     \
  49     if (a == PGEN_##state || b == PGEN_##state) return (PGEN_##state)
  50   FROB(FAIL);
  51   FROB(TRY);
  52   FROB(PASS);
  53 #undef FROB
  54   return (PGEN_DONE);
  55 }
  56
  57 /* --- @pgen_simulstep@ --- *
  58  *
  59  * Step a collection of numbers simultaneously.
  60  */
  61
  62 int pgen_simulstep(int rq, pgen_event *ev, void *p)
  63 {
  64   pgen_simulctx *ss = p;
  65   pgen_simulprime *sp;
  66   int rc = PGEN_ABORT, nrc;
  67   unsigned i;
  68   mp *q;
  69
  70   assert(ss->n);
  71   switch (rq) {
  72     case PGEN_BEGIN:
  73       rc = PGEN_DONE;
  74       q = MP_NEW;
  75       for (i = 0; i < ss->n; i++) {
  76         sp = &ss->v[i];
  77         q = mp_mul(q, ss->step, sp->mul);
  78         if (MP_LEN(q) <= 1)
  79           sp->u.step = q->v[0];
  80         else {
  81           sp->u.jump = xmalloc(sizeof(pfilt));
  82           pfilt_create(sp->u.jump, q);
  83           sp->f |= PGENF_JUMP;
  84         }
  85         q = mp_mul(q, ev->m, sp->mul);
  86         q = mp_add(q, q, sp->add);
  87         nrc = pfilt_create(&sp->p, q);
  88         rc = rcmerge(rc, nrc);
  89       }
  90       MP_DROP(q);
  91       if (rc != PGEN_FAIL)
  92         goto done;
  93     case PGEN_TRY:
  94       for (;;) {
  95         rc = PGEN_DONE;
  96         for (i = 0; i < ss->n; i++) {
  97           sp = &ss->v[i];
  98           if (sp->f & PGENF_JUMP)
  99             nrc = pfilt_jump(&sp->p, sp->u.jump);
 100           else
 101             nrc = pfilt_step(&sp->p, sp->u.step);
 102           rc = rcmerge(rc, nrc);
 103         }
 104         if (rc != PGEN_FAIL)
 105           goto done;
 106       }
 107     done:
 108       mp_drop(ev->m);
 109       ev->m = MP_COPY(ss->v[0].p.m);
 110       break;
 111     case PGEN_DONE:
 112       for (i = 0; i < ss->n; i++) {
 113         sp = &ss->v[i];
 114         if (sp->f & PGENF_JUMP) {
 115           pfilt_destroy(sp->u.jump);
 116           xfree(sp->u.jump);
 117         }
 118         if (sp->f & PGENF_KEEP)
 119           sp->u.x = MP_COPY(sp->p.m);
 120         pfilt_destroy(&sp->p);
 121       }
 122       rc = PGEN_DONE;
 123       break;
 124   }
 125   return (rc);
 126 }
 127
 128 /* --- @pgen_simultest@ --- *
 129  *
 130  * Test a collection of numbers simultaneously.
 131  */
 132
 133 int pgen_simultest(int rq, pgen_event *ev, void *p)
 134 {
 135   pgen_simulctx *ss = p;
 136   pgen_simulprime *sp;
 137   int rc = -1;
 138   unsigned i;
 139   mp *m;
 140
 141   assert(ss->n);
 142   switch (rq) {
 143     case PGEN_BEGIN:
 144       for (i = 0; i < ss->n; i++)
 145         rabin_create(&ss->v[i].r, ss->v[i].p.m);
 146       rc = PGEN_TRY;
 147       break;
 148     case PGEN_TRY:
 149       m = MP_NEW;
 150       for (i = 0; i < ss->n; i++) {
 151         sp = &ss->v[i];
 152         m = mprand_range(m, sp->p.m, ev->r, 0);
 153         rc = rabin_test(&sp->r, m);
 154         if (rc != PGEN_PASS)
 155           break;
 156       }
 157       mp_drop(m);
 158       break;
 159     case PGEN_DONE:
 160       for (i = 0; i < ss->n; i++)
 161         rabin_destroy(&ss->v[i].r);
 162       rc = PGEN_DONE;
 163       break;
 164   }
 165   return (rc);
 166 }
 167
 168 /*----- That's all, folks -------------------------------------------------*/