X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/5ee4c8931f0cd48cfe38cfc82cdd853f72e60e35..20fa0f6976d598481208c0583d72b2ccef637be9:/mpcrt.c diff --git a/mpcrt.c b/mpcrt.c index 3b1ee22..bf6459f 100644 --- a/mpcrt.c +++ b/mpcrt.c @@ -1,13 +1,13 @@ /* -*-c-*- * - * $Id: mpcrt.c,v 1.1 1999/11/22 20:50:57 mdw Exp $ + * $Id$ * * Chinese Remainder Theorem computations (Gauss's algorithm) * * (c) 1999 Straylight/Edgeware */ -/*----- Licensing notice --------------------------------------------------* +/*----- Licensing notice --------------------------------------------------* * * This file is part of Catacomb. * @@ -15,31 +15,24 @@ * 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: mpcrt.c,v $ - * Revision 1.1 1999/11/22 20:50:57 mdw - * Add support for solving Chinese Remainder Theorem problems. - * - */ - /*----- Header files ------------------------------------------------------*/ #include "mp.h" #include "mpcrt.h" -#include "mpmont.h" +#include "mpmul.h" +#include "mpbarrett.h" /*----- Main code ---------------------------------------------------------*/ @@ -64,7 +57,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 +66,57 @@ 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) { - 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; + if (n != MP_NEW) + n = MP_COPY(n); + else { + mpmul mm; + mpmul_init(&mm); + for (i = 0; i < k; i++) + mpmul_add(&mm, v[i].m); + n = mpmul_done(&mm); + } + + /* --- 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 *g = MP_NEW; + mp_gcd(&g, &v[0].ni, &v[1].ni, v[0].n, v[1].n); + assert(MP_EQ(g, MP_ONE)); + mp_drop(g); + 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@ --- */ @@ -93,20 +124,14 @@ void mpcrt_create(mpcrt *c, mpcrt_mod *v, size_t k, mp *n) if (!v[i].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); - } + v[i].ni = mp_modinv(MP_NEW, v[i].n, v[i].m); + 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 +151,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 +174,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 +211,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)) { + if (!MP_EQ(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 +240,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); + xfree(m); + xfree(r); + assert(mparena_count(MPARENA_GLOBAL) == 0); return (ok); }