From 0f5ec153dcb1c20a765d0420cc189ebac1b9682d Mon Sep 17 00:00:00 2001 From: mdw Date: Fri, 19 Nov 1999 13:17:57 +0000 Subject: [PATCH] Prime number generator and tester. --- pgen.c | 268 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pgen.h | 137 +++++++++++++++++++++++++++++++++ rabin.c | 156 +++++++++++++++++++++++++++++++++++++ rabin.h | 115 ++++++++++++++++++++++++++++ 4 files changed, 676 insertions(+) create mode 100644 pgen.c create mode 100644 pgen.h create mode 100644 rabin.c create mode 100644 rabin.h diff --git a/pgen.c b/pgen.c new file mode 100644 index 0000000..503b153 --- /dev/null +++ b/pgen.c @@ -0,0 +1,268 @@ +/* -*-c-*- + * + * $Id: pgen.c,v 1.1 1999/11/19 13:17:57 mdw Exp $ + * + * Finding and testing prime numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: pgen.c,v $ + * Revision 1.1 1999/11/19 13:17:57 mdw + * Prime number generator and tester. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpmont.h" +#include "pgen.h" +#include "ptab.h" +#include "rabin.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @pgen_create@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @mp *m@ = pointer to initial number to test + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Tests an initial number for primality by computing its + * residue modulo various small prime numbers. This is fairly + * quick, but not particularly certain. If a @PGEN_MAYBE@ + * result is returned, perform Rabin-Miller tests to confirm. + */ + +int pgen_create(pgen *p, mp *m) +{ + int rc = PGEN_MAYBE; + int i; + mp *r = MP_NEW; + mpw qw; + mp q; + + /* --- Take a copy of the number --- */ + + mp_shrink(m); + p->m = MP_COPY(m); + + /* --- Fill in the residues --- */ + + mp_build(&q, &qw, &qw + 1); + for (i = 0; i < NPRIME; i++) { + qw = ptab[i]; + mp_div(0, &r, m, &q); + p->r[i] = r->v[0]; + if (!p->r[i] && rc == PGEN_MAYBE) { + if (MP_LEN(m) == 1 && m->v[0] == ptab[i]) + rc = PGEN_PRIME; + else + rc = PGEN_COMPOSITE; + } + } + + /* --- Done --- */ + + mp_drop(r); + return (rc); +} + +/* --- @pgen_destroy@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * + * Returns: --- + * + * Use: Discards a context and all the resources it holds. + */ + +void pgen_destroy(pgen *p) +{ + mp_drop(p->m); +} + +/* --- @pgen_step@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @mpw step@ = how much to step the number + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Steps a number by a small amount. Stepping is much faster + * than initializing with a new number. The test performed is + * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ + * results should be followed up by a Rabin-Miller test. + */ + +int pgen_step(pgen *p, mpw step) +{ + mp s; + int rc = PGEN_MAYBE; + int i; + + /* --- Add the step on to the number --- */ + + mp_build(&s, &step, &step + 1); + p->m = mp_add(p->m, p->m, &s); + + /* --- Update the residue table --- */ + + for (i = 0; i < NPRIME; i++) { + p->r[i] = (p->r[i] + step) % ptab[i]; + if (!p->r[i] && rc == PGEN_MAYBE) { + if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i]) + rc = PGEN_PRIME; + else + rc = PGEN_COMPOSITE; + } + } + + /* --- Done --- */ + + return (rc); +} + +/* --- @pgen_jump@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @pgen *j@ = pointer to another generation context + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Steps a number by a large amount. Even so, jumping is much + * faster than initializing a new number. The test peformed is + * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ + * results should be followed up by a Rabin-Miller test. + * + * Note that the number stored in the @j@ context is probably + * better off being even than prime. The important thing is + * that all of the residues for the number have already been + * computed. + */ + +int pgen_jump(pgen *p, pgen *j) +{ + int rc = PGEN_MAYBE; + int i; + + /* --- Add the step on --- */ + + p->m = mp_add(p->m, p->m, j->m); + + /* --- Update the residue table --- */ + + for (i = 0; i < NPRIME; i++) { + p->r[i] = p->r[i] + j->r[i]; + if (p->r[i] > ptab[i]) + p->r[i] -= ptab[i]; + if (!p->r[i] && rc == PGEN_MAYBE) { + if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i]) + rc = PGEN_PRIME; + else + rc = PGEN_COMPOSITE; + } + } + + /* --- Done --- */ + + return (rc); +} + +/*----- Test code ---------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int verify(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mp *r = *(mp **)v[1].buf; + pgen p; + mpw gw; + mp g; + int rc; + static char baton[] = "/-\\|"; + char *b = baton; + int ok = 1; + + mp_build(&g, &gw, &gw + 1); + rc = pgen_create(&p, m); + for (;;) { + if (rc == PGEN_PRIME) + break; + if (rc == PGEN_MAYBE) { + int i; + rabin r; + rabin_create(&r, p.m); + for (i = 0; i < 5; i++) { + gw = rand(); + if (rabin_test(&r, &g) == PGEN_COMPOSITE) + break; + } + rabin_destroy(&r); + if (i == 5) + break; + putc(*b++, stderr); + putc('\b', stderr); + if (!*b) + b = baton; + } + rc = pgen_step(&p, 2); + } + + if (MP_CMP(p.m, !=, r)) { + fputs("\n*** failed.", stderr); + fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); + fputs("\np = ", stderr); mp_writefile(p.m, stderr, 10); + fputs("\nr = ", stderr); mp_writefile(r, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + mp_drop(m); + mp_drop(r); + pgen_destroy(&p); + return (ok); +} + +static test_chunk tests[] = { + { "pgen", verify, { &type_mp, &type_mp, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/pgen"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/pgen.h b/pgen.h new file mode 100644 index 0000000..536ba0a --- /dev/null +++ b/pgen.h @@ -0,0 +1,137 @@ +/* -*-c-*- + * + * $Id: pgen.h,v 1.1 1999/11/19 13:17:57 mdw Exp $ + * + * Finding and testing prime numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: pgen.h,v $ + * Revision 1.1 1999/11/19 13:17:57 mdw + * Prime number generator and tester. + * + */ + +#ifndef PGEN_H +#define PGEN_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MP_H +# include "mp.h" +#endif + +#ifndef PTAB_H +# include "ptab.h" +#endif + +/*----- Constants ---------------------------------------------------------*/ + +#define PGEN_COMPOSITE (-1) /* Number is definitely composite */ +#define PGEN_MAYBE (0) /* Number may be prime */ +#define PGEN_PRIME (1) /* Number is definitely prime */ + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct pgen { + mp *m; + unsigned char r[NPRIME]; +} pgen; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @pgen_create@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @mp *m@ = pointer to initial number to test + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Tests an initial number for primality by computing its + * residue modulo various small prime numbers. This is fairly + * quick, but not particularly certain. If a @PGEN_MAYBE@ + * result is returned, perform Rabin-Miller tests to confirm. + */ + +extern int pgen_create(pgen */*p*/, mp */*m*/); + +/* --- @pgen_destroy@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * + * Returns: --- + * + * Use: Discards a context and all the resources it holds. + */ + +extern void pgen_destroy(pgen */*p*/); + +/* --- @pgen_step@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @mpw step@ = how much to step the number + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Steps a number by a small amount. Stepping is much faster + * than initializing with a new number. The test performed is + * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ + * results should be followed up by a Rabin-Miller test. + */ + +extern int pgen_step(pgen */*p*/, mpw /*step*/); + +/* --- @pgen_jump@ --- * + * + * Arguments: @pgen *p@ = pointer to prime generation context + * @pgen *j@ = pointer to another generation context + * + * Returns: One of the @PGEN@ constants above. + * + * Use: Steps a number by a large amount. Even so, jumping is much + * faster than initializing a new number. The test peformed is + * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ + * results should be followed up by a Rabin-Miller test. + * + * Note that the number stored in the @j@ context is probably + * better off being even than prime. The important thing is + * that all of the residues for the number have already been + * computed. + */ + +extern int pgen_jump(pgen */*p*/, pgen */*j*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/rabin.c b/rabin.c new file mode 100644 index 0000000..095c21d --- /dev/null +++ b/rabin.c @@ -0,0 +1,156 @@ +/* -*-c-*- + * + * $Id: rabin.c,v 1.1 1999/11/19 13:17:57 mdw Exp $ + * + * Miller-Rabin primality test + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: rabin.c,v $ + * Revision 1.1 1999/11/19 13:17:57 mdw + * Prime number generator and tester. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpmont.h" +#include "pgen.h" +#include "rabin.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @rabin_create@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * @mp *m@ = pointer to number to test + * + * Returns: --- + * + * Use: Precomputes some useful values for performing the + * Miller-Rabin probabilistic primality test. + */ + +void rabin_create(rabin *r, mp *m) +{ + mp *m1 = mp_sub(MP_NEW, m, MP_ONE); + mpscan sc; + size_t s; + + /* --- Find @r@ and @s@ --- */ + + mpmont_create(&r->mm, m); + mp_scan(&sc, m1); + s = 0; + while (mp_step(&sc)) { + if (mp_bit(&sc)) + break; + s++; + } + r->s = s; + r->r = mp_lsr(MP_NEW, m1, s); + + /* --- Compute %$(m - 1)R \bmod m$% --- */ + + r->m1 = mpmont_mul(&r->mm, MP_NEW, m1, r->mm.r2); + mp_drop(m1); +} + +/* --- @rabin_destroy@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * + * Returns: --- + * + * Use: Disposes of a Rabin-Miller context when it's no longer + * needed. + */ + +void rabin_destroy(rabin *r) +{ + mp_drop(r->r); + mp_drop(r->m1); + mpmont_destroy(&r->mm); +} + +/* --- @rabin_test@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * @mp *g@ = base to test the number against + * + * Returns: Either @PGEN_COMPOSITE@ if the test failed, or @PGEN_MAYBE@ + * if it succeeded. + * + * Use: Performs a single iteration of the Rabin-Miller primality + * test. + */ + +int rabin_test(rabin *r, mp *g) +{ + mp *y; + mp *dd, *spare = MP_NEW; + size_t j; + int rc = PGEN_COMPOSITE; + + /* --- Calculate %$y R = g^r R \bmod m$% --- * + * + * If %$y = 1$% or %$y = m - 1$% then %$m$% is prime. If course, note that + * @y@ here has an extra factor of %$R$%. + */ + + y = mpmont_expr(&r->mm, g, r->r); + if (MP_CMP(y, ==, r->mm.r) || MP_CMP(y, ==, r->m1)) { + rc = PGEN_MAYBE; + goto done; + } + + /* --- Now for the main loop --- * + * + * If %$y^{2^j} \ne m - 1$% for any %$0 \le j < s$% then %$m$% is + * composite. Of course, %$j = 0$% has already been tested. + */ + + for (j = 1; j < r->s; j++) { + dd = mpmont_mul(&r->mm, spare, y, y); spare = y; y = dd; + if (MP_CMP(y, ==, r->mm.r)) + break; + if (MP_CMP(y, ==, r->m1)) { + rc = PGEN_MAYBE; + break; + } + } + + /* --- Done --- */ + +done: + if (spare != MP_NEW) + MP_DROP(spare); + MP_DROP(y); + return (rc); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/rabin.h b/rabin.h new file mode 100644 index 0000000..ed08aeb --- /dev/null +++ b/rabin.h @@ -0,0 +1,115 @@ +/* -*-c-*- + * + * $Id: rabin.h,v 1.1 1999/11/19 13:17:57 mdw Exp $ + * + * Miller-Rabin primality test + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: rabin.h,v $ + * Revision 1.1 1999/11/19 13:17:57 mdw + * Prime number generator and tester. + * + */ + +#ifndef RABIN_H +#define RABIN_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MP_H +# include "mp.h" +#endif + +#ifndef MPMONT_H +# include "mpmont.h" +#endif + +#ifndef PGEN_H +# include "pgen.h" +#endif + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct rabin { + mpmont mm; /* Montgomery arithmetic context */ + size_t s; /* %$m = 2^s r + 1$% */ + mp *r; /* %$m = 2^s r + 1$% */ + mp *m1; /* %$(m - 1)R \bmod m */ +} rabin; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @rabin_create@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * @mp *m@ = pointer to number to test + * + * Returns: --- + * + * Use: Precomputes some useful values for performing the + * Miller-Rabin probabilistic primality test. + */ + +extern void rabin_create(rabin */*r*/, mp */*m*/); + +/* --- @rabin_destroy@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * + * Returns: --- + * + * Use: Disposes of a Rabin-Miller context when it's no longer + * needed. + */ + +extern void rabin_destroy(rabin */*r*/); + +/* --- @rabin_test@ --- * + * + * Arguments: @rabin *r@ = pointer to Rabin-Miller context + * @mp *g@ = base to test the number against + * + * Returns: Either @PGEN_COMPOSITE@ if the test failed, or @PGEN_MAYBE@ + * if it succeeded. + * + * Use: Performs a single iteration of the Rabin-Miller primality + * test. + */ + +extern int rabin_test(rabin */*r*/, mp */*g*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif -- 2.11.0