X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/mpmont.c diff --git a/mpmont.c b/mpmont.c deleted file mode 100644 index 65b3fcc..0000000 --- a/mpmont.c +++ /dev/null @@ -1,429 +0,0 @@ -/* -*-c-*- - * - * $Id$ - * - * Montgomery reduction - * - * (c) 1999 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. - */ - -/*----- Header files ------------------------------------------------------*/ - -#include "mp.h" -#include "mpmont.h" - -/*----- Tweakables --------------------------------------------------------*/ - -/* --- @MPMONT_DISABLE@ --- * - * - * Replace all the clever Montgomery reduction with good old-fashioned long - * division. - */ - -/* #define MPMONT_DISABLE */ - -/*----- Reduction and multiplication --------------------------------------*/ - -/* --- @mpmont_create@ --- * - * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context - * @mp *m@ = modulus to use - * - * Returns: Zero on success, nonzero on error. - * - * Use: Initializes a Montgomery reduction context ready for use. - * The argument @m@ must be a positive odd integer. - */ - -#ifdef MPMONT_DISABLE - -int mpmont_create(mpmont *mm, mp *m) -{ - mp_shrink(m); - mm->m = MP_COPY(m); - mm->r = MP_ONE; - mm->r2 = MP_ONE; - mm->mi = MP_ONE; - return (0); -} - -#else - -int mpmont_create(mpmont *mm, mp *m) -{ - size_t n = MP_LEN(m); - mp *r2 = mp_new(2 * n + 1, 0); - mp r; - - /* --- Take a copy of the modulus --- */ - - if (!MP_POSP(m) || !MP_ODDP(m)) - return (-1); - mm->m = MP_COPY(m); - - /* --- Determine %$R^2$% --- */ - - mm->n = n; - MPX_ZERO(r2->v, r2->vl - 1); - r2->vl[-1] = 1; - - /* --- Find the magic value @mi@ --- */ - - mp_build(&r, r2->v + n, r2->vl); - mm->mi = mp_modinv(MP_NEW, m, &r); - mm->mi = mp_sub(mm->mi, &r, mm->mi); - - /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */ - - mm->r2 = MP_NEW; - mp_div(0, &mm->r2, r2, m); - mm->r = mpmont_reduce(mm, MP_NEW, mm->r2); - MP_DROP(r2); - return (0); -} - -#endif - -/* --- @mpmont_destroy@ --- * - * - * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context - * - * Returns: --- - * - * Use: Disposes of a context when it's no longer of any use to - * anyone. - */ - -void mpmont_destroy(mpmont *mm) -{ - MP_DROP(mm->m); - MP_DROP(mm->r); - MP_DROP(mm->r2); - MP_DROP(mm->mi); -} - -/* --- @mpmont_reduce@ --- * - * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context - * @mp *d@ = destination - * @mp *a@ = source, assumed positive - * - * Returns: Result, %$a R^{-1} \bmod m$%. - */ - -#ifdef MPMONT_DISABLE - -mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) -{ - mp_div(0, &d, a, mm->m); - return (d); -} - -#else - -mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) -{ - size_t n = mm->n; - - /* --- Check for serious Karatsuba reduction --- */ - - if (n > MPK_THRESH * 3) { - mp al; - mpw *vl; - mp *u; - - if (MP_LEN(a) >= n) - vl = a->v + n; - else - vl = a->vl; - mp_build(&al, a->v, vl); - u = mp_mul(MP_NEW, &al, mm->mi); - if (MP_LEN(u) > n) - u->vl = u->v + n; - u = mp_mul(u, u, mm->m); - d = mp_add(d, a, u); - mp_drop(u); - } - - /* --- Otherwise do it the hard way --- */ - - else { - mpw *dv, *dvl; - mpw *mv, *mvl; - mpw mi; - size_t k = n; - - /* --- Initial conditioning of the arguments --- */ - - a = MP_COPY(a); - if (d) - MP_DROP(d); - d = a; - MP_DEST(d, 2 * n + 1, a->f); - - dv = d->v; dvl = d->vl; - mv = mm->m->v; mvl = mm->m->vl; - - /* --- Let's go to work --- */ - - mi = mm->mi->v[0]; - while (k--) { - mpw u = MPW(*dv * mi); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - } - } - - /* --- Wrap everything up --- */ - - memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n)); - d->vl -= n; - if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl)) - mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl); - if (d->f & MP_NEG) { - mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl); - d->f &= ~MP_NEG; - } - MP_SHRINK(d); - return (d); -} - -#endif - -/* --- @mpmont_mul@ --- * - * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context - * @mp *d@ = destination - * @mp *a, *b@ = sources, assumed positive - * - * Returns: Result, %$a b R^{-1} \bmod m$%. - */ - -#ifdef MPMONT_DISABLE - -mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) -{ - d = mp_mul(d, a, b); - mp_div(0, &d, d, mm->m); - return (d); -} - -#else - -mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) -{ - if (mm->n > MPK_THRESH * 3) { - d = mp_mul(d, a, b); - d = mpmont_reduce(mm, d, d); - } else { - mpw *dv, *dvl; - mpw *av, *avl; - mpw *bv, *bvl; - mpw *mv, *mvl; - mpw y; - size_t n, i; - mpw mi; - - /* --- Initial conditioning of the arguments --- */ - - if (MP_LEN(a) > MP_LEN(b)) { - mp *t = a; a = b; b = t; - } - n = MP_LEN(mm->m); - - a = MP_COPY(a); - b = MP_COPY(b); - MP_DEST(d, 2 * n + 1, a->f | b->f | MP_UNDEF); - dv = d->v; dvl = d->vl; - MPX_ZERO(dv, dvl); - av = a->v; avl = a->vl; - bv = b->v; bvl = b->vl; - mv = mm->m->v; mvl = mm->m->vl; - y = *bv; - - /* --- Montgomery multiplication phase --- */ - - i = 0; - mi = mm->mi->v[0]; - while (i < n && av < avl) { - mpw x = *av++; - mpw u = MPW((*dv + x * y) * mi); - MPX_UMLAN(dv, dvl, bv, bvl, x); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - i++; - } - - /* --- Simpler Montgomery reduction phase --- */ - - while (i < n) { - mpw u = MPW(*dv * mi); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - i++; - } - - /* --- Done --- */ - - memmove(d->v, dv, MPWS(dvl - dv)); - d->vl -= dv - d->v; - if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl)) - mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl); - if ((a->f ^ b->f) & MP_NEG) - mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl); - MP_SHRINK(d); - d->f = (a->f | b->f) & MP_BURN; - MP_DROP(a); - MP_DROP(b); - } - - return (d); -} - -#endif - -/*----- Test rig ----------------------------------------------------------*/ - -#ifdef TEST_RIG - -static int tcreate(dstr *v) -{ - mp *m = *(mp **)v[0].buf; - mp *mi = *(mp **)v[1].buf; - mp *r = *(mp **)v[2].buf; - mp *r2 = *(mp **)v[3].buf; - - mpmont mm; - int ok = 1; - - mpmont_create(&mm, m); - - if (mm.mi->v[0] != mi->v[0]) { - fprintf(stderr, "\n*** bad mi: found %lu, expected %lu", - (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]); - fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - if (!MP_EQ(mm.r, r)) { - fputs("\n*** bad r", stderr); - fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); - fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10); - fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - if (!MP_EQ(mm.r2, r2)) { - fputs("\n*** bad r2", stderr); - fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); - fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10); - fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - MP_DROP(m); - MP_DROP(mi); - MP_DROP(r); - MP_DROP(r2); - mpmont_destroy(&mm); - assert(mparena_count(MPARENA_GLOBAL) == 0); - return (ok); -} - -static int tmul(dstr *v) -{ - mp *m = *(mp **)v[0].buf; - mp *a = *(mp **)v[1].buf; - mp *b = *(mp **)v[2].buf; - mp *r = *(mp **)v[3].buf; - int ok = 1; - - mpmont mm; - mpmont_create(&mm, m); - - { - mp *qr = mp_mul(MP_NEW, a, b); - mp_div(0, &qr, qr, m); - - if (!MP_EQ(qr, r)) { - fputs("\n*** classical modmul failed", stderr); - fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); - fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); - fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); - fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); - fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - mp_drop(qr); - } - - { - mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2); - mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2); - mp *mr = mpmont_mul(&mm, MP_NEW, ar, br); - mr = mpmont_reduce(&mm, mr, mr); - if (!MP_EQ(mr, r)) { - fputs("\n*** montgomery modmul failed", stderr); - fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); - fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); - fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); - fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); - fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - MP_DROP(ar); MP_DROP(br); - mp_drop(mr); - } - - - MP_DROP(m); - MP_DROP(a); - MP_DROP(b); - MP_DROP(r); - mpmont_destroy(&mm); - assert(mparena_count(MPARENA_GLOBAL) == 0); - return ok; -} - -static test_chunk tests[] = { - { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, - { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, - { 0, 0, { 0 } }, -}; - -int main(int argc, char *argv[]) -{ - sub_init(); - test_run(argc, argv, tests, SRCDIR "/tests/mpmont"); - return (0); -} - -#endif - -/*----- That's all, folks -------------------------------------------------*/