chiark / gitweb /
math/mpx-mul4-amd64-sse2.S: SSE2 multipliers for AMD64.
[catacomb] / 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 -------------------------------------------------*/