Interface changes for suggested destinations. Use Barrett reduction.

diff --git a/mpcrt.c b/mpcrt.c
index 3b1ee22..51d9b68 100644 (file)
--- a/mpcrt.c
+++ b/mpcrt.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpcrt.c,v 1.1 1999/11/22 20:50:57 mdw Exp $
+ * $Id: mpcrt.c,v 1.2 1999/12/10 23:22:32 mdw Exp $
  *
  * Chinese Remainder Theorem computations (Gauss's algorithm)
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpcrt.c,v $
+ * Revision 1.2  1999/12/10 23:22:32  mdw
+ * Interface changes for suggested destinations.  Use Barrett reduction.
+ *
  * Revision 1.1  1999/11/22 20:50:57  mdw
  * Add support for solving Chinese Remainder Theorem problems.
  *
@@ -39,7 +42,7 @@
 
 #include "mp.h"
 #include "mpcrt.h"
-#include "mpmont.h"
+#include "mpbarrett.h"
 
 /*----- Main code ---------------------------------------------------------*/
 
@@ -64,7 +67,6 @@
 
 void mpcrt_create(mpcrt *c, mpcrt_mod *v, size_t k, mp *n)
 {
-  mp *x = MP_NEW, *y = MP_NEW;
   size_t i;
 
   /* --- Simple initialization things --- */
@@ -74,18 +76,52 @@ void mpcrt_create(mpcrt *c, mpcrt_mod *v, size_t k, mp *n)
 
   /* --- Work out @n@ if I don't have it already --- */
 
-  if (n == MP_NEW) {
+  if (n != MP_NEW)
+    n = MP_COPY(n);
+  else {
     n = MP_COPY(v[0].m);
-    for (i = 1; i < k; i++) {
-      mp *d = mp_mul(x, n, v[i].m);
-      x = n;
-      n = d;
+    for (i = 1; i < k; i++)
+      n = mp_mul(n, n, v[i].m);
+  }
+
+  /* --- A quick hack if %$k = 2$% --- */
+
+  if (k == 2) {
+
+    /* --- The %$n / n_i$% values are trivial in this case --- */
+
+    if (!v[0].n)
+      v[0].n = MP_COPY(v[1].m);
+    if (!v[1].n)
+      v[1].n = MP_COPY(v[0].m);
+
+    /* --- Now sort out the inverses --- *
+     *
+     * @mp_gcd@ will ensure that the first argument is negative.
+     */
+
+    if (!v[0].ni && !v[1].ni) {
+      mp_gcd(0, &v[0].ni, &v[1].ni, v[0].n, v[1].n);
+      v[0].ni = mp_add(v[0].ni, v[0].ni, v[1].n);
+    } else {
+      int i, j;
+      mp *x;
+
+      if (!v[0].ni)
+       i = 0, j = 1;
+      else
+       i = 1, j = 0;
+      
+      x = mp_mul(MP_NEW, v[j].n, v[j].ni);
+      x = mp_sub(x, x, MP_ONE);
+      mp_div(&x, 0, x, v[i].n);
+      v[i].ni = x;
     }
   }
 
-  /* --- Set up the Montgomery context --- */
+  /* --- Set up the Barrett context --- */
 
-  mpmont_create(&c->mm, n);
+  mpbarrett_create(&c->mb, n);
 
   /* --- Walk through filling in @n@, @ni@ and @nnir@ --- */
 
@@ -94,19 +130,13 @@ void mpcrt_create(mpcrt *c, mpcrt_mod *v, size_t k, mp *n)
       mp_div(&v[i].n, 0, n, v[i].m);
     if (!v[i].ni)
       mp_gcd(0, &v[i].ni, 0, v[i].n, v[i].m);
-    if (!v[i].nnir) {
-      x = mpmont_mul(&c->mm, x, v[i].n, c->mm.r2);
-      y = mpmont_mul(&c->mm, y, v[i].ni, c->mm.r2);
-      v[i].nnir = mpmont_mul(&c->mm, MP_NEW, x, y);
-    }
+    if (!v[i].nni)
+      v[i].nni = mp_mul(MP_NEW, v[i].n, v[i].ni);
   }
 
   /* --- Done --- */
 
-  if (x)
-    mp_drop(x);
-  if (y)
-    mp_drop(y);
+  mp_drop(n);
 }
 
 /* --- @mpcrt_destroy@ --- *
@@ -126,14 +156,15 @@ void mpcrt_destroy(mpcrt *c)
     if (c->v[i].m) mp_drop(c->v[i].m);
     if (c->v[i].n) mp_drop(c->v[i].n);
     if (c->v[i].ni) mp_drop(c->v[i].ni);
-    if (c->v[i].nnir) mp_drop(c->v[i].nnir);
+    if (c->v[i].nni) mp_drop(c->v[i].nni);
   }
-  mpmont_destroy(&c->mm);
+  mpbarrett_destroy(&c->mb);
 }
 
 /* --- @mpcrt_solve@ --- *
  *
  * Arguments:  @mpcrt *c@ = pointer to CRT context
+ *             @mp *d@ = fake destination
  *             @mp **v@ = array of residues
  *
  * Returns:    The unique solution modulo the product of the individual
@@ -148,20 +179,22 @@ void mpcrt_destroy(mpcrt *c)
  *             exponentiation.
  */
 
-mp *mpcrt_solve(mpcrt *c, mp **v)
+mp *mpcrt_solve(mpcrt *c, mp *d, mp **v)
 {
   mp *a = MP_ZERO;
   mp *x = MP_NEW;
   size_t i;
 
   for (i = 0; i < c->k; i++) {
-    x = mpmont_mul(&c->mm, x, c->v[i].nnir, v[i]);
+    x = mp_mul(x, c->v[i].nni, v[i]);
+    x = mpbarrett_reduce(&c->mb, x, x);
     a = mp_add(a, a, x);
   }
   if (x)
-    mp_drop(x);
-  if (MP_CMP(a, >=, c->mm.m))
-    mp_div(0, &a, a, c->mm.m);
+    MP_DROP(x);
+  a = mpbarrett_reduce(&c->mb, a, a);
+  if (d != MP_NEW)
+    MP_DROP(d);
   return (a);
 }
 
@@ -183,17 +216,17 @@ static int verify(size_t n, dstr *v)
     m[i].m = *(mp **)v[2 * i + 1].buf;
     m[i].n = 0;
     m[i].ni = 0;
-    m[i].nnir = 0;
+    m[i].nni = 0;
   }
   a = *(mp **)v[2 * n].buf;
 
   mpcrt_create(&c, m, n, 0);
-  b = mpcrt_solve(&c, r);
+  b = mpcrt_solve(&c, MP_NEW, r);
 
   if (MP_CMP(a, !=, b)) {
     fputs("\n*** failed\n", stderr);
     fputs("n = ", stderr);
-    mp_writefile(c.mm.m, stderr, 10);
+    mp_writefile(c.mb.m, stderr, 10);
     for (i = 0; i < n; i++) {
       fprintf(stderr, "\nr[%u] = ", i);
       mp_writefile(r[i], stderr, 10);
@@ -212,11 +245,14 @@ static int verify(size_t n, dstr *v)
     ok = 0;
   }
 
+  for (i = 0; i < n; i++)
+    mp_drop(r[i]);
   mp_drop(a);
   mp_drop(b);
   mpcrt_destroy(&c);
   free(m);
   free(r);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 

diff --git a/mpcrt.h b/mpcrt.h
index 5474a92..a5dd71c 100644 (file)
--- a/mpcrt.h
+++ b/mpcrt.h
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpcrt.h,v 1.1 1999/11/22 20:50:57 mdw Exp $
+ * $Id: mpcrt.h,v 1.2 1999/12/10 23:22:32 mdw Exp $
  *
  * Chinese Remainder Theorem computations (Gauss's algorithm)
  *
@@ -30,13 +30,16 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpcrt.h,v $
+ * Revision 1.2  1999/12/10 23:22:32  mdw
+ * Interface changes for suggested destinations.  Use Barrett reduction.
+ *
  * Revision 1.1  1999/11/22 20:50:57  mdw
  * Add support for solving Chinese Remainder Theorem problems.
  *
  */
 
-#ifndef MPCRT_H
-#define MPCRT_H
+#ifndef CATACOMB_MPCRT_H
+#define CATACOMB_MPCRT_H
 
 #ifdef __cplusplus
   extern "C" {
@@ -46,12 +49,12 @@
 
 #include <stddef.h>
 
-#ifndef MP_H
+#ifndef CATACOMB_MP_H
 #  include "mp.h"
 #endif
 
-#ifndef MPMONT_H
-#  include "mpmont.h"
+#ifndef CATACOMB_MPBARRETT_H
+#  include "mpbarrett.h"
 #endif
 
 /*----- Data structures ---------------------------------------------------*/
@@ -60,12 +63,12 @@ typedef struct mpcrt_mod {
   mp *m;                               /* %$n_i$% -- the modulus */
   mp *n;                               /* %$N_i = n / n_i$% */
   mp *ni;                              /* %$M_i = N_i^{-1} \bmod n_i$% */
-  mp *nnir;                            /* %$N_i M_i R \bmod m$% */
+  mp *nni;                             /* %$N_i M_i \bmod m$% */
 } mpcrt_mod;
 
 typedef struct mpcrt {
   size_t k;                            /* Number of distinct moduli */
-  mpmont mm;                           /* Montgomery context for product */
+  mpbarrett mb;                                /* Barrett context for product */
   mpcrt_mod *v;                                /* Vector of information for each */
 } mpcrt;
 
@@ -107,6 +110,7 @@ extern void mpcrt_destroy(mpcrt */*c*/);
 /* --- @mpcrt_solve@ --- *
  *
  * Arguments:  @mpcrt *c@ = pointer to CRT context
+ *             @mp *d@ = fake destination
  *             @mp **v@ = array of residues
  *
  * Returns:    The unique solution modulo the product of the individual
@@ -121,7 +125,7 @@ extern void mpcrt_destroy(mpcrt */*c*/);
  *             exponentiation.
  */
 
-extern mp *mpcrt_solve(mpcrt */*c*/, mp **/*v*/);
+extern mp *mpcrt_solve(mpcrt */*c*/, mp */*d*/, mp **/*v*/);
 
 /*----- That's all, folks -------------------------------------------------*/