From ceb3f0c0a3b7bb3fa3250d31b04c382894095e52 Mon Sep 17 00:00:00 2001 From: mdw Date: Sun, 21 Mar 2004 22:39:46 +0000 Subject: [PATCH] Elliptic curves on binary fields work. --- .gdbinit | 14 +- Makefile.m4 | 16 +- calc/ec2.cal | 180 ++++++++++++++++ calc/gfx.cal | 35 +++- ec-bin.c | 428 +++++++++++++++++++++++++++++++++++++ ec-prime.c | 52 ++--- ec.h | 13 +- f-binpoly.c | 142 +++++++++++++ f-prime.c | 30 +-- field.h | 29 ++- gf-arith.c | 263 +++++++++++++++++++++++ gf-gcd.c | 259 +++++++++++++++++++++++ gf.h | 124 +++++++++++ gfreduce-exp.h | 84 ++++++++ gfreduce.c | 652 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ gfreduce.h | 186 ++++++++++++++++ gfx-sqr.c | 9 +- gfx.h | 18 +- mp-gcd.c | 15 +- tests/gf | 70 +++++++ tests/gfreduce | 61 ++++++ tests/gfx | 7 +- 22 files changed, 2601 insertions(+), 86 deletions(-) create mode 100644 calc/ec2.cal create mode 100644 ec-bin.c create mode 100644 f-binpoly.c create mode 100644 gf-arith.c create mode 100644 gf-gcd.c create mode 100644 gf.h create mode 100644 gfreduce-exp.h create mode 100644 gfreduce.c create mode 100644 gfreduce.h create mode 100644 tests/gf create mode 100644 tests/gfreduce diff --git a/.gdbinit b/.gdbinit index daa7ee3..c4147c3 100644 --- a/.gdbinit +++ b/.gdbinit @@ -9,22 +9,22 @@ define mp-print end define mp-printr - call (void)fputs("$arg0 = ", stdout) - if $arg0 == 0 + call (void)fputs("$arg1 = ", stdout) + if $arg1 == 0 call (void)fputs("(null)", stdout) else - if $arg1 == 16 + if $arg0 == 16 call (void)fputs("0x", stdout) else - if $arg1 == 8 + if $arg0 == 8 call (void)fputs("0", stdout) else - if $arg1 != 10 - call (void)fputs("$arg1", stdout) + if $arg0 != 10 + call (void)fputs("$arg0:", stdout) end end end - call (void)mp_writefile($arg0, stdout, $arg1) + call (void)mp_writefile($arg1, stdout, $arg0) end call (void)putchar('\n') end diff --git a/Makefile.m4 b/Makefile.m4 index 8185765..4a95044 100644 --- a/Makefile.m4 +++ b/Makefile.m4 @@ -1,6 +1,6 @@ ## -*-makefile-*- ## -## $Id: Makefile.m4,v 1.60.2.1 2003/06/10 13:43:53 mdw Exp $ +## $Id: Makefile.m4,v 1.60.2.2 2004/03/21 22:39:46 mdw Exp $ ## ## Makefile for Catacomb ## @@ -29,6 +29,9 @@ ##----- Revision history ---------------------------------------------------- ## ## $Log: Makefile.m4,v $ +## Revision 1.60.2.2 2004/03/21 22:39:46 mdw +## Elliptic curves on binary fields work. +## ## Revision 1.60.2.1 2003/06/10 13:43:53 mdw ## Simple (non-projective) curves over prime fields now seem to work. ## @@ -316,7 +319,7 @@ pkginclude_HEADERS = \ mpx.h bitops.h mpw.h mpscan.h mparena.h mp.h mptext.h mpint.h \ exp.h mpbarrett.h mpbarrett-exp.h mpmont.h mpmont-exp.h \ mpcrt.h mprand.h mpmul.h \ - gfx.h \ + gfx.h gf.h \ primetab.h pfilt.h rabin.h \ pgen.h prim.h strongprime.h limlee.h keycheck.h \ bbs.h rsa.h dh.h dsarand.h dsa.h \ @@ -342,10 +345,10 @@ define(`MP_SOURCES', GF_SOURCES PGEN_SOURCES EC_SOURCES') define(`GF_SOURCES', - `gfx.c gfx-kmul.c gfx-sqr.c') + `gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-gcd.c gfreduce.c') define(`EC_SOURCES', - `field.c f-prime.c ec.c ec-prime.c') + `field.c f-prime.c f-binpoly.c ec.c ec-prime.c ec-bin.c') define(`PGEN_SOURCES', `pfilt.c rabin.c \ @@ -511,8 +514,13 @@ CTESTRIG(mpmont-mexp) CTESTRIG(mpcrt) CTESTRIG(mpmul) CTESTRIG(gfx) +CTESTRIG(gfx-sqr) CTESTRIG(gfx-kmul) +CTESTRIG(gf-arith) +CTESTRIG(gf-gcd) +CTESTRIG(gfreduce) CTESTRIG(ec-prime) +CTESTRIG(ec-bin) CTESTRIG(pgen) CTESTRIG(dsa-gen) CTESTRIG(dsa-sign) diff --git a/calc/ec2.cal b/calc/ec2.cal new file mode 100644 index 0000000..0d7ceb9 --- /dev/null +++ b/calc/ec2.cal @@ -0,0 +1,180 @@ +/* -*-apcalc-*- + * + * $Id: ec2.cal,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Testbed for elliptic curve arithmetic over binary fields + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: ec2.cal,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + * Revision 1.1.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * + * Revision 1.1.4.1 2003/06/10 13:43:53 mdw + * Simple (non-projective) curves over prime fields now seem to work. + * + * Revision 1.1 2000/10/08 16:01:37 mdw + * Prototypes of various bits of code. + * + */ + +/*----- Object types ------------------------------------------------------*/ + +obj ec2_curve { a, b, p }; +obj ec2_pt { x, y, e }; +obj ecpp_pt { x, y, z, e }; + +/*----- Main code ---------------------------------------------------------*/ + +define ec2_curve(a, b, p) +{ + local obj ec2_curve e; + e.a = a; + e.b = b; + e.p = p; + return (e); +} + +define ec2_pt(x, y, e) +{ + local obj ec2_pt p; + p.x = x % e.p; + p.y = y % e.p; + p.e = e; + return (p); +} + +define ec2_pt_print(a) +{ + print "(" : a.x : ", " : a.y : ")" :; +} + +define ec2_pt_add(a, b) +{ + local e, alpha; + local obj ec2_pt d; + + print "> ecadd: ", a, b; + if (a == 0) + d = b; + else if (b == 0) + d = a; + else if (!istype(a, b)) + quit "bad type arguments to ec2_pt_add"; + else if (a.e != b.e) + quit "points from different curves in ec2_pt_add"; + else { + e = a.e; + if (a.x != b.x) { + alpha = ((a.y + b.y) * gf_inv(a.x + b.x, e.p)) % e.p; + d.x = (e.a + alpha^2 + alpha + a.x + b.x) % e.p; + } else if (a.y != b.y || a.x == gf(0)) + return 0; + else { + alpha = a.x + a.y * gf_inv(a.x, e.p) % e.p; + d.x = (e.a + alpha^2 + alpha) % e.p; + } + d.y = ((a.x + d.x) * alpha + d.x + a.y) % e.p; + d.e = e; + } + + print "< ecadd: ", d; + return (d); +} + +define ec2_pt_dbl(a) +{ + local e, alpha; + local obj ec2_pt d; + print "> ecdbl: ", a; + if (istype(a, 1)) + return (0); + e = a.e; + alpha = a.x + a.y * gf_inv(a.x, e.p) % e.p; + d.x = (e.a + alpha^2 + alpha) % e.p; + d.y = ((a.x + d.x) * alpha + d.x + a.y) % e.p; + d.e = e; + print "< ecdbl: ", d; + return (d); +} + +define ec2_pt_sub(a, b) { return ec2_pt_add(a, ec2_pt_neg(b)); } + +define ec2_pt_neg(a) +{ + local obj ec2_pt d; + d.x = a.x; + d.y = a.x + a.y; + d.e = a.e; + return (d); +} + +define ec2_pt_check(a) +{ + local e; + + e = a.e; + if ((a.y^2 + a.x * a.y) % e.p != (a.x^3 + e.a * a.x^2 + e.b) % e.p) + quit "bad curve point"; +} + +define ec2_pt_mul(a, b) +{ + local p, n; + local d; + + if (istype(a, 1)) { + n = a; + p = b; + } else if (istype(b, 1)) { + n = b; + p = a; + } else + return (newerror("bad arguments to ec2_pt_mul")); + + d = 0; + while (n) { + if (n & 1) + d += p; + n >>= 1; + p = ec2_pt_dbl(p); + } + return (d); +} + +/*----- FIPS186-2 standard curves -----------------------------------------*/ + +b163 = ec2_curve(gf(1),gf(0x20a601907b8c953ca1481eb10512f78744a3205fd), + gf(0x800000000000000000000000000000000000000c9)); +b163_r = 5846006549323611672814742442876390689256843201587; +b163_g = ec2_pt(0x3f0eba16286a2d57ea0991168d4994637e8343e36, + 0x0d51fbc6c71a0094fa2cdd545b11c5c0c797324f1, b163); + +/*----- That's all, folks -------------------------------------------------*/ + diff --git a/calc/gfx.cal b/calc/gfx.cal index 8d8fd00..5b19cb3 100644 --- a/calc/gfx.cal +++ b/calc/gfx.cal @@ -1,6 +1,6 @@ /* -*-apcalc-*- * - * $Id: gfx.cal,v 1.1 2000/10/08 16:01:37 mdw Exp $ + * $Id: gfx.cal,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $ * * Testbed for %$\gf{2}$% poltnomial arithmetic * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: gfx.cal,v $ + * Revision 1.1.4.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.1 2000/10/08 16:01:37 mdw * Prototypes of various bits of code. * @@ -79,7 +82,9 @@ define gf_mul(x, y) define gfx_div(rx, dx) { local r = gfint(rx), d = gfint(dx), i; - local q = 0, dbits = highbit(d), rbits = highbit(r); + local q = 0, dbits, rbits; + dbits = highbit(d); + rbits = highbit(r); for (i = rbits - dbits; i >= 0; i--) { if (bit(r, i + dbits)) { r = xor(r, d << i); @@ -91,14 +96,36 @@ define gfx_div(rx, dx) define gf_div(x, y) { - local l = gfx_div(x, y); + local l; + l = gfx_div(x, y); return gf(l[[0]]); } define gf_mod(x, y) { - local l = gfx_div(x, y); + local l; + l = gfx_div(x, y); return gf(l[[1]]); } +define gf_inv(a, b) +{ + local g, x, y, X, Y, u, v, t, q, r; + x = gf(1); X = gf(0); + y = gf(0); Y = gf(1); + + if (b == gf(0)) { g = a; } else if (a == gf(0)) { g = b; } + else { + while (b != gf(0)) { + q = gf_div(b, a); r = gf_mod(b, a); + t = X * q + x; x = X; X = t; + t = Y * q + y; y = Y; Y = t; + b = a; a = r; + } + g = a; + } + if (g != gf(1)) quit "not coprime in gf_inv"; + return Y; +} + /*----- That's all, folks -------------------------------------------------*/ diff --git a/ec-bin.c b/ec-bin.c new file mode 100644 index 0000000..4f79c3d --- /dev/null +++ b/ec-bin.c @@ -0,0 +1,428 @@ +/* -*-c-*- + * + * $Id: ec-bin.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Arithmetic for elliptic curves over binary fields + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: ec-bin.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "ec.h" + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct ecctx { + ec_curve c; + mp *a, *b; + mp *bb; +} ecctx; + +/*----- Main code ---------------------------------------------------------*/ + +static const ec_ops ec_binops, ec_binprojops; + +static ec *ecneg(ec_curve *c, ec *d, const ec *p) +{ + EC_COPY(d, p); + if (d->x) + d->y = F_ADD(c->f, d->y, d->y, d->x); + return (d); +} + +static ec *ecprojneg(ec_curve *c, ec *d, const ec *p) +{ + EC_COPY(d, p); + if (d->x) { + mp *t = F_MUL(c->f, MP_NEW, d->x, d->z); + d->y = F_ADD(c->f, d->y, d->y, t); + MP_DROP(t); + } + return (d); +} + +static ec *ecfind(ec_curve *c, ec *d, mp *x) +{ + /* write me */ + return (0); +} + +static ec *ecdbl(ec_curve *c, ec *d, const ec *a) +{ + if (EC_ATINF(a) || F_ZEROP(c->f, a->x)) + EC_SETINF(d); + else { + field *f = c->f; + ecctx *cc = (ecctx *)c; + mp *lambda; + mp *dx, *dy; + + dx = F_INV(f, MP_NEW, a->x); /* %$x^{-1}$% */ + dy = F_MUL(f, MP_NEW, dx, a->y); /* %$y/x$% */ + lambda = F_ADD(f, dy, dy, a->x); /* %$\lambda = x + y/x$% */ + + dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */ + dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */ + dx = F_ADD(f, dx, dx, cc->a); /* %$x' = a + \lambda^2 + \lambda$% */ + + dy = F_ADD(f, MP_NEW, a->x, dx); /* %$ x + x' $% */ + dy = F_MUL(f, dy, dy, lambda); /* %$ (x + x') \lambda$% */ + dy = F_ADD(f, dy, dy, a->y); /* %$ (x + x') \lambda + y$% */ + dy = F_ADD(f, dy, dy, dx); /* %$ y' = (x + x') \lambda + y + x'$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = 0; + MP_DROP(lambda); + } + return (d); +} + +static ec *ecprojdbl(ec_curve *c, ec *d, const ec *a) +{ + if (EC_ATINF(a) || F_ZEROP(c->f, a->x)) + EC_SETINF(d); + else { + field *f = c->f; + ecctx *cc = (ecctx *)c; + mp *dx, *dy, *dz, *u, *v; + + dy = F_SQR(f, MP_NEW, a->z); /* %$z^2$% */ + dx = F_MUL(f, MP_NEW, dy, cc->bb); /* %$c z^2$% */ + dx = F_ADD(f, dx, dx, a->x); /* %$x + c z^2$% */ + dz = F_SQR(f, MP_NEW, dx); /* %$(x + c z^2)^2$% */ + dx = F_SQR(f, dx, dz); /* %$x' = (x + c z^2)^4$% */ + + dz = F_MUL(f, dz, dy, a->x); /* %$z' = x z^2$% */ + + dy = F_SQR(f, dy, a->x); /* %$x^2$% */ + u = F_MUL(f, MP_NEW, a->y, a->z); /* %$y z$% */ + u = F_ADD(f, u, u, dz); /* %$z' + y z$% */ + u = F_ADD(f, u, u, dy); /* %$u = z' + x^2 + y z$% */ + + v = F_SQR(f, MP_NEW, dy); /* %$x^4$% */ + dy = F_MUL(f, dy, v, dz); /* %$x^4 z'$% */ + v = F_MUL(f, v, u, dx); /* %$u x'$% */ + dy = F_ADD(f, dy, dy, v); /* %$y' = x^4 z' + u x'$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = dz; + MP_DROP(u); + MP_DROP(v); + assert(!(d->x->f & MP_DESTROYED)); + assert(!(d->y->f & MP_DESTROYED)); + assert(!(d->z->f & MP_DESTROYED)); + } + return (d); +} + +static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b) +{ + if (a == b) + ecdbl(c, d, a); + else if (EC_ATINF(a)) + EC_COPY(d, b); + else if (EC_ATINF(b)) + EC_COPY(d, a); + else { + field *f = c->f; + ecctx *cc = (ecctx *)c; + mp *lambda; + mp *dx, *dy; + + if (!MP_EQ(a->x, b->x)) { + dx = F_ADD(f, MP_NEW, a->x, b->x); /* %$x_0 + x_1$% */ + dy = F_INV(f, MP_NEW, dx); /* %$(x_0 + x_1)^{-1}$% */ + dx = F_ADD(f, dx, a->y, b->y); /* %$y_0 + y_1$% */ + lambda = F_MUL(f, MP_NEW, dy, dx); + /* %$\lambda = (y_0 + y_1)/(x_0 + x_1)$% */ + + dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */ + dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */ + dx = F_ADD(f, dx, dx, cc->a); /* %$a + \lambda^2 + \lambda$% */ + dx = F_ADD(f, dx, dx, a->x); /* %$a + \lambda^2 + \lambda + x_0$% */ + dx = F_ADD(f, dx, dx, b->x); + /* %$x' = a + \lambda^2 + \lambda + x_0 + x_1$% */ + } else if (!MP_EQ(a->y, b->y) || F_ZEROP(f, a->x)) { + EC_SETINF(d); + return (d); + } else { + dx = F_INV(f, MP_NEW, a->x); /* %$x^{-1}$% */ + dy = F_MUL(f, MP_NEW, dx, a->y); /* %$y/x$% */ + lambda = F_ADD(f, dy, dy, a->x); /* %$\lambda = x + y/x$% */ + + dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */ + dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */ + dx = F_ADD(f, dx, dx, cc->a); /* %$x' = a + \lambda^2 + \lambda$% */ + dy = MP_NEW; + } + + dy = F_ADD(f, dy, a->x, dx); /* %$ x + x' $% */ + dy = F_MUL(f, dy, dy, lambda); /* %$ (x + x') \lambda$% */ + dy = F_ADD(f, dy, dy, a->y); /* %$ (x + x') \lambda + y$% */ + dy = F_ADD(f, dy, dy, dx); /* %$ y' = (x + x') \lambda + y + x'$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = 0; + MP_DROP(lambda); + } + return (d); +} + +static ec *ecprojadd(ec_curve *c, ec *d, const ec *a, const ec *b) +{ + if (a == b) + c->ops->dbl(c, d, a); + else if (EC_ATINF(a)) + EC_COPY(d, b); + else if (EC_ATINF(b)) + EC_COPY(d, a); + else { + field *f = c->f; + ecctx *cc = (ecctx *)c; + mp *dx, *dy, *dz, *u, *uu, *v, *t, *s, *ss, *r, *w, *l; + + dz = F_SQR(f, MP_NEW, b->z); /* %$z_1^2$% */ + u = F_MUL(f, MP_NEW, dz, a->x); /* %$u_0 = x_0 z_1^2$% */ + t = F_MUL(f, MP_NEW, dz, b->z); /* %$z_1^3$% */ + s = F_MUL(f, MP_NEW, t, a->y); /* %$s_0 = y_0 z_1^3$% */ + + dz = F_SQR(f, dz, a->z); /* %$z_0^2$% */ + uu = F_MUL(f, MP_NEW, dz, b->x); /* %$u_1 = x_1 z_0^2$% */ + t = F_MUL(f, t, dz, a->z); /* %$z_0^3$% */ + ss = F_MUL(f, MP_NEW, t, b->y); /* %$s_1 = y_1 z_0^3$% */ + + w = F_ADD(f, u, u, uu); /* %$r = u_0 + u_1$% */ + r = F_ADD(f, s, s, ss); /* %$w = s_0 + s_1$% */ + if (F_ZEROP(f, w)) { + MP_DROP(w); + MP_DROP(uu); + MP_DROP(ss); + MP_DROP(t); + MP_DROP(dz); + if (F_ZEROP(f, r)) { + MP_DROP(r); + return (c->ops->dbl(c, d, a)); + } else { + MP_DROP(r); + EC_SETINF(d); + return (d); + } + } + + l = F_MUL(f, t, a->z, w); /* %$l = z_0 w$% */ + + dz = F_MUL(f, dz, l, b->z); /* %$z' = l z_1$% */ + + ss = F_MUL(f, ss, r, b->x); /* %$r x_1$% */ + t = F_MUL(f, uu, l, b->y); /* %$l y_1$% */ + v = F_ADD(f, ss, ss, t); /* %$v = r x_1 + l y_1$% */ + + t = F_ADD(f, t, r, dz); /* %$t = r + z'$% */ + + uu = F_SQR(f, MP_NEW, dz); /* %$z'^2$% */ + dx = F_MUL(f, MP_NEW, uu, cc->a); /* %$a z'^2$% */ + uu = F_MUL(f, uu, t, r); /* %$t r$% */ + dx = F_ADD(f, dx, dx, uu); /* %$a z'^2 + t r$% */ + r = F_SQR(f, r, w); /* %$w^2$% */ + uu = F_MUL(f, uu, r, w); /* %$w^3$% */ + dx = F_ADD(f, dx, dx, uu); /* %$x' = a z'^2 + t r + w^3$% */ + + r = F_SQR(f, r, l); /* %$l^2$% */ + dy = F_MUL(f, uu, v, r); /* %$v l^2$% */ + l = F_MUL(f, l, t, dx); /* %$t x'$% */ + dy = F_ADD(f, dy, dy, l); /* %$y' = t x' + v l^2$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = dz; + MP_DROP(l); + MP_DROP(r); + MP_DROP(w); + MP_DROP(t); + MP_DROP(v); + } + return (d); +} + +static int eccheck(ec_curve *c, const ec *p) +{ + ecctx *cc = (ecctx *)c; + field *f = c->f; + int rc; + mp *u, *v; + + v = F_SQR(f, MP_NEW, p->x); + u = F_MUL(f, MP_NEW, v, p->x); + v = F_MUL(f, v, v, cc->a); + u = F_ADD(f, u, u, v); + u = F_ADD(f, u, u, cc->b); + v = F_MUL(f, v, p->x, p->y); + u = F_ADD(f, u, u, v); + v = F_SQR(f, v, p->y); + u = F_ADD(f, u, u, v); + rc = F_ZEROP(f, u); + mp_drop(u); + mp_drop(v); + return (rc); +} + +static int ecprojcheck(ec_curve *c, const ec *p) +{ + ec t = EC_INIT; + int rc; + + c->ops->fix(c, &t, p); + rc = eccheck(c, &t); + EC_DESTROY(&t); + return (rc); +} + +static void ecdestroy(ec_curve *c) +{ + ecctx *cc = (ecctx *)c; + MP_DROP(cc->a); + MP_DROP(cc->b); + if (cc->bb) MP_DROP(cc->bb); + DESTROY(cc); +} + +/* --- @ec_bin@, @ec_binproj@ --- * + * + * Arguments: @field *f@ = the underlying field for this elliptic curve + * @mp *a, *b@ = the coefficients for this curve + * + * Returns: A pointer to the curve. + * + * Use: Creates a curve structure for an elliptic curve defined over + * a binary field. The @binproj@ variant uses projective + * coordinates, which can be a win. + */ + +ec_curve *ec_bin(field *f, mp *a, mp *b) +{ + ecctx *cc = CREATE(ecctx); + cc->c.ops = &ec_binops; + cc->c.f = f; + cc->a = F_IN(f, MP_NEW, a); + cc->b = F_IN(f, MP_NEW, b); + cc->bb = 0; + return (&cc->c); +} + +ec_curve *ec_binproj(field *f, mp *a, mp *b) +{ + ecctx *cc = CREATE(ecctx); + cc->c.ops = &ec_binprojops; + cc->c.f = f; + cc->a = F_IN(f, MP_NEW, a); + cc->b = F_IN(f, MP_NEW, b); + cc->bb = F_SQRT(f, MP_NEW, b); + cc->bb = F_SQRT(f, cc->bb, cc->bb); + return (&cc->c); +} + +static const ec_ops ec_binops = { + ecdestroy, ec_idin, ec_idout, ec_idfix, + 0, ecneg, ecadd, ec_stdsub, ecdbl, eccheck +}; + +static const ec_ops ec_binprojops = { + ecdestroy, ec_projin, ec_projout, ec_projfix, + 0, ecprojneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck +}; + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#define MP(x) mp_readstring(MP_NEW, #x, 0, 0) + +int main(int argc, char *argv[]) +{ + field *f; + ec_curve *c; + ec g = EC_INIT, d = EC_INIT; + mp *p, *a, *b, *r; + int i, n = argc == 1 ? 1 : atoi(argv[1]); + + printf("ec-bin: "); + fflush(stdout); + a = MP(1); + b = MP(0x066647ede6c332c7f8c0923bb58213b333b20e9ce4281fe115f7d8f90ad); + p = MP(0x20000000000000000000000000000000000000004000000000000000001); + r = + MP(6901746346790563787434755862277025555839812737345013555379383634485462); + + f = field_binpoly(p); + c = ec_binproj(f, a, b); + + g.x = MP(0x0fac9dfcbac8313bb2139f1bb755fef65bc391f8b36f8f8eb7371fd558b); + g.y = MP(0x1006a08a41903350678e58528bebf8a0beff867a7ca36716f7e01f81052); + + for (i = 0; i < n; i++) { + ec_mul(c, &d, &g, r); + if (EC_ATINF(&d)) { + fprintf(stderr, "zero too early\n"); + return (1); + } + ec_add(c, &d, &d, &g); + if (!EC_ATINF(&d)) { + fprintf(stderr, "didn't reach zero\n"); + MP_EPRINTX("d.x", d.x); + MP_EPRINTX("d.y", d.y); + MP_EPRINTX("d.z", d.y); + return (1); + } + ec_destroy(&d); + } + + ec_destroy(&g); + ec_destroycurve(c); + F_DESTROY(f); + MP_DROP(p); MP_DROP(a); MP_DROP(b); MP_DROP(r); + assert(!mparena_count(&mparena_global)); + printf("ok\n"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/ec-prime.c b/ec-prime.c index 14e4c16..40f487e 100644 --- a/ec-prime.c +++ b/ec-prime.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec-prime.c,v 1.3.4.2 2004/03/20 00:13:31 mdw Exp $ + * $Id: ec-prime.c,v 1.3.4.3 2004/03/21 22:39:46 mdw Exp $ * * Elliptic curves over prime fields * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec-prime.c,v $ + * Revision 1.3.4.3 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.3.4.2 2004/03/20 00:13:31 mdw * Projective coordinates for prime curves * @@ -67,7 +70,8 @@ static const ec_ops ec_primeops, ec_primeprojops, ec_primeprojxops; static ec *ecneg(ec_curve *c, ec *d, const ec *p) { EC_COPY(d, p); - d->y = F_NEG(c->f, d->y, d->y); + if (d->y) + d->y = F_NEG(c->f, d->y, d->y); return (d); } @@ -254,8 +258,7 @@ static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b) dx = F_SUB(f, dx, dx, b->x); /* %$x' = \lambda^2 - x_0 - x_1$% */ dy = F_SUB(f, dy, b->x, dx); /* %$x_1 - x'$% */ dy = F_MUL(f, dy, lambda, dy); /* %$\lambda (x_1 - x')$% */ - dy = F_SUB(f, dy, dy, b->y); - /* %$y' = \lambda (x_1 - x') - y_1$% */ + dy = F_SUB(f, dy, dy, b->y); /* %$y' = \lambda (x_1 - x') - y_1$% */ EC_DESTROY(d); d->x = dx; @@ -286,17 +289,14 @@ static ec *ecprojadd(ec_curve *c, ec *d, const ec *a, const ec *b) w = F_SUB(f, p, a->x, u); /* %$w = x_0 - u$% */ r = F_SUB(f, MP_NEW, a->y, s); /* %$r = y_0 - s$% */ if (F_ZEROP(f, w)) { + MP_DROP(w); + MP_DROP(u); + MP_DROP(s); if (F_ZEROP(f, r)) { - MP_DROP(w); MP_DROP(r); - MP_DROP(u); - MP_DROP(s); return (c->ops->dbl(c, d, a)); } else { - MP_DROP(w); MP_DROP(r); - MP_DROP(u); - MP_DROP(s); EC_SETINF(d); return (d); } @@ -430,12 +430,13 @@ static const ec_ops ec_primeprojxops = { #define MP(x) mp_readstring(MP_NEW, #x, 0, 0) -int main(void) +int main(int argc, char *argv[]) { field *f; ec_curve *c; ec g = EC_INIT, d = EC_INIT; mp *p, *a, *b, *r; + int i, n = argc == 1 ? 1 : atoi(argv[1]); printf("ec-prime: "); fflush(stdout); @@ -445,25 +446,26 @@ int main(void) r = MP(6277101735386680763835789423176059013767194773182842284080); f = field_prime(p); - c = ec_prime(f, a, b); + c = ec_primeproj(f, a, b); g.x = MP(0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012); g.y = MP(0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811); - ec_mul(c, &d, &g, r); - if (EC_ATINF(&d)) { - fprintf(stderr, "zero too early\n"); - return (1); - } - ec_add(c, &d, &d, &g); - if (!EC_ATINF(&d)) { - fprintf(stderr, "didn't reach zero\n"); - MP_EPRINT("d.x", d.x); - MP_EPRINT("d.y", d.y); - return (1); + for (i = 0; i < n; i++) { + ec_mul(c, &d, &g, r); + if (EC_ATINF(&d)) { + fprintf(stderr, "zero too early\n"); + return (1); + } + ec_add(c, &d, &d, &g); + if (!EC_ATINF(&d)) { + fprintf(stderr, "didn't reach zero\n"); + MP_EPRINT("d.x", d.x); + MP_EPRINT("d.y", d.y); + return (1); + } + ec_destroy(&d); } - - ec_destroy(&d); ec_destroy(&g); ec_destroycurve(c); F_DESTROY(f); diff --git a/ec.h b/ec.h index 476d898..bb5e201 100644 --- a/ec.h +++ b/ec.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec.h,v 1.4.4.2 2004/03/20 00:13:31 mdw Exp $ + * $Id: ec.h,v 1.4.4.3 2004/03/21 22:39:46 mdw Exp $ * * Elliptic curve definitions * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec.h,v $ + * Revision 1.4.4.3 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.4.4.2 2004/03/20 00:13:31 mdw * Projective coordinates for prime curves * @@ -419,18 +422,20 @@ extern void ec_destroycurve(ec_curve */*c*/); extern ec_curve *ec_prime(field */*f*/, mp */*a*/, mp */*b*/); extern ec_curve *ec_primeproj(field */*f*/, mp */*a*/, mp */*b*/); -/* --- @ec_bin@ --- * +/* --- @ec_bin@, @ec_binproj@ --- * * * Arguments: @field *f@ = the underlying field for this elliptic curve * @mp *a, *b@ = the coefficients for this curve * * Returns: A pointer to the curve. * - * Use: Creates a curve structure for a non-supersingular elliptic - * curve defined over a binary field. + * Use: Creates a curve structure for an elliptic curve defined over + * a binary field. The @binproj@ variant uses projective + * coordinates, which can be a win. */ extern ec_curve *ec_bin(field */*f*/, mp */*a*/, mp */*b*/); +extern ec_curve *ec_binproj(field */*f*/, mp */*a*/, mp */*b*/); /*----- That's all, folks -------------------------------------------------*/ diff --git a/f-binpoly.c b/f-binpoly.c new file mode 100644 index 0000000..509efc4 --- /dev/null +++ b/f-binpoly.c @@ -0,0 +1,142 @@ +/* -*-c-*- + * + * $Id: f-binpoly.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Binary fields with polynomial basis representation + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: f-binpoly.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "field.h" +#include "gf.h" +#include "gfreduce.h" + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct fctx { + field f; + gfreduce r; +} fctx; + +/*----- Main code ---------------------------------------------------------*/ + +/* --- Field operations --- */ + +static void fdestroy(field *ff) +{ + fctx *f = (fctx *)ff; + gfreduce_destroy(&f->r); + DESTROY(f); +} + +static int fzerop(field *ff, mp *x) +{ + return (!MP_LEN(x)); +} + +static mp *fadd(field *ff, mp *d, mp *x, mp *y) +{ + return (gf_add(d, x, y)); +} + +static mp *fmul(field *ff, mp *d, mp *x, mp *y) +{ + fctx *f = (fctx *)ff; + d = gf_mul(d, x, y); + return (gfreduce_do(&f->r, d, d)); +} + +static mp *fsqr(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + d = gf_sqr(d, x); + return (gfreduce_do(&f->r, d, d)); +} + +static mp *finv(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + gf_gcd(0, 0, &d, f->r.p, x); + return (d); +} + +static mp *freduce(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + return (gfreduce_do(&f->r, d, x)); +} + +static mp *fsqrt(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + return (gfreduce_sqrt(&f->r, d, x)); +} + +static mp *fquadsolve(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + return (gfreduce_quadsolve(&f->r, d, x)); +} + +/* --- Field operations table --- */ + +static field_ops fops = { + fdestroy, + freduce, field_id, + fzerop, field_id, fadd, fadd, fmul, fsqr, finv, freduce, fsqrt, + fquadsolve, + 0, 0, 0, 0 +}; + +/* --- @field_binpoly@ --- * + * + * Arguments: @mp *p@ = the reduction polynomial + * + * Returns: A pointer to the field. + * + * Use: Creates a field structure for a binary field mod @p@. + */ + +field *field_binpoly(mp *p) +{ + fctx *f = CREATE(fctx); + f->f.ops = &fops; + f->f.zero = MP_ZERO; + f->f.one = MP_ONE; + gfreduce_create(&f->r, p); + return (&f->f); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/f-prime.c b/f-prime.c index 7c1dae5..7215ec8 100644 --- a/f-prime.c +++ b/f-prime.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: f-prime.c,v 1.3.4.2 2004/03/20 00:13:31 mdw Exp $ + * $Id: f-prime.c,v 1.3.4.3 2004/03/21 22:39:46 mdw Exp $ * * Prime fields with Montgomery arithmetic * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: f-prime.c,v $ + * Revision 1.3.4.3 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.3.4.2 2004/03/20 00:13:31 mdw * Projective coordinates for prime curves * @@ -146,6 +149,16 @@ static mp *freduce(field *ff, mp *d, mp *x) return (d); } +static mp *fsqrt(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + d = mpmont_reduce(&f->mm, d, x); + d = mp_modsqrt(d, d, f->mm.m); + if (!d) + return (d); + return (mpmont_mul(&f->mm, d, d, f->mm.r2)); +} + static mp *fdbl(field *ff, mp *d, mp *x) { fctx *f = (fctx *)ff; @@ -189,23 +202,14 @@ static mp *fhlv(field *ff, mp *d, mp *x) return (mp_lsr(d, x, 1)); } -static mp *fsqrt(field *ff, mp *d, mp *x) -{ - fctx *f = (fctx *)ff; - d = mpmont_reduce(&f->mm, d, x); - d = mp_modsqrt(d, d, f->mm.m); - if (!d) - return (d); - return (mpmont_mul(&f->mm, d, d, f->mm.r2)); -} - /* --- Field operations table --- */ static field_ops fops = { fdestroy, fin, fout, - fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, - fdbl, ftpl, fqdl, fhlv, fsqrt + fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, fsqrt, + 0, + fdbl, ftpl, fqdl, fhlv }; /* --- @field_prime@ --- * diff --git a/field.h b/field.h index 8096700..ea019c5 100644 --- a/field.h +++ b/field.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: field.h,v 1.3.4.1 2004/03/20 00:13:31 mdw Exp $ + * $Id: field.h,v 1.3.4.2 2004/03/21 22:39:46 mdw Exp $ * * Definitions for field arithmetic * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: field.h,v $ + * Revision 1.3.4.2 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.3.4.1 2004/03/20 00:13:31 mdw * Projective coordinates for prime curves * @@ -81,6 +84,11 @@ typedef struct field_ops { mp *(*sqr)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*inv)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*reduce)(field */*f*/, mp */*d*/, mp */*x*/); + mp *(*sqrt)(field */*f*/, mp */*d*/, mp */*x*/); + + /* --- Operations for binary fields only --- */ + + mp *(*quadsolve)(field */*f*/, mp */*d*/, mp */*x*/); /* --- Operations for prime fields only --- */ @@ -88,7 +96,6 @@ typedef struct field_ops { mp *(*tpl)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*qdl)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*hlv)(field */*f*/, mp */*d*/, mp */*x*/); - mp *(*sqrt)(field */*f*/, mp */*d*/, mp */*x*/); } field_ops; @@ -105,12 +112,14 @@ typedef struct field_ops { #define F_SQR(f, d, x) (f)->ops->sqr((f), (d), (x)) #define F_INV(f, d, x) (f)->ops->inv((f), (d), (x)) #define F_REDUCE(f, d, x) (f)->ops->reduce((f), (d), (x)) +#define F_SQRT(f, d, x) (f)->ops->sqrt((f), (d), (x)) + +#define F_QUADSOLVE(f, d, x) (f)->ops->quadsolve((f), (d), (x)) #define F_DBL(f, d, x) (f)->ops->dbl((f), (d), (x)) #define F_TPL(f, d, x) (f)->ops->tpl((f), (d), (x)) #define F_QDL(f, d, x) (f)->ops->qdl((f), (d), (x)) #define F_HLV(f, d, x) (f)->ops->hlv((f), (d), (x)) -#define F_SQRT(f, d, x) (f)->ops->sqrt((f), (d), (x)) /*----- Helpful field operations ------------------------------------------*/ @@ -154,20 +163,6 @@ extern field *field_prime(mp */*p*/); extern field *field_binpoly(mp */*p*/); -/* --- @field_binsparsepoly@ --- * - * - * Arguments: @unsigned deg@ = the degree of the field polynomial - * @const unsigned *c@ = degrees of nonzero coefficients - * - * Returns: A pointer to the field. - * - * Use: Creates a field structure for a binary field taking advantage - * of the sparseness of the given polynomial to improve - * reduction performance. - */ - -extern field *field_binsparsepoly(unsigned /*deg*/, const unsigned */*c*/); - /*----- That's all, folks -------------------------------------------------*/ #ifdef __cplusplus diff --git a/gf-arith.c b/gf-arith.c new file mode 100644 index 0000000..6838e44 --- /dev/null +++ b/gf-arith.c @@ -0,0 +1,263 @@ +/* -*-c-*- + * + * $Id: gf-arith.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Basic arithmetic on binary polynomials + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gf-arith.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "gf.h" + +/*----- Macros ------------------------------------------------------------*/ + +#define MAX(x, y) ((x) >= (y) ? (x) : (y)) + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gf_add@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a, *b@ = sources + * + * Returns: Result, @a@ added to @b@. + */ + +mp *gf_add(mp *d, mp *a, mp *b) +{ + MP_DEST(d, MAX(MP_LEN(a), MP_LEN(b)), (a->f | b->f) & MP_BURN); + gfx_add(d->v, d->vl, a->v, a->vl, b->v, b->vl); + d->f = (a->f | b->f) & MP_BURN; + MP_SHRINK(d); + return (d); +} + +/* --- @gf_mul@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a, *b@ = sources + * + * Returns: Result, @a@ multiplied by @b@. + */ + +mp *gf_mul(mp *d, mp *a, mp *b) +{ + a = MP_COPY(a); + b = MP_COPY(b); + + if (MP_LEN(a) <= MPK_THRESH || MP_LEN(b) <= GFK_THRESH) { + MP_DEST(d, MP_LEN(a) + MP_LEN(b), a->f | b->f | MP_UNDEF); + gfx_mul(d->v, d->vl, a->v, a->vl, b->v, b->vl); + } else { + size_t m = MAX(MP_LEN(a), MP_LEN(b)); + mpw *s; + MP_DEST(d, 2 * m, a->f | b->f | MP_UNDEF); + s = mpalloc(d->a, 2 * m); + gfx_kmul(d->v, d->vl, a->v, a->vl, b->v, b->vl, s, s + 2 * m); + mpfree(d->a, s); + } + + d->f = (a->f | b->f) & MP_BURN; + MP_SHRINK(d); + MP_DROP(a); + MP_DROP(b); + return (d); +} + +/* --- @gf_sqr@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a@ = source + * + * Returns: Result, @a@ squared. + */ + +mp *gf_sqr(mp *d, mp *a) +{ + MP_COPY(a); + MP_DEST(d, 2 * MP_LEN(a), a->f & MP_BURN); + gfx_sqr(d->v, d->vl, a->v, a->vl); + d->f = a->f & MP_BURN; + MP_SHRINK(d); + MP_DROP(a); + return (d); +} + +/* --- @gf_div@ --- * + * + * Arguments: @mp **qq, **rr@ = destination, quotient and remainder + * @mp *a, *b@ = sources + * + * Use: Calculates the quotient and remainder when @a@ is divided by + * @b@. The destinations @*qq@ and @*rr@ must be distinct. + * Either of @qq@ or @rr@ may be null to indicate that the + * result is irrelevant. (Discarding both results is silly.) + * There is a performance advantage if @a == *rr@. + */ + +void gf_div(mp **qq, mp **rr, mp *a, mp *b) + { + mp *r = rr ? *rr : MP_NEW; + mp *q = qq ? *qq : MP_NEW; + + /* --- Set the remainder up right --- */ + + b = MP_COPY(b); + a = MP_COPY(a); + if (r) + MP_DROP(r); + r = a; + MP_DEST(r, MP_LEN(b) + 2, a->f | b->f); + + /* --- Fix up the quotient too --- */ + + r = MP_COPY(r); + MP_DEST(q, MP_LEN(r), r->f | MP_UNDEF); + MP_DROP(r); + + /* --- Perform the calculation --- */ + + gfx_div(q->v, q->vl, r->v, r->vl, b->v, b->vl); + + /* --- Sort out the sign of the results --- * + * + * If the signs of the arguments differ, and the remainder is nonzero, I + * must add one to the absolute value of the quotient and subtract the + * remainder from @b@. + */ + + q->f = (r->f | b->f) & MP_BURN; + r->f = (r->f | b->f) & MP_BURN; + + /* --- Store the return values --- */ + + MP_DROP(b); + + if (!qq) + MP_DROP(q); + else { + MP_SHRINK(q); + *qq = q; + } + + if (!rr) + MP_DROP(r); + else { + MP_SHRINK(r); + *rr = r; + } +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b) +{ + if (!MP_EQ(expect, result)) { + fprintf(stderr, "\n*** %s failed", op); + fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 16); + fputs("\n*** b = ", stderr); mp_writefile(b, stderr, 16); + fputs("\n*** result = ", stderr); mp_writefile(result, stderr, 16); + fputs("\n*** expect = ", stderr); mp_writefile(expect, stderr, 16); + fputc('\n', stderr); + return (0); + } + return (1); +} + +#define RIG(name, op) \ + static int t##name(dstr *v) \ + { \ + mp *a = *(mp **)v[0].buf; \ + mp *b = *(mp **)v[1].buf; \ + mp *r = *(mp **)v[2].buf; \ + mp *c = op(MP_NEW, a, b); \ + int ok = verify(#name, r, c, a, b); \ + mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(r); \ + assert(mparena_count(MPARENA_GLOBAL) == 0); \ + return (ok); \ + } + +RIG(add, gf_add) +RIG(mul, gf_mul) + +#undef RIG + +static int tsqr(dstr *v) +{ + mp *a = *(mp **)v[0].buf; + mp *r = *(mp **)v[1].buf; + mp *c = MP_NEW; + int ok = 1; + c = gf_sqr(MP_NEW, a); + ok &= verify("sqr", r, c, a, MP_ZERO); + mp_drop(a); mp_drop(r); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int tdiv(dstr *v) +{ + mp *a = *(mp **)v[0].buf; + mp *b = *(mp **)v[1].buf; + mp *q = *(mp **)v[2].buf; + mp *r = *(mp **)v[3].buf; + mp *c = MP_NEW, *d = MP_NEW; + int ok = 1; + gf_div(&c, &d, a, b); + ok &= verify("div(quotient)", q, c, a, b); + ok &= verify("div(remainder)", r, d, a, b); + mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(d); mp_drop(r); mp_drop(q); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk tests[] = { + { "add", tadd, { &type_mp, &type_mp, &type_mp, 0 } }, + { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } }, + { "sqr", tsqr, { &type_mp, &type_mp, 0 } }, + { "div", tdiv, { &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/gf"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gf-gcd.c b/gf-gcd.c new file mode 100644 index 0000000..64d61d4 --- /dev/null +++ b/gf-gcd.c @@ -0,0 +1,259 @@ +/* -*-c-*- + * + * $Id: gf-gcd.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Euclidian algorithm on binary polynomials + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gf-gcd.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "gf.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gf_gcd@ --- * + * + * Arguments: @mp **gcd, **xx, **yy@ = where to write the results + * @mp *a, *b@ = sources (must be nonzero) + * + * + * Returns: --- + * + * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that + * @ax + by = gcd(a, b)@. This is useful for computing modular + * inverses. + */ + +void gf_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) +{ + mp *x = MP_ONE, *X = MP_ZERO; + mp *y = MP_ZERO, *Y = MP_ONE; + mp *u, *v; + mp *q = MP_NEW; + unsigned f = 0; + +#define f_swap 1u +#define f_ext 2u + + /* --- Sort out some initial flags --- */ + + if (xx || yy) + f |= f_ext; + + /* --- Ensure that @a@ is larger than @b@ --- * + * + * If they're the same length we don't care which order they're in, so this + * unsigned comparison is fine. + */ + + if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) { + { mp *t = a; a = b; b = t; } + f |= f_swap; + } + + /* --- Check for zeroness --- */ + + if (MP_EQ(b, MP_ZERO)) { + + /* --- Store %$|a|$% as the GCD --- */ + + if (gcd) { + if (*gcd) MP_DROP(*gcd); + a = MP_COPY(a); + *gcd = a; + } + + /* --- Store %$1$% and %$0$% in the appropriate bins --- */ + + if (f & f_ext) { + if (f & f_swap) { + mp **t = xx; xx = yy; yy = t; + } + if (xx) { + if (*xx) MP_DROP(*xx); + if (MP_EQ(a, MP_ZERO)) + *xx = MP_ZERO; + else + *xx = MP_ONE; + } + if (yy) { + if (*yy) MP_DROP(*yy); + *yy = MP_ZERO; + } + } + return; + } + + /* --- Take a reference to the arguments --- */ + + a = MP_COPY(a); + b = MP_COPY(b); + + /* --- Make sure @a@ and @b@ are not both even --- */ + + MP_SPLIT(a); a->f &= ~MP_NEG; + MP_SPLIT(b); b->f &= ~MP_NEG; + + u = MP_COPY(a); + v = MP_COPY(b); + + while (MP_LEN(v)) { + mp *t; + gf_div(&q, &u, u, v); + if (f & f_ext) { + t = gf_mul(MP_NEW, X, q); + t = gf_add(t, x, t); + MP_DROP(x); x = X; X = t; + t = gf_mul(MP_NEW, Y, q); + t = gf_add(t, y, t); + MP_DROP(y); y = Y; Y = t; + } + t = u; u = v; v = t; + } + + MP_DROP(q); + if (!gcd) + MP_DROP(u); + else { + if (*gcd) MP_DROP(*gcd); + u->f &= ~MP_NEG; + *gcd = u; + } + + /* --- Perform a little normalization --- */ + + if (f & f_ext) { + + /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */ + + if (f & f_swap) { + mp *t = x; x = y; y = t; + t = a; a = b; b = t; + } + + /* --- Store the results --- */ + + if (!xx) + MP_DROP(x); + else { + if (*xx) MP_DROP(*xx); + *xx = x; + } + + if (!yy) + MP_DROP(y); + else { + if (*yy) MP_DROP(*yy); + *yy = y; + } + } + + MP_DROP(v); + MP_DROP(X); MP_DROP(Y); + MP_DROP(a); MP_DROP(b); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int gcd(dstr *v) +{ + int ok = 1; + mp *a = *(mp **)v[0].buf; + mp *b = *(mp **)v[1].buf; + mp *g = *(mp **)v[2].buf; + mp *x = *(mp **)v[3].buf; + mp *y = *(mp **)v[4].buf; + + mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW; + gf_gcd(&gg, &xx, &yy, a, b); + if (!MP_EQ(x, xx)) { + fputs("\n*** mp_gcd(x) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + if (!MP_EQ(y, yy)) { + fputs("\n*** mp_gcd(y) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + + if (!ok) { + mp *ax = gf_mul(MP_NEW, a, xx); + mp *by = gf_mul(MP_NEW, b, yy); + ax = gf_add(ax, ax, by); + if (MP_EQ(ax, gg)) + fputs("\n*** (Alternative result found.)\n", stderr); + MP_DROP(ax); + MP_DROP(by); + } + + if (!MP_EQ(g, gg)) { + fputs("\n*** mp_gcd(gcd) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y); + MP_DROP(gg); MP_DROP(xx); MP_DROP(yy); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk tests[] = { + { "gcd", gcd, { &type_mp, &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/gf"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gf.h b/gf.h new file mode 100644 index 0000000..889cd9b --- /dev/null +++ b/gf.h @@ -0,0 +1,124 @@ +/* -*-c-*- + * + * $Id: gf.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Arithmetic on binary polynomials + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gf.h,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +#ifndef CATACOMB_GF_H +#define CATACOMB_GF_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef CATACOMB_MP_H +# include "mp.h" +#endif + +#ifndef CATACOMB_GFX_H +# include "gfx.h" +#endif + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @gf_add@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a, *b@ = sources + * + * Returns: Result, @a@ added to @b@. + */ + +extern mp *gf_add(mp */*d*/, mp */*a*/, mp */*b*/); +#define gf_sub gf_add + +/* --- @gf_mul@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a, *b@ = sources + * + * Returns: Result, @a@ multiplied by @b@. + */ + +extern mp *gf_mul(mp */*d*/, mp */*a*/, mp */*b*/); + +/* --- @gf_sqr@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a@ = source + * + * Returns: Result, @a@ squared. + */ + +extern mp *gf_sqr(mp */*d*/, mp */*a*/); + +/* --- @gf_div@ --- * + * + * Arguments: @mp **qq, **rr@ = destination, quotient and remainder + * @mp *a, *b@ = sources + * + * Use: Calculates the quotient and remainder when @a@ is divided by + * @b@. The destinations @*qq@ and @*rr@ must be distinct. + * Either of @qq@ or @rr@ may be null to indicate that the + * result is irrelevant. (Discarding both results is silly.) + * There is a performance advantage if @a == *rr@. + */ + +extern void gf_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/); + +/* --- @gf_gcd@ --- * + * + * Arguments: @mp **gcd, **xx, **yy@ = where to write the results + * @mp *a, *b@ = sources (must be nonzero) + * + * + * Returns: --- + * + * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that + * @ax + by = gcd(a, b)@. This is useful for computing modular + * inverses. + */ + +extern void gf_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/, + mp */*a*/, mp */*b*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/gfreduce-exp.h b/gfreduce-exp.h new file mode 100644 index 0000000..0393145 --- /dev/null +++ b/gfreduce-exp.h @@ -0,0 +1,84 @@ +/* -*-c-*- + * + * $Id: gfreduce-exp.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Exponentiation operations for binary field reduction + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gfreduce-exp.h,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +#ifndef CATACOMB_GFREDUCE_EXP_H +#define CATACOMB_GFREDUCE_EXP_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Exponentiation definitions ----------------------------------------*/ + +#define EXP_TYPE mp * + +#define EXP_COPY(d, x) d = MP_COPY(x) +#define EXP_DROP(x) MP_DROP(x) + +#define EXP_MUL(a, x) do { \ + mp *t = gf_mul(spare, a, x); \ + spare = a; \ + a = gfreduce_do(gr, t, t); \ +} while (0) + +#define EXP_SQR(a) do { \ + mp *t = gf_sqr(spare, a); \ + spare = a; \ + a = gfreduce_do(gr, t, t); \ +} while (0) + +#define EXP_FIX(x) + +#define EXP_SETMUL(d, x, y) do { \ + d = gf_mul(MP_NEW, x, y); \ + d = gfreduce_do(gr, d, d); \ +} while (0) + +#define EXP_SETSQR(d, x) do { \ + d = gf_sqr(MP_NEW, x); \ + d = gfreduce_do(gr, d, d); \ +} while (0) + +#include "exp.h" + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/gfreduce.c b/gfreduce.c new file mode 100644 index 0000000..3969f11 --- /dev/null +++ b/gfreduce.c @@ -0,0 +1,652 @@ +/* -*-c-*- + * + * $Id: gfreduce.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Efficient reduction modulo sparse binary polynomials + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gfreduce.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include + +#include "gf.h" +#include "gfreduce.h" +#include "gfreduce-exp.h" +#include "fibrand.h" +#include "mprand.h" + +/*----- Data structures ---------------------------------------------------*/ + +DA_DECL(instr_v, gfreduce_instr); + +/*----- Main code ---------------------------------------------------------*/ + +/* --- What's going on here? --- * + * + * Let's face it, @gfx_div@ sucks. It works (I hope), but it's not in any + * sense fast. Here, we do efficient reduction modulo sparse polynomials. + * + * Suppose we have a polynomial @X@ we're trying to reduce mod @P@. If we + * take the topmost nonzero word of @X@, call it @w@, then we can eliminate + * it by subtracting off @w P x^{k}@ for an appropriate value of @k@. The + * trick is in observing that if @P@ is sparse we can do this multiplication + * and subtraction efficiently, just by XORing appropriate shifts of @w@ into + * @X@. + * + * The first tricky bit is in working out when to stop. I'll use eight-bit + * words to demonstrate what I'm talking about. + * + * xxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx + * 001ppppp pppppppp pppppppp pppppppp + * || + * |<------------ bp ------------->| + * |<------------ nw --------------->| + * + * The trick of taking whole words off of @X@ stops working when there are + * only @nw@ words left. Then we have to mask off the bottom bits of @w@ + * before continuing. + */ + +/* --- @gfreduce_create@ --- * + * + * Arguments: @gfreduce *r@ = structure to fill in + * @mp *x@ = a (hopefully sparse) polynomial + * + * Returns: --- + * + * Use: Initializes a context structure for reduction. + */ + +void gfreduce_create(gfreduce *r, mp *p) +{ + instr_v iv = DA_INIT; + unsigned long d, dw; + mpscan sc; + unsigned long i; + gfreduce_instr *ip; + unsigned f = 0; + size_t w, ww, wi, wl, ll; + + /* --- Sort out the easy stuff --- */ + + d = mp_bits(p); assert(d); d--; + r->lim = d/MPW_BITS; + dw = d%MPW_BITS; + if (!dw) + r->mask = 0; + else { + r->mask = MPW(((mpw)-1) << dw); + r->lim++; + } + r->p = mp_copy(p); + + /* --- Stash a new instruction --- */ + +#define INSTR(op_, arg_) do { \ + DA_ENSURE(&iv, 1); \ + DA(&iv)[DA_LEN(&iv)].op = (op_); \ + DA(&iv)[DA_LEN(&iv)].arg = (arg_); \ + DA_EXTEND(&iv, 1); \ +} while (0) + +#define f_lsr 1u + + w = (d + MPW_BITS - 1)/MPW_BITS; + INSTR(GFRI_LOAD, w); + wi = DA_LEN(&iv); + f = 0; + ll = 0; + for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) { + if (!mp_bit(&sc)) + continue; + ww = (d - i + MPW_BITS - 1)/MPW_BITS; + if (ww != w) { + wl = DA_LEN(&iv); + INSTR(GFRI_STORE, w); + if (!ll) + ll = DA_LEN(&iv); + if (!(f & f_lsr)) + INSTR(GFRI_LOAD, ww); + else { + INSTR(GFRI_LOAD, w - 1); + for (; wi < wl; wi++) { + ip = &DA(&iv)[wi]; + assert(ip->op == GFRI_LSL); + if (ip->arg) + INSTR(GFRI_LSR, MPW_BITS - ip->arg); + } + if (w - 1 != ww) { + INSTR(GFRI_STORE, w - 1); + INSTR(GFRI_LOAD, ww); + } + f &= ~f_lsr; + } + w = ww; + wi = DA_LEN(&iv); + } + INSTR(GFRI_LSL, (i - d)%MPW_BITS); + if ((i - d)%MPW_BITS) + f |= f_lsr; + } + wl = DA_LEN(&iv); + INSTR(GFRI_STORE, w); + if (!ll) + ll = DA_LEN(&iv); + if (f & f_lsr) { + INSTR(GFRI_LOAD, w - 1); + for (; wi < wl; wi++) { + ip = &DA(&iv)[wi]; + assert(ip->op == GFRI_LSL); + if (ip->arg) + INSTR(GFRI_LSR, MPW_BITS - ip->arg); + } + INSTR(GFRI_STORE, w - 1); + } + +#undef INSTR + + r->in = DA_LEN(&iv); + r->iv = xmalloc(r->in * sizeof(gfreduce_instr)); + r->liv = r->iv + ll; + memcpy(r->iv, DA(&iv), r->in * sizeof(gfreduce_instr)); + DA_DESTROY(&iv); +} + +/* --- @gfreduce_destroy@ --- * + * + * Arguments: @gfreduce *r@ = structure to free + * + * Returns: --- + * + * Use: Reclaims the resources from a reduction context. + */ + +void gfreduce_destroy(gfreduce *r) +{ + mp_drop(r->p); + xfree(r->iv); +} + +/* --- @gfreduce_dump@ --- * + * + * Arguments: @gfreduce *r@ = structure to dump + * @FILE *fp@ = file to dump on + * + * Returns: --- + * + * Use: Dumps a reduction context. + */ + +void gfreduce_dump(gfreduce *r, FILE *fp) +{ + size_t i; + + fprintf(fp, "poly = "); mp_writefile(r->p, fp, 16); + fprintf(fp, "\n lim = %lu; mask = %lx\n", + (unsigned long)r->lim, (unsigned long)r->mask); + for (i = 0; i < r->in; i++) { + static const char *opname[] = { "load", "lsl", "lsr", "store" }; + assert(r->iv[i].op < N(opname)); + fprintf(fp, " %s %lu\n", + opname[r->iv[i].op], + (unsigned long)r->iv[i].arg); + } +} + +/* --- @gfreduce_do@ --- * + * + * Arguments: @gfreduce *r@ = reduction context + * @mp *d@ = destination + * @mp *x@ = source + * + * Returns: Destination, @x@ reduced modulo the reduction poly. + */ + +static void run(const gfreduce_instr *i, const gfreduce_instr *il, + mpw *v, mpw z) +{ + mpw w = 0; + + for (; i < il; i++) { + switch (i->op) { + case GFRI_LOAD: w = *(v - i->arg); break; + case GFRI_LSL: w ^= z << i->arg; break; + case GFRI_LSR: w ^= z >> i->arg; break; + case GFRI_STORE: *(v - i->arg) = MPW(w); break; + default: abort(); + } + } +} + +mp *gfreduce_do(gfreduce *r, mp *d, mp *x) +{ + mpw *v, *vl; + const gfreduce_instr *il; + mpw z; + + /* --- Try to reuse the source's space --- */ + + MP_COPY(x); + if (d) MP_DROP(d); + MP_DEST(x, MP_LEN(x), x->f); + + /* --- Do the reduction --- */ + + il = r->iv + r->in; + if (MP_LEN(x) >= r->lim) { + v = x->v + r->lim; + vl = x->vl; + while (vl-- > v) { + while (*vl) { + z = *vl; + *vl = 0; + run(r->iv, il, vl, z); + } + } + if (r->mask) { + while (*vl & r->mask) { + z = *vl & r->mask; + *vl &= ~r->mask; + run(r->iv, il, vl, z); + } + } + } + + /* --- Done --- */ + + MP_SHRINK(x); + return (x); +} + +/* --- @gfreduce_sqrt@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: The square root of @x@ modulo @r->p@, or null. + */ + +mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x) +{ + mp *y = MP_COPY(x); + mp *z, *spare = MP_NEW; + unsigned long m = mp_bits(r->p) - 1; + unsigned long i; + + for (i = 0; i < m - 1; i++) { + mp *t = gf_sqr(spare, y); + spare = y; + y = gfreduce_do(r, t, t); + } + z = gf_sqr(spare, y); + z = gfreduce_do(r, z, z); + if (!MP_EQ(x, z)) { + mp_drop(y); + y = 0; + } + mp_drop(z); + mp_drop(d); + return (y); +} + +/* --- @gfreduce_trace@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *x@ = some polynomial + * + * Returns: The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$% + * if %$x \in \gf{2^m}$%). + */ + +int gfreduce_trace(gfreduce *r, mp *x) +{ + mp *y = MP_COPY(x); + mp *spare = MP_NEW; + unsigned long m = mp_bits(r->p) - 1; + unsigned long i; + int rc; + + for (i = 0; i < m - 1; i++) { + mp *t = gf_sqr(spare, y); + spare = y; + y = gfreduce_do(r, t, t); + y = gf_add(y, y, x); + } + rc = !MP_ISZERO(y); + mp_drop(spare); + mp_drop(y); + return (rc); +} + +/* --- @gfreduce_halftrace@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: The half-trace of @x@. + * (%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$% + * if %$x \in \gf{2^m}$% with %$m$% odd). + */ + +mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x) +{ + mp *y = MP_COPY(x); + mp *spare = MP_NEW; + unsigned long m = mp_bits(r->p) - 1; + unsigned long i; + + mp_drop(d); + for (i = 0; i < m - 1; i += 2) { + mp *t = gf_sqr(spare, y); + spare = y; + y = gfreduce_do(r, t, t); + t = gf_sqr(spare, y); + spare = y; + y = gfreduce_do(r, t, t); + y = gf_add(y, y, x); + } + mp_drop(spare); + return (y); +} + +/* --- @gfreduce_quadsolve@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: A polynomial @y@ such that %$y^2 + y = x$%, or null. + */ + +mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x) +{ + unsigned long m = mp_bits(r->p) - 1; + mp *t; + + MP_COPY(x); + if (m & 1) + d = gfreduce_halftrace(r, d, x); + else { + mp *z, *w, *rho = MP_NEW; + mp *spare = MP_NEW; + grand *fr = fibrand_create(0); + unsigned long i; + + for (;;) { + rho = mprand(rho, m, fr, 0); + z = MP_ZERO; + w = MP_COPY(rho); + for (i = 0; i < m - 1; i++) { + t = gf_sqr(spare, z); spare = z; z = gfreduce_do(r, t, t); + t = gf_sqr(spare, w); spare = w; w = gfreduce_do(r, t, t); + t = gf_mul(spare, w, x); t = gfreduce_do(r, t, t); spare = t; + z = gf_add(z, z, t); + w = gf_add(w, w, rho); + } + if (!MP_ISZERO(w)) + break; + MP_DROP(z); + MP_DROP(w); + } + if (d) MP_DROP(d); + MP_DROP(w); + MP_DROP(spare); + MP_DROP(rho); + fr->ops->destroy(fr); + d = z; + } + + t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d); + if (!MP_EQ(t, x)) { + MP_DROP(d); + d = 0; + } + MP_DROP(t); + MP_DROP(x); + d->v[0] &= ~(mpw)1; + return (d); +} + +/* --- @gfreduce_exp@ --- * + * + * Arguments: @gfreduce *gr@ = pointer to reduction context + * @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e \bmod m$%. + */ + +mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e) +{ + mp *x = MP_ONE; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + + MP_SHRINK(e); + if (!MP_LEN(e)) + ; + else if (MP_LEN(e) < EXP_THRESH) + EXP_SIMPLE(x, a, e); + else + EXP_WINDOW(x, a, e); + mp_drop(d); + mp_drop(spare); + return (x); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#define MP(x) mp_readstring(MP_NEW, #x, 0, 0) + +static int vreduce(dstr *v) +{ + mp *d = *(mp **)v[0].buf; + mp *n = *(mp **)v[1].buf; + mp *r = *(mp **)v[2].buf; + mp *c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, d); + c = gfreduce_do(&rr, MP_NEW, n); + if (!MP_EQ(c, r)) { + fprintf(stderr, "\n*** reduction failed\n*** "); + gfreduce_dump(&rr, stderr); + fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 16); + fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16); + fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int vmodexp(dstr *v) +{ + mp *p = *(mp **)v[0].buf; + mp *g = *(mp **)v[1].buf; + mp *x = *(mp **)v[2].buf; + mp *r = *(mp **)v[3].buf; + mp *c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, p); + c = gfreduce_exp(&rr, MP_NEW, g, x); + if (!MP_EQ(c, r)) { + fprintf(stderr, "\n*** modexp failed\n*** "); + fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16); + fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 16); + fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16); + fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16); + fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int vsqrt(dstr *v) +{ + mp *p = *(mp **)v[0].buf; + mp *x = *(mp **)v[1].buf; + mp *r = *(mp **)v[2].buf; + mp *c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, p); + c = gfreduce_sqrt(&rr, MP_NEW, x); + if (!MP_EQ(c, r)) { + fprintf(stderr, "\n*** sqrt failed\n*** "); + fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16); + fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16); + fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16); + fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int vtr(dstr *v) +{ + mp *p = *(mp **)v[0].buf; + mp *x = *(mp **)v[1].buf; + int r = *(int *)v[2].buf, c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, p); + c = gfreduce_trace(&rr, x); + if (c != r) { + fprintf(stderr, "\n*** trace failed\n*** "); + fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16); + fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16); + fprintf(stderr, "\n*** c = %d", c); + fprintf(stderr, "\n*** r = %d", r); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(p); mp_drop(x); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int vhftr(dstr *v) +{ + mp *p = *(mp **)v[0].buf; + mp *x = *(mp **)v[1].buf; + mp *r = *(mp **)v[2].buf; + mp *c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, p); + c = gfreduce_halftrace(&rr, MP_NEW, x); + if (!MP_EQ(c, r)) { + fprintf(stderr, "\n*** halftrace failed\n*** "); + fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16); + fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16); + fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16); + fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int vquad(dstr *v) +{ + mp *p = *(mp **)v[0].buf; + mp *x = *(mp **)v[1].buf; + mp *r = *(mp **)v[2].buf; + mp *c; + int ok = 1; + gfreduce rr; + + gfreduce_create(&rr, p); + c = gfreduce_quadsolve(&rr, MP_NEW, x); + if (!MP_EQ(c, r)) { + fprintf(stderr, "\n*** quadsolve failed\n*** "); + fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16); + fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16); + fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16); + fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16); + fprintf(stderr, "\n"); + ok = 0; + } + gfreduce_destroy(&rr); + mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk defs[] = { + { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } }, + { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { "sqrt", vsqrt, { &type_mp, &type_mp, &type_mp, 0 } }, + { "trace", vtr, { &type_mp, &type_mp, &type_int, 0 } }, + { "halftrace", vhftr, { &type_mp, &type_mp, &type_mp, 0 } }, + { "quadsolve", vquad, { &type_mp, &type_mp, &type_mp, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + test_run(argc, argv, defs, SRCDIR"/tests/gfreduce"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gfreduce.h b/gfreduce.h new file mode 100644 index 0000000..2fc4c0a --- /dev/null +++ b/gfreduce.h @@ -0,0 +1,186 @@ +/* -*-c-*- + * + * $Id: gfreduce.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Reduction modulo sparse binary polynomials + * + * (c) 2004 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. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gfreduce.h,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +#ifndef CATACOMB_GFREDUCE_H +#define CATACOMB_GFREDUCE_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef CATACOMB_GF_H +# include "gf.h" +#endif + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct gfreduce_instr { + unsigned op; /* Instruction opcode */ + size_t arg; /* Immediate argument */ +} gfreduce_instr; + +enum { + GFRI_LOAD, /* Load @p[arg]@ */ + GFRI_LSL, /* XOR with @w << arg@ */ + GFRI_LSR, /* XOR with @w >> arg@ */ + GFRI_STORE, /* Store @p[arg]@ */ + GFRI_MAX +}; + +typedef struct gfreduce { + size_t lim; /* Word of degree bit */ + mpw mask; /* Mask for degree word */ + mp *p; /* Copy of the polynomial */ + size_t in; /* Number of instruction words */ + gfreduce_instr *iv, *liv; /* Vector of instructions */ +} gfreduce; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @gfreduce_create@ --- * + * + * Arguments: @gfreduce *r@ = structure to fill in + * @mp *x@ = a (hopefully sparse) polynomial + * + * Returns: --- + * + * Use: Initializes a context structure for reduction. + */ + +extern void gfreduce_create(gfreduce */*r*/, mp */*p*/); + +/* --- @gfreduce_destroy@ --- * + * + * Arguments: @gfreduce *r@ = structure to free + * + * Returns: --- + * + * Use: Reclaims the resources from a reduction context. + */ + +extern void gfreduce_destroy(gfreduce */*r*/); + +/* --- @gfreduce_dump@ --- * + * + * Arguments: @gfreduce *r@ = structure to dump + * @FILE *fp@ = file to dump on + * + * Returns: --- + * + * Use: Dumps a reduction context. + */ + +extern void gfreduce_dump(gfreduce */*r*/, FILE */*fp*/); + +/* --- @gfreduce_do@ --- * + * + * Arguments: @gfreduce *r@ = reduction context + * @mp *d@ = destination + * @mp *x@ = source + * + * Returns: Destination, @x@ reduced modulo the reduction poly. + */ + +extern mp *gfreduce_do(gfreduce */*r*/, mp */*d*/, mp */*x*/); + +/* --- @gfreduce_sqrt@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: The square root of @x@ modulo @r->p@, or null. + */ + +extern mp *gfreduce_sqrt(gfreduce */*r*/, mp */*d*/, mp */*x*/); + +/* --- @gfreduce_trace@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *x@ = some polynomial + * + * Returns: The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$% + * if %$x \in \gf{2^m}$%). + */ + +extern int gfreduce_trace(gfreduce */*r*/, mp */*x*/); + +/* --- @gfreduce_halftrace@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: The half-trace of @x@. + * (%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$% + * if %$x \in \gf{2^m}$% with %$m$% odd). + */ + +extern mp *gfreduce_halftrace(gfreduce */*r*/, mp */*d*/, mp */*x*/); + +/* --- @gfreduce_quadsolve@ --- * + * + * Arguments: @gfreduce *r@ = pointer to reduction context + * @mp *d@ = destination + * @mp *x@ = some polynomial + * + * Returns: A polynomial @y@ such that %$y^2 + y = x$%, or null. + */ + +extern mp *gfreduce_quadsolve(gfreduce */*r*/, mp */*d*/, mp */*x*/); + +/* --- @gfreduce_exp@ --- * + * + * Arguments: @gfreduce *gr@ = pointer to reduction context + * @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent + * + * Returns: Result, %$a^e \bmod m$%. + */ + +extern mp *gfreduce_exp(gfreduce */*gr*/, mp */*d*/, mp */*a*/, mp */*e*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/gfx-sqr.c b/gfx-sqr.c index b253a32..778f85a 100644 --- a/gfx-sqr.c +++ b/gfx-sqr.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gfx-sqr.c,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * $Id: gfx-sqr.c,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $ * * Sqaring binary polynomials * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: gfx-sqr.c,v $ + * Revision 1.1.4.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.1 2000/10/08 15:49:37 mdw * First glimmerings of binary polynomial arithmetic. * @@ -38,7 +41,7 @@ /*----- Header files ------------------------------------------------------*/ #include "mpx.h" -/* #include "gfx.h" */ +#include "gfx.h" #include "gfx-sqr-tab.h" /*----- Static variables --------------------------------------------------*/ @@ -99,7 +102,7 @@ void gfx_sqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl) /* --- Output buffering --- */ - if (bb > MPW_BITS) { + if (bb >= MPW_BITS) { *dv++ = MPW(aa); if (dv >= dvl) return; diff --git a/gfx.h b/gfx.h index 778baee..d525650 100644 --- a/gfx.h +++ b/gfx.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: gfx.h,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * $Id: gfx.h,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $ * * Low-level arithmetic on binary polynomials * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: gfx.h,v $ + * Revision 1.1.4.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.1 2000/10/08 15:49:37 mdw * First glimmerings of binary polynomial arithmetic. * @@ -112,6 +115,19 @@ extern void gfx_mul(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, const mpw */*bv*/, const mpw */*bvl*/); +/* --- @gfx_sqr@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = argument vector base and limit + * + * Returns: --- + * + * Use: Performs squaring of binary polynomials. + */ + +extern void gfx_sqr(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/); + /* --- @gfx_div@ --- * * * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit diff --git a/mp-gcd.c b/mp-gcd.c index a17a502..f55d0aa 100644 --- a/mp-gcd.c +++ b/mp-gcd.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mp-gcd.c,v 1.5 2000/10/08 12:02:41 mdw Exp $ + * $Id: mp-gcd.c,v 1.5.4.1 2004/03/21 22:39:46 mdw Exp $ * * Extended GCD calculation * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: mp-gcd.c,v $ + * Revision 1.5.4.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * * Revision 1.5 2000/10/08 12:02:41 mdw * Use Euclid's algorithm rather than the binary one. * @@ -74,12 +77,10 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) mp *q = MP_NEW; unsigned f = 0; - enum { - f_swap = 1u, - f_aneg = 2u, - f_bneg = 4u, - f_ext = 8u - }; +#define f_swap 1u +#define f_aneg 2u +#define f_bneg 4u +#define f_ext 8u /* --- Sort out some initial flags --- */ diff --git a/tests/gf b/tests/gf new file mode 100644 index 0000000..0c3987f --- /dev/null +++ b/tests/gf @@ -0,0 +1,70 @@ +# $Id: gf,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ +# +# Test cases for higher-level binary poly arithmetic. + +add { + 0 0 0; + 1 1 0; + 1 2 3; + 4 5 1; + 0x7fb838a8a0a95046b9d9d9fb4440f7bb + 0xc1a7bd3b4e853fc92d4e1588719986aa + 0xbe1f8593ee2c6f8f9497cc7335d97111; + 0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b21 + 0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b22 3; +} + +mul { + 0 0 0; + 1 0 0; + 0 1 0; + 1 1 1; + 0x7fb838a8a0a95046b9d9d9fb4440f7bb + 0xc1a7bd3b4e853fc92d4e1588719986aa + 0x207ccad257b4ed64447158315bfb9aca5cbc5622cfb8fcbb1380eea1bc5c624e; + 0xc1a7bd3b4e853fc92d4e1588719986aa + 0x283ed59f1226dcefa7ff0ef87ceff5d5 + 0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b22; + 0xbe1f8593ee2c6f8f9497cc7335d97111 + 0x35a8e33503b3695be00528f8b82db931 + 0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b21; +} + +sqr { + 0 0; + 1 1; + 3 5; + 0x7fb838a8a0a95046b9d9d9fb4440f7bb + 0x1555454005404440440044411100101445415141514155451010100055154545; + 0x01f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3 + 0x0155004001541441551011510504111011050015141444454140511111554414010450154545041554440501455004140401514041104554415000450141144505; +} + +div { + 0 1 0 0; + 0x207ccad257b4ed64447158315bfb9aca5cbc5622cfb8fcbb1380eea1bc5c624e + 0x7fb838a8a0a95046b9d9d9fb4440f7bb + 0xc1a7bd3b4e853fc92d4e1588719986aa 0; + 0x6e0e2a282a5411ae76767ed1103deef069ef4ed3a14ff24b + 0x5385621c6661aaa35a24150d2c08332e + 0x01c2334cc957151dc7 + 0x398c4111da6d06cdf3d83704ee403101; +} + +gcd { + 0xc1a7bd3b4e853fc92d4e1588719986aa + 0xbe1f8593ee2c6f8f9497cc7335d97111 + 3 + 0x283ed59f1226dcefa7ff0ef87ceff5d5 + 0x35a8e33503b3695be00528f8b82db931; + 0xbe1f8593ee2c6f8f9497cc7335d97111 + 0xc1a7bd3b4e853fc92d4e1588719986aa + 3 + 0x35a8e33503b3695be00528f8b82db931 + 0x283ed59f1226dcefa7ff0ef87ceff5d5; + 0x800000000000000000000000000000000000000c9 + 0x3f0eba16286a2d57ea0991168d4994637e8343e36 + 1 + 0xa17e704470d80cb5a78f295db0ce543dda16a169 + 0x3c8c172e24598e90b9542e6b8f6571f54be572b50; +} diff --git a/tests/gfreduce b/tests/gfreduce new file mode 100644 index 0000000..aec5318 --- /dev/null +++ b/tests/gfreduce @@ -0,0 +1,61 @@ +# $Id: gfreduce,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ +# +# Test efficient polynomial reduction + +reduce { + 0x10000000 + 0x4509823098098435 + 0x8098435; + 0x100000000000000050002 + 0x4509823098098435 + 0x4509823098098435; + 0x100000000000000050002 + 0x450982309809843545609843098560803495 + 0x144f98a2f5cbc4773cfd; + 0xb2ca471b0867d5fae2e4f27a2d2706da + 0xf254423fef93d5d7a76ecf22c656c1352c53257875945d33 + 0x582f783fc210f72814780e69b0bd29ff; +} + +modexp { + 0x20000000000000000000000000000000000000000000000000000000000001001 + 0x02 + 0x1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff + 1; + 0x20000000000000000000000000000000000000000000000000000000000001001 + 0x435932098459080438094509845 + 0x1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff + 1; + 0x10000000000000000000000000000000000000000003 + 0x02 + 0x0fffffffffffffffffffffffffffffffffffffffffff + 1; + 0x10000000000000000000000000000000000000000003 + 0x34235950984598345900983409845690805680985 + 0x0fffffffffffffffffffffffffffffffffffffffffff + 1; +} + +sqrt { + 0x20000000000000000000000000000000000000000000000000000000000001001 + 0x1f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3 + 0x7fb838a8a0a95046b9d9d9fb4440f7bbc1a7bd3b4e853fc92d4e1588719986aa; + 0x10000000000000000000000000000000000000000003 + 0x4594094509835690805698083560980459903450984 + 0x820291881a244a02840a2f8ece3f23f88f38bf0b3a; +} + +halftrace { + 0x20000000000000000000000000000000000000000000000000000000000001001 + 0x174e65c7d14a8ec286df8c7df17662f13f1d3563f13c8c63f23f5d0bd5d1b45cd + 0x8d68905434b020ccb849e17a03a5c441d2a104aaf523699c1cc7a93174d21d9d; +} + +quadsolve { + 0x20000000000000000000000000000000000000000000000000000000000001001 + 0x174e65c7d14a8ec286df8c7df17662f13f1d3563f13c8c63f23f5d0bd5d1b45cd + 0x8d68905434b020ccb849e17a03a5c441d2a104aaf523699c1cc7a93174d21d9c; + 0x10000000000000000000000000000000000000000003 + 0x3b818b447e90713da04f13c3b07cb5e2681d08e4700 + 0x27aa17c97dfa80bbdef9f91b243c6e6ddba1a223cac; +} diff --git a/tests/gfx b/tests/gfx index 17bdbf7..866bd39 100644 --- a/tests/gfx +++ b/tests/gfx @@ -1,6 +1,6 @@ # Test vectors for low-level GF functions # -# $Id: gfx,v 1.1 2000/10/08 16:01:48 mdw Exp $ +# $Id: gfx,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $ # --- Addition (and subtraction) --- @@ -356,6 +356,11 @@ sqr { 1500551145554555115501404405515454551004050045401511450544004455; 93c9240145720297b96a404941792a21 4105504104100001101115040004411545411444100010411001154104440401; + + # --- Regression --- + + 01f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3 + 0155004001541441551011510504111011050015141444454140511111554414010450154545041554440501455004140401514041104554415000450141144505; } div { -- 2.11.0