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
## -*-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
##
##----- 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.
##
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 \
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 \
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)
--- /dev/null
+/* -*-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 -------------------------------------------------*/
+
/* -*-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
*
/*----- 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.
*
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);
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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <mLib/sub.h>
+
+#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 -------------------------------------------------*/
/* -*-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
*
/*----- 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
*
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);
}
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;
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);
}
#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);
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);
/* -*-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
*
/*----- 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
*
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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <mLib/sub.h>
+
+#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 -------------------------------------------------*/
/* -*-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
*
/*----- 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
*
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;
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@ --- *
/* -*-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
*
/*----- 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
*
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 --- */
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;
#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 ------------------------------------------*/
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
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
--- /dev/null
+/* -*-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
--- /dev/null
+/* -*-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 <mLib/alloc.h>
+#include <mLib/darray.h>
+#include <mLib/macros.h>
+
+#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
+ * |<rp>|
+ * |<------------ 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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
/* -*-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
*
/*----- 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.
*
/*----- Header files ------------------------------------------------------*/
#include "mpx.h"
-/* #include "gfx.h" */
+#include "gfx.h"
#include "gfx-sqr-tab.h"
/*----- Static variables --------------------------------------------------*/
/* --- Output buffering --- */
- if (bb > MPW_BITS) {
+ if (bb >= MPW_BITS) {
*dv++ = MPW(aa);
if (dv >= dvl)
return;
/* -*-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
*
/*----- 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.
*
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
/* -*-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
*
/*----- 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.
*
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 --- */
--- /dev/null
+# $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;
+}
--- /dev/null
+# $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;
+}
# 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) ---
1500551145554555115501404405515454551004050045401511450544004455;
93c9240145720297b96a404941792a21
4105504104100001101115040004411545411444100010411001154104440401;
+
+ # --- Regression ---
+
+ 01f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3
+ 0155004001541441551011510504111011050015141444454140511111554414010450154545041554440501455004140401514041104554415000450141144505;
}
div {
projects / u / mdw / catacomb / commitdiff