+ return (d);
+ } else {
+ ecctx *cc = (ecctx *)c;
+ dx = F_SQR(f, MP_NEW, a->x); /* %$x_0^2$% */
+ dx = F_TPL(f, dx, dx); /* %$3 x_0^2$% */
+ dx = F_ADD(f, dx, dx, cc->a); /* %$3 x_0^2 + A$% */
+ dy = F_DBL(f, MP_NEW, a->y); /* %$2 y_0$% */
+ dy = F_INV(f, dy, dy); /* %$(2 y_0)^{-1}$% */
+ lambda = F_MUL(f, MP_NEW, dx, dy);
+ /* %$\lambda = (3 x_0^2 + A)/(2 y_0)$% */
+ }
+
+ dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */
+ dx = F_SUB(f, dx, dx, a->x); /* %$\lambda^2 - x_0$% */
+ 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$% */
+
+ 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;
+ mp *p, *q, *r, *w, *u, *s, *dx, *dy, *dz;
+
+ q = F_SQR(f, MP_NEW, a->z); /* %$z_0^2$% */
+ u = F_MUL(f, MP_NEW, q, b->x); /* %$u = x_1 z_0^2$% */
+ p = F_MUL(f, MP_NEW, q, b->y); /* %$y_1 z_0^2$% */
+ s = F_MUL(f, q, p, a->z); /* %$s = y_1 z_0^3$% */
+
+ 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(r);
+ return (c->ops->dbl(c, d, a));
+ } else {
+ MP_DROP(r);
+ EC_SETINF(d);
+ return (d);
+ }
+ }
+ u = F_ADD(f, u, u, a->x); /* %$t = x_0 + u$% */
+ s = F_ADD(f, s, s, a->y); /* %$m = y_0 + r$% */
+
+ dz = F_MUL(f, MP_NEW, a->z, w); /* %$z' = z_0 w$% */
+
+ p = F_SQR(f, MP_NEW, w); /* %$w^2$% */
+ q = F_MUL(f, MP_NEW, p, u); /* %$t w^2$% */
+ u = F_MUL(f, u, p, w); /* %$w^3$% */
+ p = F_MUL(f, p, u, s); /* %$m w^3$% */
+
+ dx = F_SQR(f, u, r); /* %$r^2$% */
+ dx = F_SUB(f, dx, dx, q); /* %$x' = r^2 - t w^2$% */
+
+ s = F_DBL(f, s, dx); /* %$2 x'$% */
+ q = F_SUB(f, q, q, s); /* %$v = t w^2 - 2 x'$% */
+ dy = F_MUL(f, s, q, r); /* %$v r$% */
+ dy = F_SUB(f, dy, dy, p); /* %$v r - m w^3$% */
+ dy = F_HLV(f, dy, dy); /* %$y' = (v r - m w^3)/2$% */
+
+ EC_DESTROY(d);
+ d->x = dx;
+ d->y = dy;
+ d->z = dz;
+ MP_DROP(p);
+ MP_DROP(q);
+ MP_DROP(r);
+ MP_DROP(w);
+ }
+ return (d);
+}
+
+static int eccheck(ec_curve *c, const ec *p)
+{
+ ecctx *cc = (ecctx *)c;
+ field *f = c->f;
+ int rc;
+ mp *l = F_SQR(f, MP_NEW, p->y);
+ mp *x = F_SQR(f, MP_NEW, p->x);
+ mp *r = F_MUL(f, MP_NEW, x, p->x);
+ x = F_MUL(f, x, cc->a, p->x);
+ r = F_ADD(f, r, r, x);
+ r = F_ADD(f, r, r, cc->b);
+ rc = MP_EQ(l, r) ? 0 : -1;
+ mp_drop(l);
+ mp_drop(x);
+ mp_drop(r);
+ 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);
+ DESTROY(cc);
+}
+
+/* --- @ec_prime@, @ec_primeproj@ --- *
+ *
+ * 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 prime field. The @primeproj@ variant uses projective
+ * coordinates, which can be a win.
+ */
+
+extern ec_curve *ec_prime(field *f, mp *a, mp *b)
+{
+ ecctx *cc = CREATE(ecctx);
+ cc->c.ops = &ec_primeops;
+ cc->c.f = f;
+ cc->a = F_IN(f, MP_NEW, a);
+ cc->b = F_IN(f, MP_NEW, b);
+ return (&cc->c);
+}
+
+extern ec_curve *ec_primeproj(field *f, mp *a, mp *b)
+{
+ ecctx *cc = CREATE(ecctx);
+ mp *ax;
+
+ ax = mp_add(MP_NEW, a, MP_THREE);
+ ax = F_IN(f, ax, ax);
+ if (F_ZEROP(f, ax))
+ cc->c.ops = &ec_primeprojxops;
+ else
+ cc->c.ops = &ec_primeprojops;
+ MP_DROP(ax);
+ cc->c.f = f;
+ cc->a = F_IN(f, MP_NEW, a);
+ cc->b = F_IN(f, MP_NEW, b);
+ return (&cc->c);
+}
+
+static const ec_ops ec_primeops = {
+ ecdestroy, ec_idin, ec_idout, ec_idfix,
+ 0, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
+};
+
+static const ec_ops ec_primeprojops = {
+ ecdestroy, ec_projin, ec_projout, ec_projfix,
+ 0, ecneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck
+};
+
+static const ec_ops ec_primeprojxops = {
+ ecdestroy, ec_projin, ec_projout, ec_projfix,
+ 0, ecneg, ecprojadd, ec_stdsub, ecprojxdbl, ecprojcheck
+};
+
+/*----- Test rig ----------------------------------------------------------*/
+
+#ifdef TEST_RIG