From 5bda60bd9209028aa3db06bccd7c708a1d3a9c5e Mon Sep 17 00:00:00 2001 From: mdw Date: Fri, 10 Dec 1999 23:22:32 +0000 Subject: [PATCH] Interface changes for suggested destinations. Use Barrett reduction. --- mpcrt.c | 94 +++++++++++++++++++++++++++++++++++++++++++++-------------------- mpcrt.h | 22 ++++++++------- 2 files changed, 78 insertions(+), 38 deletions(-) diff --git a/mpcrt.c b/mpcrt.c index 3b1ee22..51d9b68 100644 --- 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 --- 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 -#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 -------------------------------------------------*/ -- 2.11.0