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 \
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 \
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',
*
* 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.
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);
}
{
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;
*
* 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.
* isn't finalized.
*/
-void gdsa_endhash(gdsa *c, ghash *h) { ; }
+void gdsa_endhash(const gdsa *c, ghash *h) { ; }
/* --- @gdsa_sign@ --- *
*
/* -*-c-*-
*
- * $Id: gdsa.h,v 1.2 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Generalized version of DSA
*
* isn't finalized.
*/
-extern void gdsa_endhash(gdsa */*c*/, ghash */*h*/);
+extern void gdsa_endhash(const gdsa */*c*/, ghash */*h*/);
/* --- @gdsa_sign@ --- *
*
/* -*-c-*-
*
- * $Id: gf-arith.c,v 1.4 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Basic arithmetic on binary polynomials
*
RIG(add, gf_add)
RIG(mul, gf_mul)
+RIG(exp, gf_exp)
#undef RIG
{ "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 } },
};
--- /dev/null
+/* -*-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 <assert.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
/* -*-c-*-
*
- * $Id: gf.h,v 1.4 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Arithmetic on binary polynomials
*
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
/*----- 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.
*
* isn't finalized.
*/
-void gkcdsa_endhash(gkcdsa *c, ghash *h) { ; }
+void gkcdsa_endhash(const gkcdsa *c, ghash *h) { ; }
/* --- @hashge@ --- *
*
/* -*-c-*-
*
- * $Id: gkcdsa.h,v 1.2 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Generalized version of KCDSA
*
/*----- 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.
*
* isn't finalized.
*/
-extern void gkcdsa_endhash(gkcdsa */*c*/, ghash */*h*/);
+extern void gkcdsa_endhash(const gkcdsa */*c*/, ghash */*h*/);
/* --- @gkcdsa_sign@ --- *
*
RIG(add, mp_add)
RIG(sub, mp_sub)
RIG(mul, mp_mul)
+RIG(exp, mp_exp)
#undef RIG
{ "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 } },
--- /dev/null
+/* -*-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 <assert.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
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
/* -*-c-*-
*
- * $Id: mpbarrett.c,v 1.10 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Barrett modular reduction
*
* @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$% --- */
b->vl[-1] = 1;
mp_div(&b, 0, b, m);
mb->mu = b;
+ return (0);
}
/* --- @mpbarrett_destroy@ --- *
/* -*-c-*-
*
- * $Id: mpbarrett.h,v 1.5 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Barrett modular reduction
*
* 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@ --- *
*
* 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.
#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);
/* --- 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$% --- */
mp_div(0, &mm->r2, r2, m);
mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
MP_DROP(r2);
+ return (0);
}
#endif
/* -*-c-*-
*
- * $Id: mpmont.h,v 1.8 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Montgomery reduction
*
* 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@ --- *
*
* 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 };
/* --- 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;
}
}
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
r->iv[i + r->in].argy = b;
}
}
- DA_DESTROY(&iv);
+ DA_DESTROY(&iv);
+ return (0);
}
/* --- @mpreduce_destroy@ --- *
/* -*-c-*-
*
- * $Id: mpreduce.h,v 1.2 2004/04/08 01:36:15 mdw Exp $
+ * $Id$
*
* Efficient reduction modulo nice primes
*
* 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@ --- *
*
-# $Id: gf,v 1.3 2004/03/27 17:54:12 mdw Exp $
+# $Id$
#
# Test cases for higher-level binary poly arithmetic.
0x398c4111da6d06cdf3d83704ee403101;
}
+exp {
+ 4 0 1;
+ 4 1 4;
+ 0x7 2 0x15;
+ 3 563 0xf000f000f000f0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000f000f000f000f;
+}
+
irred {
0xc1a7bd3b4e853fc92d4e1588719986aa 0;
0x800000000000000000000000000000000000000c9 1;
# 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;
2 1;
}
+exp {
+ 4 0 1;
+ 4 1 4;
+ 7 2 49;
+ 3 564 124849745640593184256214502788000232711984346194239284918599169775251467106591187580476305077269760425019686159071753053924227569816588462643229463821875763427430576080998505780547826368760514503807579784278708008217584939464444237989070811887584423210788916656247499281;
+}
+
bin2c {
and 5 3 1;
or 5 3 7;
--- /dev/null
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/alloc.h>
+#include <mLib/darray.h>
+#include <mLib/dstr.h>
+#include <mLib/quis.h>
+#include <mLib/report.h>
+
+#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);
+}