From f4535c6454395e6d56ce0091a07b6d4f7d54a47f Mon Sep 17 00:00:00 2001 From: mdw Date: Tue, 9 Nov 2004 11:26:04 +0000 Subject: [PATCH] A variety of small tweaks and fixes. Make mpmont etc. return errors rather than exploding messily. Add program for finding primitive polynomials (includes a poor version of ECM factoring!). Add exponentiation for integers and binary polynomials. --- Makefile.m4 | 8 +- f-niceprime.c | 7 +- f-prime.c | 7 +- field.h | 2 +- gdsa.c | 2 +- gdsa.h | 4 +- gf-arith.c | 4 +- gf-exp.c | 67 ++++++++++ gf-exp.h | 69 ++++++++++ gf.h | 13 +- gkcdsa.c | 6 +- gkcdsa.h | 8 +- mp-arith.c | 2 + mp-exp.c | 67 ++++++++++ mp-exp.h | 69 ++++++++++ mp.h | 11 ++ mpbarrett.c | 10 +- mpbarrett.h | 7 +- mpmont.c | 11 +- mpmont.h | 6 +- mpreduce.c | 16 ++- mpreduce.h | 6 +- tests/gf | 9 +- tests/mp | 9 +- utils/prim.c | 412 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 25 files changed, 783 insertions(+), 49 deletions(-) create mode 100644 gf-exp.c create mode 100644 gf-exp.h create mode 100644 mp-exp.c create mode 100644 mp-exp.h create mode 100644 utils/prim.c diff --git a/Makefile.m4 b/Makefile.m4 index 0fac1c4..b269aee 100644 --- a/Makefile.m4 +++ b/Makefile.m4 @@ -156,9 +156,9 @@ pkginclude_HEADERS = \ lcrand.h fibrand.h rc4.h seal.h rand.h noise.h fipstest.h maurer.h \ key.h key-data.h passphrase.h pixie.h lmem.h \ mpx.h bitops.h mpw.h mpscan.h mparena.h mp.h mptext.h mpint.h \ - exp.h mpbarrett.h mpmont.h mpreduce.h \ + exp.h mpbarrett.h mpmont.h mpreduce.h mp-exp.h \ mpcrt.h mprand.h mpmul.h \ - gfx.h gf.h gfreduce.h gfn.h \ + gfx.h gf.h gfreduce.h gfn.h gf-exp.h \ primetab.h pfilt.h rabin.h \ pgen.h prim.h strongprime.h limlee.h keycheck.h \ bbs.h rsa.h dh.h dsarand.h dsa.h gdsa.h gkcdsa.h \ @@ -176,7 +176,7 @@ pkginclude_HEADERS = \ define(`MP_SOURCES', `mpx.c mpx-kmul.c mpx-ksqr.c mpscan.c mparena.c qdparse.c \ mp-misc.c mp-mem.c mp-const.c mp-io.c mp-arith.c mp-test.c \ - mp-sqrt.c mp-gcd.c mp-jacobi.c mp-modsqrt.c \ + mp-sqrt.c mp-gcd.c mp-jacobi.c mp-modsqrt.c mp-exp.c \ mpint.c mptext.c mptext-file.c mptext-string.c mptext-dstr.c \ mptext-len.c \ exp.c mpcrt.c mpmul.c mprand.c \ @@ -190,7 +190,7 @@ define(`MP_SOURCES', GF_SOURCES PGEN_SOURCES EC_SOURCES') define(`GF_SOURCES', - `gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-gcd.c \ + `gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-exp.c gf-gcd.c \ gfreduce.c gfreduce-exp.h gfn.c') define(`EC_SOURCES', diff --git a/f-niceprime.c b/f-niceprime.c index a6191da..ce96977 100644 --- a/f-niceprime.c +++ b/f-niceprime.c @@ -138,7 +138,7 @@ static const field_ops fops = { * * Arguments: @mp *p@ = the characteristic of the field * - * Returns: A pointer to the field. + * Returns: A pointer to the field, or null. * * Use: Creates a field structure for a prime field of size %$p$%, * using efficient reduction for nice primes. @@ -152,7 +152,10 @@ field *field_niceprime(mp *p) f->f.one = MP_ONE; f->f.nbits = mp_bits(p); f->f.noctets = (f->f.nbits + 7) >> 3; - mpreduce_create(&f->r, p); + if (mpreduce_create(&f->r, p)) { + DESTROY(f); + return (0); + } f->f.m = f->r.p; return (&f->f); } diff --git a/f-prime.c b/f-prime.c index b1c4367..f2fd7eb 100644 --- a/f-prime.c +++ b/f-prime.c @@ -160,11 +160,12 @@ field *field_prime(mp *p) { fctx_prime *f; - if (!MP_POSP(p) || !MP_ODDP(p)) - return (0); f = CREATE(fctx_prime); f->f.ops = &fops; - mpmont_create(&f->mm, p); + if (mpmont_create(&f->mm, p)) { + DESTROY(f); + return (0); + } f->f.zero = MP_ZERO; f->f.one = f->mm.r; f->f.m = f->mm.m; diff --git a/field.h b/field.h index 4164bce..c2ad26e 100644 --- a/field.h +++ b/field.h @@ -208,7 +208,7 @@ extern field *field_prime(mp */*p*/); * * Arguments: @mp *p@ = the characteristic of the field * - * Returns: A pointer to the field. + * Returns: A pointer to the field, or null. * * Use: Creates a field structure for a prime field of size %$p$%, * using efficient reduction for nice primes. diff --git a/gdsa.c b/gdsa.c index 751bcce..f60fb2d 100644 --- a/gdsa.c +++ b/gdsa.c @@ -61,7 +61,7 @@ ghash *gdsa_beginhash(const gdsa *c) { return (GH_INIT(c->h)); } * isn't finalized. */ -void gdsa_endhash(gdsa *c, ghash *h) { ; } +void gdsa_endhash(const gdsa *c, ghash *h) { ; } /* --- @gdsa_sign@ --- * * diff --git a/gdsa.h b/gdsa.h index 479ef3c..660a465 100644 --- a/gdsa.h +++ b/gdsa.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gdsa.h,v 1.2 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Generalized version of DSA * @@ -91,7 +91,7 @@ extern ghash *gdsa_beginhash(const gdsa */*c*/); * isn't finalized. */ -extern void gdsa_endhash(gdsa */*c*/, ghash */*h*/); +extern void gdsa_endhash(const gdsa */*c*/, ghash */*h*/); /* --- @gdsa_sign@ --- * * diff --git a/gf-arith.c b/gf-arith.c index c23aa19..f380a35 100644 --- a/gf-arith.c +++ b/gf-arith.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gf-arith.c,v 1.4 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Basic arithmetic on binary polynomials * @@ -230,6 +230,7 @@ static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b) RIG(add, gf_add) RIG(mul, gf_mul) +RIG(exp, gf_exp) #undef RIG @@ -285,6 +286,7 @@ static test_chunk tests[] = { { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } }, { "sqr", tsqr, { &type_mp, &type_mp, 0 } }, { "div", tdiv, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { "exp", texp, { &type_mp, &type_mp, &type_mp, 0 } }, { "irred", tirred, { &type_mp, &type_int, 0 } }, { 0, 0, { 0 } }, }; diff --git a/gf-exp.c b/gf-exp.c new file mode 100644 index 0000000..61e2836 --- /dev/null +++ b/gf-exp.c @@ -0,0 +1,67 @@ +/* -*-c-*- + * + * $Id$ + * + * Exponentiation for binary polynomials + * + * (c) 2004 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. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "gf.h" +#include "gf-exp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gf_exp@ --- * + * + * Arguments: @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e$%. + */ + +mp *gf_exp(mp *d, mp *a, mp *e) +{ + mp *x = MP_ONE; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + assert(!MP_NEGP(e)); + + MP_COPY(a); + if (MP_ZEROP(e)) + ; + else if (MP_LEN(e) < EXP_THRESH) + EXP_SIMPLE(x, a, e); + else + EXP_WINDOW(x, a, e); + mp_drop(d); + mp_drop(spare); + mp_drop(a); + return (x); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gf-exp.h b/gf-exp.h new file mode 100644 index 0000000..5130bd5 --- /dev/null +++ b/gf-exp.h @@ -0,0 +1,69 @@ +/* -*-c-*- + * + * $Id$ + * + * Exponentiation for binary polynomials + * + * (c) 2004 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. + */ + +#ifndef CATACOMB_GF_EXP_H +#define CATACOMB_GF_EXP_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Exponentiation ----------------------------------------------------*/ + +#define EXP_TYPE mp * + +#define EXP_COPY(d, x) d = MP_COPY(x) +#define EXP_DROP(x) MP_DROP(x) + +#define EXP_MUL(a, x) do { \ + mp *t = gf_mul(spare, a, x); \ + spare = a; \ + a = t; \ +} while (0) + +#define EXP_SQR(a) do { \ + mp *t = gf_sqr(spare, a); \ + spare = a; \ + a = t; \ +} while (0) + +#define EXP_FIX(x) + +#define EXP_SETMUL(d, x, y) d = gf_mul(MP_NEW, x, y) +#define EXP_SETSQR(d, x) d = gf_sqr(MP_NEW, x) + +#include "exp.h" + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/gf.h b/gf.h index 2429e95..eea204f 100644 --- a/gf.h +++ b/gf.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gf.h,v 1.4 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Arithmetic on binary polynomials * @@ -91,6 +91,17 @@ extern mp *gf_sqr(mp */*d*/, mp */*a*/); extern void gf_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/); +/* --- @gf_exp@ --- * + * + * Arguments: @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e$%. + */ + +extern mp *gf_exp(mp */*d*/, mp */*a*/, mp */*e*/); + /* --- @gf_irreduciblep@ --- * * * Arguments: @mp *f@ = a polynomial diff --git a/gkcdsa.c b/gkcdsa.c index 433a158..a67ab0d 100644 --- a/gkcdsa.c +++ b/gkcdsa.c @@ -37,9 +37,9 @@ /*----- Main code ---------------------------------------------------------*/ -/* --- @gdsa_beginhash@ --- * +/* --- @gkcdsa_beginhash@ --- * * - * Arguments: @const gdsa *c@ = pointer to the context structure + * Arguments: @const gkcdsa *c@ = pointer to the context structure * * Returns: A hashing context for you to hash the message. * @@ -72,7 +72,7 @@ ghash *gkcdsa_beginhash(const gkcdsa *c) * isn't finalized. */ -void gkcdsa_endhash(gkcdsa *c, ghash *h) { ; } +void gkcdsa_endhash(const gkcdsa *c, ghash *h) { ; } /* --- @hashge@ --- * * diff --git a/gkcdsa.h b/gkcdsa.h index ec27eab..5392158 100644 --- a/gkcdsa.h +++ b/gkcdsa.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gkcdsa.h,v 1.2 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Generalized version of KCDSA * @@ -67,9 +67,9 @@ typedef struct gkcdsa_sig { /*----- Functions provided ------------------------------------------------*/ -/* --- @gdsa_beginhash@ --- * +/* --- @gkcdsa_beginhash@ --- * * - * Arguments: @const gdsa *c@ = pointer to the context structure + * Arguments: @const gkcdsa *c@ = pointer to the context structure * * Returns: A hashing context for you to hash the message. * @@ -91,7 +91,7 @@ extern ghash *gkcdsa_beginhash(const gkcdsa */*c*/); * isn't finalized. */ -extern void gkcdsa_endhash(gkcdsa */*c*/, ghash */*h*/); +extern void gkcdsa_endhash(const gkcdsa */*c*/, ghash */*h*/); /* --- @gkcdsa_sign@ --- * * diff --git a/mp-arith.c b/mp-arith.c index 4bd0976..0172981 100644 --- a/mp-arith.c +++ b/mp-arith.c @@ -723,6 +723,7 @@ RIG(lsr2c, mp_lsr2c) RIG(add, mp_add) RIG(sub, mp_sub) RIG(mul, mp_mul) +RIG(exp, mp_exp) #undef RIG @@ -908,6 +909,7 @@ static test_chunk tests[] = { { "sub", tsub, { &type_mp, &type_mp, &type_mp, 0 } }, { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } }, { "div", tdiv, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { "exp", texp, { &type_mp, &type_mp, &type_mp, 0 } }, { "bin2c", tbin, { &type_string, &type_mp, &type_mp, &type_mp, 0 } }, { "odd", todd, { &type_mp, &type_uint32, &type_mp, 0 } }, { "neg", tneg, { &type_mp, &type_mp, 0 } }, diff --git a/mp-exp.c b/mp-exp.c new file mode 100644 index 0000000..d4a45e2 --- /dev/null +++ b/mp-exp.c @@ -0,0 +1,67 @@ +/* -*-c-*- + * + * $Id$ + * + * Exponentiation for large integers + * + * (c) 2004 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. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "mp.h" +#include "mp-exp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_exp@ --- * + * + * Arguments: @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e$%. + */ + +mp *mp_exp(mp *d, mp *a, mp *e) +{ + mp *x = MP_ONE; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + assert(!MP_NEGP(e)); + + MP_COPY(a); + if (MP_ZEROP(e)) + ; + else if (MP_LEN(e) < EXP_THRESH) + EXP_SIMPLE(x, a, e); + else + EXP_WINDOW(x, a, e); + mp_drop(d); + mp_drop(spare); + mp_drop(a); + return (x); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-exp.h b/mp-exp.h new file mode 100644 index 0000000..13f5d28 --- /dev/null +++ b/mp-exp.h @@ -0,0 +1,69 @@ +/* -*-c-*- + * + * $Id$ + * + * Exponentiation for large integers + * + * (c) 2004 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. + */ + +#ifndef CATACOMB_MP_EXP_H +#define CATACOMB_MP_EXP_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Exponentiation ----------------------------------------------------*/ + +#define EXP_TYPE mp * + +#define EXP_COPY(d, x) d = MP_COPY(x) +#define EXP_DROP(x) MP_DROP(x) + +#define EXP_MUL(a, x) do { \ + mp *t = mp_mul(spare, a, x); \ + spare = a; \ + a = t; \ +} while (0) + +#define EXP_SQR(a) do { \ + mp *t = mp_sqr(spare, a); \ + spare = a; \ + a = t; \ +} while (0) + +#define EXP_FIX(x) + +#define EXP_SETMUL(d, x, y) d = mp_mul(MP_NEW, x, y) +#define EXP_SETSQR(d, x) d = mp_sqr(MP_NEW, x) + +#include "exp.h" + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/mp.h b/mp.h index f6fa14f..aa9bbac 100644 --- a/mp.h +++ b/mp.h @@ -866,6 +866,17 @@ extern mp *mp_sqr(mp */*d*/, mp */*a*/); extern void mp_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/); +/* --- @mp_exp@ --- * + * + * Arguments: @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e$%. + */ + +extern mp *mp_exp(mp */*d*/, mp */*a*/, mp */*e*/); + /* --- @mp_odd@ --- * * * Arguments: @mp *d@ = pointer to destination integer diff --git a/mpbarrett.c b/mpbarrett.c index ffd2649..b5e604e 100644 --- a/mpbarrett.c +++ b/mpbarrett.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpbarrett.c,v 1.10 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Barrett modular reduction * @@ -40,18 +40,19 @@ * @mp *m@ = modulus to work to * * - * Returns: --- + * Returns: Zero on success, nonzero on error. * * Use: Initializes a Barrett reduction context ready for use. */ -void mpbarrett_create(mpbarrett *mb, mp *m) +int mpbarrett_create(mpbarrett *mb, mp *m) { mp *b; /* --- Validate the arguments --- */ - assert(((void)"Barrett modulus must be positive", (m->f & MP_NEG) == 0)); + if (!MP_POSP(m)) + return (-1); /* --- Compute %$\mu$% --- */ @@ -63,6 +64,7 @@ void mpbarrett_create(mpbarrett *mb, mp *m) b->vl[-1] = 1; mp_div(&b, 0, b, m); mb->mu = b; + return (0); } /* --- @mpbarrett_destroy@ --- * diff --git a/mpbarrett.h b/mpbarrett.h index 85c4a14..bb8b36a 100644 --- a/mpbarrett.h +++ b/mpbarrett.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpbarrett.h,v 1.5 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Barrett modular reduction * @@ -72,13 +72,12 @@ typedef struct mpbarrett { * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context * @mp *m@ = modulus to work to * - * - * Returns: --- + * Returns: Zero on success, nonzero on error. * * Use: Initializes a Barrett reduction context ready for use. */ -extern void mpbarrett_create(mpbarrett */*mb*/, mp */*m*/); +extern int mpbarrett_create(mpbarrett */*mb*/, mp */*m*/); /* --- @mpbarrett_destroy@ --- * * diff --git a/mpmont.c b/mpmont.c index 24fc9a3..e62678d 100644 --- a/mpmont.c +++ b/mpmont.c @@ -49,7 +49,7 @@ * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context * @mp *m@ = modulus to use * - * Returns: --- + * Returns: Zero on success, nonzero on error. * * Use: Initializes a Montgomery reduction context ready for use. * The argument @m@ must be a positive odd integer. @@ -57,18 +57,19 @@ #ifdef MPMONT_DISABLE -void mpmont_create(mpmont *mm, mp *m) +int mpmont_create(mpmont *mm, mp *m) { mp_shrink(m); mm->m = MP_COPY(m); mm->r = MP_ONE; mm->r2 = MP_ONE; mm->mi = MP_ONE; + return (0); } #else -void mpmont_create(mpmont *mm, mp *m) +int mpmont_create(mpmont *mm, mp *m) { size_t n = MP_LEN(m); mp *r2 = mp_new(2 * n + 1, 0); @@ -76,7 +77,8 @@ void mpmont_create(mpmont *mm, mp *m) /* --- Take a copy of the modulus --- */ - assert(MP_POSP(m) && MP_ODDP(m)); + if (!MP_POSP(m) || !MP_ODDP(m)) + return (-1); mm->m = MP_COPY(m); /* --- Determine %$R^2$% --- */ @@ -97,6 +99,7 @@ void mpmont_create(mpmont *mm, mp *m) mp_div(0, &mm->r2, r2, m); mm->r = mpmont_reduce(mm, MP_NEW, mm->r2); MP_DROP(r2); + return (0); } #endif diff --git a/mpmont.h b/mpmont.h index b931162..44be367 100644 --- a/mpmont.h +++ b/mpmont.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpmont.h,v 1.8 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Montgomery reduction * @@ -92,13 +92,13 @@ typedef struct mpmont { * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context * @mp *m@ = modulus to use * - * Returns: --- + * Returns: Zero on success, nonzero on error. * * Use: Initializes a Montgomery reduction context ready for use. * The argument @m@ must be a positive odd integer. */ -extern void mpmont_create(mpmont */*mm*/, mp */*m*/); +extern int mpmont_create(mpmont */*mm*/, mp */*m*/); /* --- @mpmont_destroy@ --- * * diff --git a/mpreduce.c b/mpreduce.c index 13e705e..d99ff16 100644 --- a/mpreduce.c +++ b/mpreduce.c @@ -47,12 +47,12 @@ DA_DECL(instr_v, mpreduce_instr); * Arguments: @gfreduce *r@ = structure to fill in * @mp *x@ = an integer * - * Returns: --- + * Returns: Zero if successful; nonzero on failure. * * Use: Initializes a context structure for reduction. */ -void mpreduce_create(mpreduce *r, mp *p) +int mpreduce_create(mpreduce *r, mp *p) { mpscan sc; enum { Z = 0, Z1 = 2, X = 4, X0 = 6 }; @@ -64,7 +64,8 @@ void mpreduce_create(mpreduce *r, mp *p) /* --- Fill in the easy stuff --- */ - assert(MP_POSP(p)); + if (!MP_POSP(p)) + return (-1); d = mp_bits(p); r->lim = d/MPW_BITS; r->s = d%MPW_BITS; @@ -118,9 +119,9 @@ void mpreduce_create(mpreduce *r, mp *p) } } if ((DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) { - fprintf(stderr, - "mpreduce can't handle moduli of the form 2^m + x\n"); - abort(); + mp_drop(r->p); + DA_DESTROY(&iv); + return (-1); } #undef INSTR @@ -149,7 +150,8 @@ void mpreduce_create(mpreduce *r, mp *p) r->iv[i + r->in].argy = b; } } - DA_DESTROY(&iv); + DA_DESTROY(&iv); + return (0); } /* --- @mpreduce_destroy@ --- * diff --git a/mpreduce.h b/mpreduce.h index c885ccb..48c9e1d 100644 --- a/mpreduce.h +++ b/mpreduce.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpreduce.h,v 1.2 2004/04/08 01:36:15 mdw Exp $ + * $Id$ * * Efficient reduction modulo nice primes * @@ -66,12 +66,12 @@ typedef struct mpreduce { * Arguments: @gfreduce *r@ = structure to fill in * @mp *x@ = an integer * - * Returns: --- + * Returns: Zero for success, nonzero on error. * * Use: Initializes a context structure for reduction. */ -extern void mpreduce_create(mpreduce */*r*/, mp */*p*/); +extern int mpreduce_create(mpreduce */*r*/, mp */*p*/); /* --- @mpreduce_destroy@ --- * * diff --git a/tests/gf b/tests/gf index df5b793..801993d 100644 --- a/tests/gf +++ b/tests/gf @@ -1,4 +1,4 @@ -# $Id: gf,v 1.3 2004/03/27 17:54:12 mdw Exp $ +# $Id$ # # Test cases for higher-level binary poly arithmetic. @@ -51,6 +51,13 @@ div { 0x398c4111da6d06cdf3d83704ee403101; } +exp { + 4 0 1; + 4 1 4; + 0x7 2 0x15; + 3 563 0xf000f000f000f0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000f000f000f000f; +} + irred { 0xc1a7bd3b4e853fc92d4e1588719986aa 0; 0x800000000000000000000000000000000000000c9 1; diff --git a/tests/mp b/tests/mp index 34f16f8..c8be382 100644 --- a/tests/mp +++ b/tests/mp @@ -1,6 +1,6 @@ # Test vectors for MP functions # -# $Id: mp,v 1.17 2004/04/01 12:50:41 mdw Exp $ +# $Id$ add { 5 4 9; 5 -4 1; -5 4 -1; -5 -4 -9; @@ -32,6 +32,13 @@ div { 2 1; } +exp { + 4 0 1; + 4 1 4; + 7 2 49; + 3 564 124849745640593184256214502788000232711984346194239284918599169775251467106591187580476305077269760425019686159071753053924227569816588462643229463821875763427430576080998505780547826368760514503807579784278708008217584939464444237989070811887584423210788916656247499281; +} + bin2c { and 5 3 1; or 5 3 7; diff --git a/utils/prim.c b/utils/prim.c new file mode 100644 index 0000000..b279304 --- /dev/null +++ b/utils/prim.c @@ -0,0 +1,412 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "fibrand.h" +#include "field.h" +#include "ec.h" +#include "gf.h" +#include "gfreduce.h" +#include "mp.h" +#include "mpbarrett.h" +#include "mpint.h" +#include "mprand.h" +#include "pgen.h" +#include "primetab.h" + +DA_DECL(mp_v, mp *); + +typedef struct fbe { + unsigned long p; + unsigned e; + mp *n; +} fbe; +DA_DECL(fbe_v, fbe); + +typedef struct fact { + mp *p; + unsigned e; + mp *n; + mp *r; +} fact; +DA_DECL(fact_v, fact); + +static grand *rng; + +/*----- Factoring ---------------------------------------------------------*/ + +static void mpv_free(mp_v *v) +{ + size_t i; + + for (i = 0; i < DA_LEN(v); i++) + mp_drop(DA(v)[i]); + DA_DESTROY(v); +} + +static mp *smallfactors(mp *x, mp_v *v) +{ + mp f; + mpw fw; + mp *y = MP_NEW, *z = MP_NEW; + unsigned i; + + MP_COPY(x); + mp_build(&f, &fw, &fw + 1); +#ifdef DEBUG + MP_PRINT("target", x); +#endif + for (i = 0; i < NPRIME; i++) { + again: + fw = primetab[i]; + mp_div(&y, &z, x, &f); + if (MP_ZEROP(z)) { +#ifdef DEBUG + MP_PRINT("factor", &f); +#endif + MP_DROP(x); + x = y; + y = MP_NEW; +#ifdef DEBUG + MP_PRINT("remainder", x); +#endif + DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i])); + goto again; + } + } + mp_drop(y); + mp_drop(z); + return (x); +} + +#ifdef POLLARDRHO +static mp *pollardrho(mp *x) +{ + unsigned i; + mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW; + mpbarrett mb; + + mpbarrett_create(&mb, x); + y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y); + for (i = 0; i < 1024; i++) { + y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y); + z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); + z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z); + t = mp_sub(t, y, z); + mp_gcd(&g, 0, 0, x, t); + if (!MP_EQ(g, MP_ONE)) + goto done; + } + MP_DROP(g); + g = 0; + +done: + mpbarrett_destroy(&mb); + MP_DROP(y); + MP_DROP(z); + MP_DROP(t); +#ifdef DEBUG + MP_PRINT("rho factor", g); +#endif + return (g); +} +#endif + +static void addfbe(fbe_v *v, mp *B, unsigned long p) +{ + fbe e; + unsigned w; + mp *pp = mp_fromulong(MP_NEW, p); + mp *qq = MP_NEW, *rr = MP_NEW; + mp *t; + + w = 0; + qq = MP_COPY(pp); + for (;;) { + w++; + rr = mp_mul(rr, qq, pp); + if (MP_CMP(rr, >=, B)) break; + t = rr; rr = qq; qq = t; + } + e.p = p; + e.e = w; + e.n = qq; + MP_DROP(rr); + MP_DROP(pp); + DA_PUSH(v, e); +} + +static void makefactorbase(mp *x, fbe_v *v) +{ + mp *B; + unsigned nb; + unsigned long max; + unsigned long i; + size_t j; + + nb = sqrt(mp_bits(x)) * 1.5; + B = mp_lsl(MP_NEW, MP_ONE, nb); + if (nb < 32) + max = 1 << nb; + else + max = 0xfffffff0; + + addfbe(v, B, 2); + for (i = 3; i < max; i++) { + for (j = 0; j < DA_LEN(v); j++) { + if (i % DA(v)[j].p == 0) + goto composite; + } + addfbe(v, B, i); + composite:; + } + +#ifdef DEBUG + for (j = 0; j < DA_LEN(v); j++) { + printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e); + mp_writefile(DA(v)[j].n, stdout, 10); + putchar('\n'); + } +#endif +} + +static void freefactorbase(fbe_v *v) +{ + size_t i; + + for (i = 0; i < DA_LEN(v); i++) + mp_drop(DA(v)[i].n); + DA_DESTROY(v); +} + +static mp *ecm(mp *x) +{ + fbe_v fb = DA_INIT; + field *f = field_prime(x); + ec_curve *c; + ec p; + mp *a = MP_NEW, *b = MP_NEW; + mp *g = MP_NEW; + unsigned i; + size_t j; + + makefactorbase(x, &fb); + for (;;) { + a = mprand_range(a, x, rng, 0); + b = mprand_range(b, x, rng, 0); + c = ec_primeproj(f, a, b); + for (i = 0; i < 1024; i++) { + EC_CREATE(&p); + p.x = mprand_range(p.x, x, rng, 0); + p.y = mprand_range(p.y, x, rng, 0); + p.z = MP_ONE; + for (j = 0; j < DA_LEN(&fb); j++) + ec_imul(c, &p, &p, DA(&fb)[j].n); + if (EC_ATINF(&p)) + continue; + mp_gcd(&g, 0, 0, x, p.z); + if (!MP_EQ(g, MP_ONE)) + goto done; + EC_DESTROY(&p); + } + ec_destroycurve(c); + } + +done: + MP_DROP(a); MP_DROP(b); + EC_DESTROY(&p); + ec_destroycurve(c); + F_DESTROY(f); + freefactorbase(&fb); +#ifdef DEBUG + MP_PRINT("ecm factor", g); +#endif + return (g); +} + +static void dofactor(mp *x, mp_v *v) +{ + mp *f; + + x = smallfactors(x, v); + if (MP_EQ(x, MP_ONE)) return; + +#ifdef POLLARDRHO + for (;;) { + if (pgen_primep(x, rng)) + goto done; + f = pollardrho(x); + if (!f) break; + dofactor(f, v); + mp_div(&x, 0, x, f); + MP_DROP(f); + } +#endif + + for (;;) { + if (pgen_primep(x, rng)) + goto done; + f = ecm(x); + if (!f) break; + dofactor(f, v); + mp_div(&x, 0, x, f); + MP_DROP(f); + } + +done: + DA_PUSH(v, x); +} + +static int cmpmp(const void *a, const void *b) +{ + mp *const *aa = a, *const *bb = b; + return (mp_cmp(*aa, *bb)); +} + +static void factor(mp *x, fact_v *v) +{ + fact f; + mp_v fv = DA_INIT; + mp *e = MP_NEW; + size_t i; + + dofactor(x, &fv); + qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp); + for (i = 0; i < DA_LEN(&fv); ) { + f.p = MP_COPY(DA(&fv)[i]); + for (f.e = 1, i++; + i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]); + f.e++, i++) ; + e = mp_fromuint(e, f.e); + f.n = mp_exp(MP_NEW, f.p, e); + f.r = MP_NEW; + mp_div(&f.r, 0, x, f.p); + DA_PUSH(v, f); + } + mp_drop(e); + mpv_free(&fv); +} + +static void freefactors(fact_v *v) +{ + size_t i; + + for (i = 0; i < DA_LEN(v); i++) { + mp_drop(DA(v)[i].p); + mp_drop(DA(v)[i].n); + mp_drop(DA(v)[i].r); + } + DA_DESTROY(v); +} + + +/*----- Primitive polynomials ---------------------------------------------*/ + +static int primitivep(fact_v *f, mp *p) +{ + gfreduce r; + unsigned i; + mp *x = MP_NEW; + int rc = 1; + + if (!gf_irreduciblep(p)) + return (0); +#ifdef DEBUG + MP_PRINTX("p", p); +#endif + gfreduce_create(&r, p); + for (i = 0; i < DA_LEN(f); i++) { + x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r); +#ifdef DEBUG + MP_PRINT(" r", DA(f)[i].r); + MP_PRINTX(" x", x); +#endif + if (MP_EQ(x, MP_ONE)) { + rc = 0; + break; + } + } + gfreduce_destroy(&r); + MP_DROP(x); + return (rc); +} + +static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x) +{ + unsigned i; + + for (i = n + 1; i < x; i++) { + *p = mp_setbit(*p, *p, i); + if (n ? dofip(f, m, p, n - 1, i) : primitivep(f, *p)) + return (1); + *p = mp_clearbit(*p, *p, i); + } + return (0); +} + +static mp *fip(fact_v *f, unsigned m) +{ + mp *p = MP_ONE; + unsigned n; + + p = mp_setbit(p, p, m); + n = 0; + while (!dofip(f, m, &p, n, m)) + n += 2; + return (p); +} + +static void findprim(unsigned m) +{ + fact_v f = DA_INIT; + mp *q = mp_lsl(MP_NEW, MP_ONE, m); + mp *p; + unsigned i; + + q = mp_sub(q, q, MP_ONE); + factor(q, &f); + +#ifdef DEBUG + { + size_t i; + for (i = 0; i < DA_LEN(&f); i++) { + mp_writefile(DA(&f)[i].p, stdout, 10); + printf("^%u = ", DA(&f)[i].e); + mp_writefile(DA(&f)[i].n, stdout, 10); + fputs(" (", stdout); + mp_writefile(DA(&f)[i].r, stdout, 10); + fputs(")\n", stdout); + } + } +#endif + + p = fip(&f, m); + MP_PRINTX("p", p); + for (i = m + 1; i--; ) { + if (mp_testbit(p, i)) + printf("%u ", i); + } + putchar('\n'); + mp_drop(p); + mp_drop(q); + freefactors(&f); +} + +int main(int argc, char *argv[]) +{ + ego(argv[0]); + rng = fibrand_create(0); + if (argc != 2) { + fprintf(stderr, "Usage: %s M\n", QUIS); + exit(1); + } + findprim(atoi(argv[1])); + return (0); +} -- 2.11.0