/* -*-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.
*
* 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 <mLib/darray.h>
* 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 };
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;
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
/* --- 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 {
r->iv[i + r->in].argy = b;
}
}
- DA_DESTROY(&iv);
+ DA_DESTROY(&iv);
+
+#ifdef DEBUG
+ mpreduce_dump(r, stdout);
+#endif
+ return (0);
}
/* --- @mpreduce_destroy@ --- *
void mpreduce_destroy(mpreduce *r)
{
mp_drop(r->p);
- xfree(r->iv);
+ if (r->iv) xfree(r->iv);
}
/* --- @mpreduce_dump@ --- *
/* --- If source is negative, divide --- */
- if (MP_ISNEG(x)) {
+ if (MP_NEGP(x)) {
mp_div(0, &d, x, r->p);
return (d);
}
*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
}
}
*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
}
}
/* --- @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)
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);
projects / u / mdw / catacomb / blobdiff