From f46efa79cd2bb9adc81541f1218965f85a6b2eef Mon Sep 17 00:00:00 2001
From: mdw <mdw>
Date: Sat, 27 Mar 2004 00:04:46 +0000
Subject: [PATCH] Implement efficient reduction for pleasant-looking primes.

---
 Makefile.m4    |  17 ++-
 ec-prime.c     |   7 +-
 ec-test.c      |  10 +-
 f-niceprime.c  | 210 ++++++++++++++++++++++++++++
 field.h        |  17 ++-
 gfreduce.c     |  12 +-
 mpreduce-exp.h |  90 ++++++++++++
 mpreduce.c     | 432 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 mpreduce.h     | 136 ++++++++++++++++++
 mpx.c          |  56 +++++++-
 mpx.h          |  38 ++++-
 tests/ec       | 150 +++++++++++++++++++-
 tests/mpreduce |  21 +++
 13 files changed, 1178 insertions(+), 18 deletions(-)
 create mode 100644 f-niceprime.c
 create mode 100644 mpreduce-exp.h
 create mode 100644 mpreduce.c
 create mode 100644 mpreduce.h
 create mode 100644 tests/mpreduce

diff --git a/Makefile.m4 b/Makefile.m4
index ee3872a..cc500de 100644
--- a/Makefile.m4
+++ b/Makefile.m4
@@ -1,6 +1,6 @@
 ## -*-m4-*-
 ##
-## $Id: Makefile.m4,v 1.69 2004/03/23 15:19:32 mdw Exp $
+## $Id: Makefile.m4,v 1.70 2004/03/27 00:04:46 mdw Exp $
 ##
 ## Makefile for Catacomb
 ##
@@ -29,6 +29,9 @@
 ##----- Revision history ----------------------------------------------------
 ##
 ## $Log: Makefile.m4,v $
+## Revision 1.70  2004/03/27 00:04:46  mdw
+## Implement efficient reduction for pleasant-looking primes.
+##
 ## Revision 1.69  2004/03/23 15:19:32  mdw
 ## Test elliptic curves more thoroughly.
 ##
@@ -346,9 +349,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 mpbarrett-exp.h mpmont.h mpmont-exp.h \
+	exp.h mpbarrett.h mpmont.h mpreduce.h \
 	mpcrt.h mprand.h mpmul.h \
-	gfx.h gf.h gfreduce.h gfreduce-exp.h \
+	gfx.h gf.h gfreduce.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 \
@@ -370,14 +373,17 @@ define(`MP_SOURCES',
 	exp.c mpcrt.c mpmul.c mprand.c \
 	mpbarrett.c mpbarrett-mexp.c mpbarrett-exp.h \
 	mpmont.c mpmont-mexp.c mpmont-exp.h \
+	mpreduce.c mpreduce-exp.h \
 	rho.c buf.c \
 	GF_SOURCES PGEN_SOURCES EC_SOURCES')
 
 define(`GF_SOURCES',
-	`gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-gcd.c gfreduce.c')
+	`gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-gcd.c \
+	gfreduce.c gfreduce-exp.h ')
 
 define(`EC_SOURCES',
-	`field.c f-prime.c f-binpoly.c ec.c ec-prime.c ec-bin.c ec-test.c')
+	`field.c f-prime.c f-niceprime.c f-binpoly.c \
+	ec.c ec-prime.c ec-bin.c ec-test.c')
 
 define(`PGEN_SOURCES',
 	`pfilt.c rabin.c \
@@ -543,6 +549,7 @@ CTESTRIG(mpbarrett)
 CTESTRIG(mpbarrett-mexp)
 CTESTRIG(mpmont)
 CTESTRIG(mpmont-mexp)
+CTESTRIG(mpreduce)
 CTESTRIG(mpcrt)
 CTESTRIG(mpmul)
 CTESTRIG(gfx)
diff --git a/ec-prime.c b/ec-prime.c
index 5c6297f..7712bbf 100644
--- a/ec-prime.c
+++ b/ec-prime.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: ec-prime.c,v 1.6 2004/03/23 15:19:32 mdw Exp $
+ * $Id: ec-prime.c,v 1.7 2004/03/27 00:04:46 mdw Exp $
  *
  * Elliptic curves over prime fields
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: ec-prime.c,v $
+ * Revision 1.7  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.6  2004/03/23 15:19:32  mdw
  * Test elliptic curves more thoroughly.
  *
@@ -464,7 +467,7 @@ int main(int argc, char *argv[])
   p = MP(6277101735386680763835789423207666416083908700390324961279);
   r = MP(6277101735386680763835789423176059013767194773182842284080);
 
-  f = field_prime(p);
+  f = field_niceprime(p);
   c = ec_primeproj(f, a, b);
   
   g.x = MP(0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012);
diff --git a/ec-test.c b/ec-test.c
index 6d42a70..59acf4e 100644
--- a/ec-test.c
+++ b/ec-test.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: ec-test.c,v 1.1 2004/03/23 15:19:32 mdw Exp $
+ * $Id: ec-test.c,v 1.2 2004/03/27 00:04:46 mdw Exp $
  *
  * Code for testing elliptic-curve stuff
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: ec-test.c,v $
+ * Revision 1.2  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.1  2004/03/23 15:19:32  mdw
  * Test elliptic curves more thoroughly.
  *
@@ -189,7 +192,7 @@ static void ecvcvt(const char *buf, dstr *d)
   int i;
 
   static const char *fnames[] = {
-    "prime", "binpoly", 0
+    "prime", "niceprime", "binpoly", 0
   };
   static const char *ecnames[] = {
     "prime", "primeproj", "bin", "binproj", 0
@@ -197,7 +200,8 @@ static void ecvcvt(const char *buf, dstr *d)
 
   switch (i = ckstring(&p, fnames), ckchar(&p, ':'), i) {
     case 0: m = getmp(&p); f = field_prime(m); mp_drop(m); break;
-    case 1: m = getmp(&p); f = field_binpoly(m); mp_drop(m); break;
+    case 1: m = getmp(&p); f = field_niceprime(m); mp_drop(m); break;
+    case 2: m = getmp(&p); f = field_binpoly(m); mp_drop(m); break;
     default: abort();
   }
   ckchar(&p, '/');
diff --git a/f-niceprime.c b/f-niceprime.c
new file mode 100644
index 0000000..c34e24e
--- /dev/null
+++ b/f-niceprime.c
@@ -0,0 +1,210 @@
+/* -*-c-*-
+ *
+ * $Id: f-niceprime.c,v 1.1 2004/03/27 00:04:46 mdw Exp $
+ *
+ * Prime fields with efficient reduction for special-form primes
+ *
+ * (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.
+ */
+
+/*----- Revision history --------------------------------------------------* 
+ *
+ * $Log: f-niceprime.c,v $
+ * Revision 1.1  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include <mLib/sub.h>
+
+#include "field.h"
+#include "mpreduce.h"
+#include "mprand.h"
+
+/*----- Data structures ---------------------------------------------------*/
+
+typedef struct fctx {
+  field f;
+  mpreduce r;
+} fctx;
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- Field operations --- */
+
+static void fdestroy(field *ff)
+{
+  fctx *f = (fctx *)ff;
+  mpreduce_destroy(&f->r);
+  DESTROY(f);
+}
+
+static mp *frand(field *ff, mp *d, grand *r)
+{
+  fctx *f = (fctx *)ff;
+  return (mprand_range(d, f->r.p, r, 0));
+}
+
+static int fzerop(field *ff, mp *x)
+{
+  return (!MP_LEN(x));
+}
+
+static mp *fneg(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  return (mp_sub(d, f->r.p, x));
+}
+
+static mp *fadd(field *ff, mp *d, mp *x, mp *y)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_add(d, x, y);
+  if (d->f & MP_NEG)
+    d = mp_add(d, d, f->r.p);
+  else if (MP_CMP(d, >, f->r.p))
+    d = mp_sub(d, d, f->r.p);
+  return (d);
+}
+
+static mp *fsub(field *ff, mp *d, mp *x, mp *y)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_sub(d, x, y);
+  if (d->f & MP_NEG)
+    d = mp_add(d, d, f->r.p);
+  else if (MP_CMP(d, >, f->r.p))
+    d = mp_sub(d, d, f->r.p);
+  return (d);
+}
+
+static mp *fmul(field *ff, mp *d, mp *x, mp *y)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_mul(d, x, y);
+  return (mpreduce_do(&f->r, d, d));
+}
+
+static mp *fsqr(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_sqr(d, x);
+  return (mpreduce_do(&f->r, d, d));
+}
+
+static mp *finv(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  mp_gcd(0, 0, &d, f->r.p, x);
+  return (d);
+}
+
+static mp *freduce(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  return (mpreduce_do(&f->r, d, x));
+}
+
+static mp *fsqrt(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  return (mp_modsqrt(d, x, f->r.p));
+}
+
+static mp *fdbl(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_lsl(d, x, 1);
+  if (MP_CMP(d, >, f->r.p))
+    d = mp_sub(d, d, f->r.p);
+  return (d);
+}
+
+static mp *ftpl(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  MP_DEST(d, MP_LEN(x) + 1, x->f);
+  MPX_UMULN(d->v, d->vl, x->v, x->vl, 3);
+  while (MP_CMP(d, >, f->r.p))
+    d = mp_sub(d, d, f->r.p);
+  return (d);
+}
+
+static mp *fqdl(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  d = mp_lsl(d, x, 2);
+  while (MP_CMP(d, >, f->r.p))
+    d = mp_sub(d, d, f->r.p);
+  return (d);
+}
+
+static mp *fhlv(field *ff, mp *d, mp *x)
+{
+  fctx *f = (fctx *)ff;
+  if (!MP_LEN(x)) {
+    MP_COPY(x);
+    MP_DROP(d);
+    return (x);
+  }
+  if (x->v[0] & 1) {
+    d = mp_add(d, x, f->r.p);
+    x = d;
+  }
+  return (mp_lsr(d, x, 1));
+}
+
+/* --- Field operations table --- */
+
+static field_ops fops = {
+  FTY_PRIME, "niceprime",
+  fdestroy, frand,
+  freduce, field_id,
+  fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, fsqrt,
+  0,
+  fdbl, ftpl, fqdl, fhlv
+};
+
+/* --- @field_niceprime@ --- *
+ *
+ * Arguments:	@mp *p@ = the characteristic of the field
+ *
+ * Returns:	A pointer to the field.
+ *
+ * Use:		Creates a field structure for a prime field of size %$p$%,
+ *		using efficient reduction for nice primes.
+ */
+
+field *field_niceprime(mp *p)
+{
+  fctx *f = CREATE(fctx);
+  f->f.ops = &fops;
+  f->f.zero = MP_ZERO;
+  f->f.one = MP_ONE;
+  mpreduce_create(&f->r, p);
+  return (&f->f);
+}
+
+/*----- That's all, folks -------------------------------------------------*/
diff --git a/field.h b/field.h
index db27e63..21ac7dd 100644
--- a/field.h
+++ b/field.h
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: field.h,v 1.6 2004/03/23 15:19:32 mdw Exp $
+ * $Id: field.h,v 1.7 2004/03/27 00:04:46 mdw Exp $
  *
  * Definitions for field arithmetic
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: field.h,v $
+ * Revision 1.7  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.6  2004/03/23 15:19:32  mdw
  * Test elliptic curves more thoroughly.
  *
@@ -179,6 +182,18 @@ extern mp *field_id(field */*f*/, mp */*d*/, mp */*x*/);
 
 extern field *field_prime(mp */*p*/);
 
+/* --- @field_niceprime@ --- *
+ *
+ * Arguments:	@mp *p@ = the characteristic of the field
+ *
+ * Returns:	A pointer to the field.
+ *
+ * Use:		Creates a field structure for a prime field of size %$p$%,
+ *		using efficient reduction for nice primes.
+ */
+
+extern field *field_niceprime(mp */*p*/);
+
 /* --- @field_binpoly@ --- *
  *
  * Arguments:	@mp *p@ = an irreducible polynomial over %$\gf{2}$%
diff --git a/gfreduce.c b/gfreduce.c
index 929c46c..4f07ccf 100644
--- a/gfreduce.c
+++ b/gfreduce.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: gfreduce.c,v 1.3 2004/03/23 15:19:32 mdw Exp $
+ * $Id: gfreduce.c,v 1.4 2004/03/27 00:04:46 mdw Exp $
  *
  * Efficient reduction modulo sparse binary polynomials
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: gfreduce.c,v $
+ * Revision 1.4  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.3  2004/03/23 15:19:32  mdw
  * Test elliptic curves more thoroughly.
  *
@@ -98,7 +101,8 @@ DA_DECL(instr_v, gfreduce_instr);
 void gfreduce_create(gfreduce *r, mp *p)
 {
   instr_v iv = DA_INIT;
-  unsigned long d, dw;
+  unsigned long d;
+  unsigned dw;
   mpscan sc;
   unsigned long i;
   gfreduce_instr *ip;
@@ -162,8 +166,8 @@ void gfreduce_create(gfreduce *r, mp *p)
       w = ww;
       wi = DA_LEN(&iv);
     }
-    INSTR(GFRI_LSL, (i - d)%MPW_BITS);
-    if ((i - d)%MPW_BITS)
+    INSTR(GFRI_LSL, (MPW_BITS + i - d)%MPW_BITS);
+    if ((MPW_BITS + i - d)%MPW_BITS)
       f |= f_lsr;
   }
   wl = DA_LEN(&iv);
diff --git a/mpreduce-exp.h b/mpreduce-exp.h
new file mode 100644
index 0000000..781cfab
--- /dev/null
+++ b/mpreduce-exp.h
@@ -0,0 +1,90 @@
+/* -*-c-*-
+ *
+ * $Id: mpreduce-exp.h,v 1.1 2004/03/27 00:04:46 mdw Exp $
+ *
+ * Exponentiation operations for binary field reduction
+ *
+ * (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.
+ */
+
+/*----- Revision history --------------------------------------------------* 
+ *
+ * $Log: mpreduce-exp.h,v $
+ * Revision 1.1  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
+ * Revision 1.2  2004/03/21 22:52:06  mdw
+ * Merge and close elliptic curve branch.
+ *
+ * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+ * Elliptic curves on binary fields work.
+ *
+ */
+
+#ifndef CATACOMB_MPREDUCE_EXP_H
+#define CATACOMB_MPREDUCE_EXP_H
+
+#ifdef __cplusplus
+  extern "C" {
+#endif
+
+/*----- Exponentiation definitions ----------------------------------------*/
+
+#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 = mpreduce_do(mr, t, t);						\
+} while (0)
+
+#define EXP_SQR(a) do {							\
+  mp *t = mp_sqr(spare, a);						\
+  spare = a;								\
+  a = mpreduce_do(mr, t, t);						\
+} while (0)
+
+#define EXP_FIX(x)
+
+#define EXP_SETMUL(d, x, y) do {					\
+  d = mp_mul(MP_NEW, x, y);						\
+  d = mpreduce_do(mr, d, d);						\
+} while (0)
+
+#define EXP_SETSQR(d, x) do {						\
+  d = mp_sqr(MP_NEW, x);						\
+  d = mpreduce_do(mr, d, d);						\
+} while (0)
+
+#include "exp.h"
+
+/*----- That's all, folks -------------------------------------------------*/
+
+#ifdef __cplusplus
+  }
+#endif
+
+#endif
diff --git a/mpreduce.c b/mpreduce.c
new file mode 100644
index 0000000..857549a
--- /dev/null
+++ b/mpreduce.c
@@ -0,0 +1,432 @@
+/* -*-c-*-
+ *
+ * $Id: mpreduce.c,v 1.1 2004/03/27 00:04:46 mdw Exp $
+ *
+ * Efficient reduction modulo nice primes
+ *
+ * (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.
+ */
+
+/*----- Revision history --------------------------------------------------* 
+ *
+ * $Log: mpreduce.c,v $
+ * Revision 1.1  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include <mLib/darray.h>
+#include <mLib/macros.h>
+
+#include "mp.h"
+#include "mpreduce.h"
+#include "mpreduce-exp.h"
+
+/*----- Data structures ---------------------------------------------------*/
+
+DA_DECL(instr_v, mpreduce_instr);
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @mpreduce_create@ --- *
+ *
+ * Arguments:	@gfreduce *r@ = structure to fill in
+ *		@mp *x@ = an integer
+ *
+ * Returns:	---
+ *
+ * Use:		Initializes a context structure for reduction.
+ */
+
+void mpreduce_create(mpreduce *r, mp *p)
+{
+  mpscan sc;
+  enum { Z = 0, Z1 = 2, X = 4, X0 = 6 };
+  unsigned st = Z;
+  instr_v iv = DA_INIT;
+  unsigned long d, i;
+  unsigned op;
+  size_t w, b;
+
+  /* --- Fill in the easy stuff --- */
+
+  assert(MP_ISPOS(p));
+  d = mp_bits(p);
+  r->lim = d/MPW_BITS;
+  r->s = d%MPW_BITS;
+  if (r->s)
+    r->lim++;
+  r->p = mp_copy(p);
+
+  /* --- Stash a new instruction --- */
+
+#define INSTR(op_, argx_, argy_) do {					\
+  DA_ENSURE(&iv, 1);							\
+  DA(&iv)[DA_LEN(&iv)].op = (op_);					\
+  DA(&iv)[DA_LEN(&iv)].argx = (argx_);					\
+  DA(&iv)[DA_LEN(&iv)].argy = (argy_);					\
+  DA_EXTEND(&iv, 1);							\
+} while (0)
+
+  /* --- Main loop --- *
+   *
+   * A simple state machine decomposes @p@ conveniently into positive and
+   * negative powers of 2.  The pure form of the state machine is left below
+   * for reference (and in case I need inspiration for a NAF exponentiator).
+   */
+
+#ifdef DEBUG
+  for (i = 0, mp_scan(&sc, p); mp_step(&sc); i++) {
+    switch (st | mp_bit(&sc)) {
+      case  Z | 1: st = Z1; break;
+      case Z1 | 0: st =  Z; printf("+ %lu\n", i - 1); break;
+      case Z1 | 1: st =  X; printf("- %lu\n", i - 1); break;
+      case  X | 0: st = X0; break;
+      case X0 | 1: st =  X; printf("- %lu\n", i - 1); break;
+      case X0 | 0: st =  Z; printf("+ %lu\n", i - 1); break;
+    }
+  }
+  if (st >= X) printf("+ %lu\n", i);
+#endif
+
+  for (i = 0, mp_scan(&sc, p); i < d  - 1 && mp_step(&sc); i++) {
+    switch (st | mp_bit(&sc)) {
+      case  Z | 1: st = Z1; break;
+      case Z1 | 0: st =  Z; op = MPRI_SUB; goto instr;
+      case Z1 | 1: st =  X; op = MPRI_ADD; goto instr;
+      case  X | 0: st = X0; break;
+      case X0 | 1: st =  X; op = MPRI_ADD; goto instr;
+      case X0 | 0: st =  Z; op = MPRI_SUB; goto instr;
+      instr:
+	w = (d - i)/MPW_BITS + 1;
+	b = (MPW_BITS + i - d - 1)%MPW_BITS;
+	INSTR(op | !!b, w, b);
+    }
+  }
+  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();
+  }
+
+#undef INSTR
+
+  /* --- Wrap up --- */
+
+  r->in = DA_LEN(&iv);
+  if (!r->s) {
+    r->iv = xmalloc(r->in * sizeof(mpreduce_instr));
+    memcpy(r->iv, DA(&iv), r->in * sizeof(mpreduce_instr));
+  } else {
+    r->iv = xmalloc(r->in * 2 * sizeof(mpreduce_instr));
+    for (i = 0; i < r->in; i++) {
+      r->iv[i] = DA(&iv)[i];
+      op = r->iv[i].op & ~1u;
+      w = r->iv[i].argx;
+      b = r->iv[i].argy;
+      b += r->s;
+      if (b >= MPW_BITS) {
+	b -= MPW_BITS;
+	w--;
+      }
+      if (b) op |= 1;
+      r->iv[i + r->in].op = op;
+      r->iv[i + r->in].argx = w;
+      r->iv[i + r->in].argy = b;
+    }
+  }
+  DA_DESTROY(&iv);  
+}
+
+/* --- @mpreduce_destroy@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = structure to free
+ *
+ * Returns:	---
+ *
+ * Use:		Reclaims the resources from a reduction context.
+ */
+
+void mpreduce_destroy(mpreduce *r)
+{
+  mp_drop(r->p);
+  xfree(r->iv);
+}
+
+/* --- @mpreduce_dump@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = structure to dump
+ *		@FILE *fp@ = file to dump on
+ *
+ * Returns:	---
+ *
+ * Use:		Dumps a reduction context.
+ */
+
+void mpreduce_dump(mpreduce *r, FILE *fp)
+{
+  size_t i;
+  static const char *opname[] = { "add", "addshift", "sub", "subshift" };
+
+  fprintf(fp, "mod = "); mp_writefile(r->p, fp, 16);
+  fprintf(fp, "\n  lim = %lu; s = %d\n", (unsigned long)r->lim, r->s);
+  for (i = 0; i < r->in; i++) {
+    assert(r->iv[i].op < N(opname));
+    fprintf(fp, "  %s %lu %lu\n",
+	    opname[r->iv[i].op],
+	    (unsigned long)r->iv[i].argx,
+	    (unsigned long)r->iv[i].argy);
+  }
+  if (r->s) {
+    fprintf(fp, "tail end charlie\n");
+    for (i = r->in; i < 2 * r->in; i++) {
+      assert(r->iv[i].op < N(opname));
+      fprintf(fp, "  %s %lu %lu\n",
+	      opname[r->iv[i].op],
+	      (unsigned long)r->iv[i].argx,
+	      (unsigned long)r->iv[i].argy);
+    }
+  }
+}
+
+/* --- @mpreduce_do@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = reduction context
+ *		@mp *d@ = destination
+ *		@mp *x@ = source
+ *
+ * Returns:	Destination, @x@ reduced modulo the reduction poly.
+ */
+
+static void run(const mpreduce_instr *i, const mpreduce_instr *il,
+		mpw *v, mpw z)
+{
+  for (; i < il; i++) {
+#ifdef DEBUG
+    mp vv;
+    mp_build(&vv, v - i->argx, v + 1);
+    printf("  0x"); mp_writefile(&vv, stdout, 16);
+    printf(" %c (0x%lx << %u) == 0x",
+	   (i->op & ~1u) == MPRI_ADD ? '+' : '-',
+	   (unsigned long)z,
+	   i->argy);
+#endif
+    switch (i->op) {
+      case MPRI_ADD: MPX_UADDN(v - i->argx, v + 1, z); break;
+      case MPRI_ADDLSL: mpx_uaddnlsl(v - i->argx, v + 1, z, i->argy); break;
+      case MPRI_SUB: MPX_USUBN(v - i->argx, v + 1, z); break;
+      case MPRI_SUBLSL: mpx_usubnlsl(v - i->argx, v + 1, z, i->argy); break;
+      default:
+	abort();
+    }
+#ifdef DEBUG
+    mp_build(&vv, v - i->argx, v + 1);
+    mp_writefile(&vv, stdout, 16);
+    printf("\n");
+#endif
+  }
+}
+
+mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
+{
+  mpw *v, *vl;
+  const mpreduce_instr *il;
+  mpw z;
+
+#ifdef DEBUG
+  mp *_r = 0, *_rr = 0;
+#endif
+
+  /* --- If source is negative, divide --- */
+
+  if (MP_ISNEG(x)) {
+    mp_div(0, &d, x, r->p);
+    return (d);
+  }
+
+  /* --- Try to reuse the source's space --- */
+
+  MP_COPY(x);
+  if (d) MP_DROP(d);
+  MP_DEST(x, MP_LEN(x), x->f);
+
+  /* --- Do the reduction --- */
+
+#ifdef DEBUG
+  _r = MP_NEW;
+  mp_div(0, &_r, x, r->p);
+  MP_PRINTX("x", x);
+  _rr = 0;
+#endif
+
+  il = r->iv + r->in;
+  if (MP_LEN(x) >= r->lim) {
+    v = x->v + r->lim;
+    vl = x->vl;
+    while (vl-- > v) {
+      while (*vl) {
+	z = *vl;
+	*vl = 0;
+	run(r->iv, il, vl, z);
+#ifdef DEBUG
+  MP_PRINTX("x", x);
+  mp_div(0, &_rr, x, r->p);
+  assert(MP_EQ(_r, _rr));
+#endif
+      }
+    }
+    if (r->s) {
+      while (*vl >> r->s) {
+	z = *vl >> r->s;
+	*vl &= ((1 << r->s) - 1);
+	run(r->iv + r->in, il + r->in, vl, z);
+#ifdef DEBUG
+  MP_PRINTX("x", x);
+  mp_div(0, &_rr, x, r->p);
+  assert(MP_EQ(_r, _rr));
+#endif
+      }
+    }
+  }
+
+  /* --- Finishing touches --- */
+
+  MP_SHRINK(x);
+  if (MP_CMP(x, >=, r->p))
+    x = mp_sub(x, x, r->p);
+
+  /* --- Done --- */
+
+#ifdef DEBUG
+  assert(MP_EQ(_r, x));
+  mp_drop(_r);
+  mp_drop(_rr);
+#endif
+  return (x);
+}
+
+/* --- @mpreduce_exp@ --- *
+ *
+ * Arguments:	@mpreduce *mr@ = pointer to reduction context
+ *              @mp *d@ = fake destination
+ *              @mp *a@ = base
+ *              @mp *e@ = exponent
+ *
+ * Returns:     Result, %$a^e \bmod m$%.
+ */
+
+mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e)
+{
+  mp *x = MP_ONE;
+  mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
+
+  MP_SHRINK(e);
+  if (!MP_LEN(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);
+  return (x);
+}
+
+/*----- Test rig ----------------------------------------------------------*/
+
+
+#ifdef TEST_RIG
+
+#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
+
+static int vreduce(dstr *v)
+{
+  mp *d = *(mp **)v[0].buf;
+  mp *n = *(mp **)v[1].buf;
+  mp *r = *(mp **)v[2].buf;
+  mp *c;
+  int ok = 1;
+  mpreduce rr;
+
+  mpreduce_create(&rr, d);
+  c = mpreduce_do(&rr, MP_NEW, n);
+  if (!MP_EQ(c, r)) {
+    fprintf(stderr, "\n*** reduction failed\n*** ");
+    mpreduce_dump(&rr, stderr);
+    fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 10);
+    fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 10);
+    fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 10);
+    fprintf(stderr, "\n");
+    ok = 0;
+  }
+  mpreduce_destroy(&rr);
+  mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
+  return (ok);
+}
+
+static int vmodexp(dstr *v)
+{
+  mp *p = *(mp **)v[0].buf;
+  mp *g = *(mp **)v[1].buf;
+  mp *x = *(mp **)v[2].buf;
+  mp *r = *(mp **)v[3].buf;
+  mp *c;
+  int ok = 1;
+  mpreduce rr;
+
+  mpreduce_create(&rr, p);
+  c = mpreduce_exp(&rr, MP_NEW, g, x);
+  if (!MP_EQ(c, r)) {
+    fprintf(stderr, "\n*** modexp failed\n*** ");
+    fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 10);
+    fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 10);
+    fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 10);
+    fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 10);
+    fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 10);
+    fprintf(stderr, "\n");
+    ok = 0;
+  }
+  mpreduce_destroy(&rr);
+  mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
+  return (ok);
+}
+
+static test_chunk defs[] = {
+  { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } },
+  { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+  { 0, 0, { 0 } }
+};
+
+int main(int argc, char *argv[])
+{
+  test_run(argc, argv, defs, SRCDIR"/tests/mpreduce");
+  return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/
diff --git a/mpreduce.h b/mpreduce.h
new file mode 100644
index 0000000..72f27fa
--- /dev/null
+++ b/mpreduce.h
@@ -0,0 +1,136 @@
+/* -*-c-*-
+ *
+ * $Id: mpreduce.h,v 1.1 2004/03/27 00:04:46 mdw Exp $
+ *
+ * Efficient reduction modulo nice primes
+ *
+ * (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.
+ */
+
+/*----- Revision history --------------------------------------------------* 
+ *
+ * $Log: mpreduce.h,v $
+ * Revision 1.1  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
+ */
+
+#ifndef CATACOMB_MPREDUCE_H
+#define CATACOMB_MPREDUCE_H
+
+#ifdef __cplusplus
+  extern "C" {
+#endif
+
+/*----- Header files ------------------------------------------------------*/
+
+/*----- Data structures ---------------------------------------------------*/
+
+typedef struct mpreduce_instr {
+  unsigned op;				/* Instruction opcode */
+  size_t argx, argy;			/* Immediate arguments */
+} mpreduce_instr;
+
+enum {
+  MPRI_ADD,				/* Add @p@ offset by @x@ words */
+  MPRI_ADDLSL,				/* Add @p << y@ offset by @x@ */
+  MPRI_SUB,				/* Sub @p@ offset by @x@ words */
+  MPRI_SUBLSL,				/* Sub @p << y@ offset by @x@ */
+  MPRI_MAX
+};
+
+typedef struct mpreduce {
+  size_t lim;				/* Word containing top bit */
+  unsigned s;				/* Shift for top word */
+  mp *p;				/* Copy of the modulus */
+  size_t in;				/* Number of instruction words */
+  mpreduce_instr *iv;			/* Vector of instructions */
+} mpreduce;
+
+/*----- Functions provided ------------------------------------------------*/
+
+/* --- @mpreduce_create@ --- *
+ *
+ * Arguments:	@gfreduce *r@ = structure to fill in
+ *		@mp *x@ = an integer
+ *
+ * Returns:	---
+ *
+ * Use:		Initializes a context structure for reduction.
+ */
+
+extern void mpreduce_create(mpreduce */*r*/, mp */*p*/);
+
+/* --- @mpreduce_destroy@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = structure to free
+ *
+ * Returns:	---
+ *
+ * Use:		Reclaims the resources from a reduction context.
+ */
+
+extern void mpreduce_destroy(mpreduce */*r*/);
+
+/* --- @mpreduce_dump@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = structure to dump
+ *		@FILE *fp@ = file to dump on
+ *
+ * Returns:	---
+ *
+ * Use:		Dumps a reduction context.
+ */
+
+extern void mpreduce_dump(mpreduce */*r*/, FILE */*fp*/);
+
+/* --- @mpreduce_do@ --- *
+ *
+ * Arguments:	@mpreduce *r@ = reduction context
+ *		@mp *d@ = destination
+ *		@mp *x@ = source
+ *
+ * Returns:	Destination, @x@ reduced modulo the reduction poly.
+ */
+
+extern mp *mpreduce_do(mpreduce */*r*/, mp */*d*/, mp */*x*/);
+
+/* --- @mpreduce_exp@ --- *
+ *
+ * Arguments:	@mpreduce *mr@ = pointer to reduction context
+ *              @mp *d@ = fake destination
+ *              @mp *a@ = base
+ *              @mp *e@ = exponent
+ *
+ * Returns:     Result, %$a^e \bmod m$%.
+ */
+
+extern mp *mpreduce_exp(mpreduce */*mr*/, mp */*d*/, mp */*a*/, mp */*e*/);
+
+/*----- That's all, folks -------------------------------------------------*/
+
+#ifdef __cplusplus
+  }
+#endif
+
+#endif
diff --git a/mpx.c b/mpx.c
index f40d7d9..f1cbbd9 100644
--- a/mpx.c
+++ b/mpx.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.c,v 1.17 2004/03/27 00:04:46 mdw Exp $
  *
  * Low-level multiprecision arithmetic
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.c,v $
+ * Revision 1.17  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.16  2003/05/16 09:09:24  mdw
  * Fix @mp_lsl2c@.  Turns out to be surprisingly tricky.
  *
@@ -877,6 +880,30 @@ void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 
 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
 
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+  mpd x = (mpd)a << o;
+
+  while (x && dv < dvl) {
+    x += *dv;
+    *dv++ = MPW(x);
+    x >>= MPW_BITS;
+  }
+}
+
 /* --- @mpx_usub@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -931,6 +958,33 @@ void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 
 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
 
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+  mpw b = a >> (MPW_BITS - o);
+  a <<= o;
+
+  if (dv < dvl) {
+    mpd x = (mpd)*dv - (mpd)a;
+    *dv++ = MPW(x);
+    if (x >> MPW_BITS)
+      b++;
+    MPX_USUBN(dv, dvl, b);
+  }
+}
+
 /* --- @mpx_umul@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
diff --git a/mpx.h b/mpx.h
index bf73912..13b63cc 100644
--- a/mpx.h
+++ b/mpx.h
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.h,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.h,v 1.17 2004/03/27 00:04:46 mdw Exp $
  *
  * Low level multiprecision arithmetic
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.h,v $
+ * Revision 1.17  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.16  2003/05/16 09:09:24  mdw
  * Fix @mp_lsl2c@.  Turns out to be surprisingly tricky.
  *
@@ -560,6 +563,22 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
 
 extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
 
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+extern void mpx_uaddnlsl(mpw */*dv*/, mpw */*dvl*/,
+			 mpw /*a*/, unsigned /*o*/);
+
 /* --- @mpx_usub@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -607,6 +626,23 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
 
 extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
 
+/* --- @mpx_usubnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d - 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual, so you get two's
+ *		complement.  Which might be what you wanted...  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+extern void mpx_usubnlsl(mpw */*dv*/, mpw */*dvl*/,
+			 mpw /*a*/, unsigned /*o*/);
+
 /* --- @mpx_umul@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
diff --git a/tests/ec b/tests/ec
index 3831c53..1ea76ff 100644
--- a/tests/ec
+++ b/tests/ec
@@ -1,4 +1,4 @@
-# $Id: ec,v 1.1 2004/03/23 15:19:32 mdw Exp $
+# $Id: ec,v 1.2 2004/03/27 00:04:46 mdw Exp $
 #
 # Elliptic curve tests
 
@@ -39,6 +39,42 @@ check {
        0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee"
     0;
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    0;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794810"
+    -1;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee"
+    0;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18801f4ff0afd82ff1411, 
+      0xdccf19d3e76abfa05d529c07575f54c94fa5fc9f3decc246"
+    0;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    0;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794810"
+    -1;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee"
+    0;
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x3f0eba16286a2d57ea0991168d4994637e8343e36,
@@ -91,6 +127,20 @@ find {
     "0x188da80eb03090f67cbf20eb43a18801f4ff0afd82ff1411, 
       0xdccf19d3e76abfa05d529c07575f54c94fa5fc9f3decc246";
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1011 inf;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    0x188da80eb03090f67cbf20eb43a18801f4ff0afd82ff1411
+    "0x188da80eb03090f67cbf20eb43a18801f4ff0afd82ff1411, 
+      0xdccf19d3e76abfa05d529c07575f54c94fa5fc9f3decc246";
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     0x3f0eba16286a2d57ea0991168d4994637e8343e36
@@ -120,6 +170,19 @@ neg {
     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
        0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee";
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee";
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x3f0eba16286a2d57ea0991168d4994637e8343e36,
@@ -148,6 +211,19 @@ dbl {
     "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
        0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x3f0eba16286a2d57ea0991168d4994637e8343e36,
@@ -180,6 +256,23 @@ add {
     "0x76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da,
        0x782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd";
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab"
+    "0x76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da,
+       0x782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab"
+    "0x76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da,
+       0x782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd";
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x3f0eba16286a2d57ea0991168d4994637e8343e36,
@@ -216,6 +309,23 @@ sub {
     "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
        0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da,
+       0x782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+    "0x76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da,
+       0x782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+    "0xdafebf5828783f2ad35534631588a3f629a70fb16982a888,
+       0xdd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab";
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x634000577f86aa315009d6f9b906691f6edd691fe,
@@ -262,6 +372,33 @@ mul {
      6277101735386680763835789423176059013767194773182842284081
      inf;
 
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+     6277101735386680763835789423176059013767194773182842284080
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+        0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     prime: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+     6277101735386680763835789423176059013767194773182842284081
+     inf;
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+     6277101735386680763835789423176059013767194773182842284080
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+        0xf8e6d46a003725879cefee1294db32298c06885ee186b7ee";
+  "niceprime: 6277101735386680763835789423207666416083908700390324961279
+     primeproj: -3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1"
+     "0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+       0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811"
+     6277101735386680763835789423176059013767194773182842284081
+     inf;
+
   "binpoly: 0x800000000000000000000000000000000000000c9
     bin: 1, 0x20a601907b8c953ca1481eb10512f78744a3205fd"
     "0x3f0eba16286a2d57ea0991168d4994637e8343e36,
@@ -288,4 +425,15 @@ mul {
       0x325f41d0ef702dc310254c42d65851a3b91471ac7"
     5846006549323611672814742442876390689256843201587
     inf;
+
+  "binpoly: 0x800000000000000000000000000000000000000c9/bin: 1, 1"
+    "0x2fe13c0537bbc11acaa07d793de4e6d5e5c94eee8,
+      0x289070fb05d38ff58321f2e800536d538ccdaa3d9"
+    5846006549323611672814741753598448348329118574063
+    inf;
+  "binpoly: 0x800000000000000000000000000000000000000c9/binproj: 1, 1"
+    "0x2fe13c0537bbc11acaa07d793de4e6d5e5c94eee8,
+      0x289070fb05d38ff58321f2e800536d538ccdaa3d9"
+    5846006549323611672814741753598448348329118574063
+    inf;
 }
diff --git a/tests/mpreduce b/tests/mpreduce
new file mode 100644
index 0000000..2816b29
--- /dev/null
+++ b/tests/mpreduce
@@ -0,0 +1,21 @@
+# $Id: mpreduce,v 1.1 2004/03/27 00:04:46 mdw Exp $
+#
+# Tests for efficient reduction
+
+reduce {
+  0x72e2c37447f8bca34c4a39b130ea8e5c9a7d8b54564aa88ea773
+  0x367aa8f5ba9ac4e8e2ea198b8af2c3b3081deab392ffc05715783b245a62a6fa
+  0x08e8c03ebf398c63d71d8fd7ca4ece12367a8dde180ca650afb6;
+
+  0xfffffffdffffffffffffffffffffffff
+  0x7fb838a8a0a95046b9d9d9fb4440f7bbc1a7bd3b
+  0xa019c198b9d9d9fb4440f7bc415ff5e4;
+}
+
+modexp {
+  0xfffffffdffffffffffffffffffffffff 2 
+    0xfffffffdfffffffffffffffffffffffe 1;
+  0xfffffffdffffffffffffffffffffffff 2 
+    0xfffffffdfffffffffffffffffffffffd 
+    0x7fffffff000000000000000000000000;
+}
-- 
2.11.0