X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/f46efa79cd2bb9adc81541f1218965f85a6b2eef..95d9246390251adba7e6e9f0cc70bf0ebe0b2e60:/mpreduce.c diff --git a/mpreduce.c b/mpreduce.c index 857549a..5016f29 100644 --- a/mpreduce.c +++ b/mpreduce.c @@ -1,13 +1,13 @@ /* -*-c-*- * - * $Id: mpreduce.c,v 1.1 2004/03/27 00:04:46 mdw Exp $ + * $Id$ * * Efficient reduction modulo nice primes * * (c) 2004 Straylight/Edgeware */ -/*----- Licensing notice --------------------------------------------------* +/*----- Licensing notice --------------------------------------------------* * * This file is part of Catacomb. * @@ -15,26 +15,18 @@ * 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: mpreduce.c,v $ - * Revision 1.1 2004/03/27 00:04:46 mdw - * Implement efficient reduction for pleasant-looking primes. - * - */ - /*----- Header files ------------------------------------------------------*/ #include @@ -55,12 +47,12 @@ DA_DECL(instr_v, mpreduce_instr); * Arguments: @gfreduce *r@ = structure to fill in * @mp *x@ = an integer * - * Returns: --- + * Returns: Zero if successful; nonzero on failure. * * Use: Initializes a context structure for reduction. */ -void mpreduce_create(mpreduce *r, mp *p) +int mpreduce_create(mpreduce *r, mp *p) { mpscan sc; enum { Z = 0, Z1 = 2, X = 4, X0 = 6 }; @@ -68,11 +60,12 @@ void mpreduce_create(mpreduce *r, mp *p) instr_v iv = DA_INIT; unsigned long d, i; unsigned op; - size_t w, b; + size_t w, b, bb; /* --- Fill in the easy stuff --- */ - assert(MP_ISPOS(p)); + if (!MP_POSP(p)) + return (-1); d = mp_bits(p); r->lim = d/MPW_BITS; r->s = d%MPW_BITS; @@ -101,34 +94,36 @@ void mpreduce_create(mpreduce *r, mp *p) for (i = 0, mp_scan(&sc, p); mp_step(&sc); i++) { switch (st | mp_bit(&sc)) { case Z | 1: st = Z1; break; - case Z1 | 0: st = Z; printf("+ %lu\n", i - 1); break; - case Z1 | 1: st = X; printf("- %lu\n", i - 1); break; + case Z1 | 0: st = Z; printf("+ %lu\n", i - 1); break; + case Z1 | 1: st = X; printf("- %lu\n", i - 1); break; case X | 0: st = X0; break; - case X0 | 1: st = X; printf("- %lu\n", i - 1); break; - case X0 | 0: st = Z; printf("+ %lu\n", i - 1); break; + case X0 | 1: st = X; printf("- %lu\n", i - 1); break; + case X0 | 0: st = Z; printf("+ %lu\n", i - 1); break; } } - if (st >= X) printf("+ %lu\n", i); + if (st >= X) printf("+ %lu\n", i - 1); + st = Z; #endif - for (i = 0, mp_scan(&sc, p); i < d - 1 && mp_step(&sc); i++) { + bb = MPW_BITS - (d + 1)%MPW_BITS; + for (i = 0, mp_scan(&sc, p); i < d && mp_step(&sc); i++) { switch (st | mp_bit(&sc)) { case Z | 1: st = Z1; break; - case Z1 | 0: st = Z; op = MPRI_SUB; goto instr; - case Z1 | 1: st = X; op = MPRI_ADD; goto instr; + case Z1 | 0: st = Z; op = MPRI_SUB; goto instr; + case Z1 | 1: st = X; op = MPRI_ADD; goto instr; case X | 0: st = X0; break; - case X0 | 1: st = X; op = MPRI_ADD; goto instr; - case X0 | 0: st = Z; op = MPRI_SUB; goto instr; + case X0 | 1: st = X; op = MPRI_ADD; goto instr; + case X0 | 0: st = Z; op = MPRI_SUB; goto instr; instr: w = (d - i)/MPW_BITS + 1; - b = (MPW_BITS + i - d - 1)%MPW_BITS; + b = (bb + i)%MPW_BITS; INSTR(op | !!b, w, b); } } - if ((DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) { - fprintf(stderr, - "mpreduce can't handle moduli of the form 2^m + x\n"); - abort(); + if (DA_LEN(&iv) && (DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) { + mp_drop(r->p); + DA_DESTROY(&iv); + return (-1); } #undef INSTR @@ -136,7 +131,9 @@ void mpreduce_create(mpreduce *r, mp *p) /* --- Wrap up --- */ r->in = DA_LEN(&iv); - if (!r->s) { + if (!r->in) + r->iv = 0; + else if (!r->s) { r->iv = xmalloc(r->in * sizeof(mpreduce_instr)); memcpy(r->iv, DA(&iv), r->in * sizeof(mpreduce_instr)); } else { @@ -157,7 +154,12 @@ void mpreduce_create(mpreduce *r, mp *p) r->iv[i + r->in].argy = b; } } - DA_DESTROY(&iv); + DA_DESTROY(&iv); + +#ifdef DEBUG + mpreduce_dump(r, stdout); +#endif + return (0); } /* --- @mpreduce_destroy@ --- * @@ -172,7 +174,7 @@ void mpreduce_create(mpreduce *r, mp *p) void mpreduce_destroy(mpreduce *r) { mp_drop(r->p); - xfree(r->iv); + if (r->iv) xfree(r->iv); } /* --- @mpreduce_dump@ --- * @@ -261,7 +263,7 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x) /* --- If source is negative, divide --- */ - if (MP_ISNEG(x)) { + if (MP_NEGP(x)) { mp_div(0, &d, x, r->p); return (d); } @@ -291,9 +293,9 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x) *vl = 0; run(r->iv, il, vl, z); #ifdef DEBUG - MP_PRINTX("x", x); - mp_div(0, &_rr, x, r->p); - assert(MP_EQ(_r, _rr)); + MP_PRINTX("x", x); + mp_div(0, &_rr, x, r->p); + assert(MP_EQ(_r, _rr)); #endif } } @@ -303,9 +305,9 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x) *vl &= ((1 << r->s) - 1); run(r->iv + r->in, il + r->in, vl, z); #ifdef DEBUG - MP_PRINTX("x", x); - mp_div(0, &_rr, x, r->p); - assert(MP_EQ(_r, _rr)); + MP_PRINTX("x", x); + mp_div(0, &_rr, x, r->p); + assert(MP_EQ(_r, _rr)); #endif } } @@ -330,11 +332,11 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x) /* --- @mpreduce_exp@ --- * * * Arguments: @mpreduce *mr@ = pointer to reduction context - * @mp *d@ = fake destination - * @mp *a@ = base - * @mp *e@ = exponent + * @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent * - * Returns: Result, %$a^e \bmod m$%. + * Returns: Result, %$a^e \bmod m$%. */ mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e) @@ -343,12 +345,18 @@ mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e) mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; MP_SHRINK(e); - if (!MP_LEN(e)) + MP_COPY(a); + if (MP_ZEROP(e)) ; - else if (MP_LEN(e) < EXP_THRESH) - EXP_SIMPLE(x, a, e); - else - EXP_WINDOW(x, a, e); + else { + if (MP_NEGP(e)) + a = mp_modinv(a, a, mr->p); + if (MP_LEN(e) < EXP_THRESH) + EXP_SIMPLE(x, a, e); + else + EXP_WINDOW(x, a, e); + } + mp_drop(a); mp_drop(d); mp_drop(spare); return (x);