/* -*-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.
*
* 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 ---------------------------------------------------------*/
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 --- */
/* --- 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@ --- */
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@ --- *
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
* 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);
}
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);
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);
}