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

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