+/**************************************************************
+ *
+ * ellproj.c
+ *
+ Fast algorithms for fundamental elliptic curve arithmetic,
+ projective format. Such algorithms apply in domains such as:
+ -- factoring
+ -- primality studies (e.g. rigorous primality proofs)
+ -- elliptic curve cryptography (ECC)
+
+ PROJECTIVE FORMAT
+
+ Functions are supplied herein for projective format
+ of points. Alternative formats differ in their
+ range of applicability, efficiency, and so on.
+ Primary advantages of the projective format herein are:
+ -- No explicit inversions (until perhaps one such at the end of
+ an elliptic multiply operation)
+ -- Fairly low operation count (~11 muls for point doubling,
+ ~16 muls for point addition)
+
+ The elliptic curve is over F_p, with p > 3 prime, and Weierstrass
+ parameterization:
+
+ y^2 = x^3 + a x + b
+
+ The projective-format coordinates are actually stored in
+ the form {X, Y, Z}, with true x,y
+ coordinates on the curve given by {x,y} = {X/Z^2, Y/Z^3}.
+ The function normalize_proj() performs the
+ transformation from projective->true.
+ (The other direction is trivial, i.e. {x,y} -> {x,y,1} will do.)
+ The basic point multiplication function is
+
+ ell_mul_proj()
+
+ which obtains the result k * P for given point P and integer
+ multiplier k. If true {x,y} are required for a multiple, one
+ passes a point P = {X, Y, 1} to ell_mul_proj(), then afterwards
+ calls normalize_proj(),
+
+ Projective format is an answer to the typical sluggishness of
+ standard elliptic arithmetic, whose explicit inversion in the
+ field is, depending of course on the machinery and programmer,
+ costly. Projective format is thus especially interesting for
+ cryptography.
+
+ REFERENCES
+
+ Crandall R and Pomerance C 1998, "Prime numbers: a computational
+ perspective," Springer-Verlag, manuscript
+
+ Solinas J 1998, IEEE P1363 Annex A (draft standard)
+
+ LEGAL AND PATENT NOTICE
+
+ This and related PSI library source code is intended solely for
+ educational and research applications, and should not be used
+ for commercial purposes without explicit permission from PSI
+ (not to mention proper clearance of legal restrictions beyond
+ the purview of PSI).
+ The code herein will not perform cryptography per se,
+ although the number-theoretical algorithms herein -- all of which
+ are in the public domain -- can be used in principle to effect
+ what is known as elliptic curve cryptography (ECC). Note that
+ there are strict constraints on how cryptography may be used,
+ especially in regard to exportability.
+ Therefore one should avoid any casual distribution of actual
+ cryptographic software. Along similar lines, because of various
+ patents, proprietary to Apple Computer, Inc., and perhaps to other
+ organizations, one should not tacitly assume that an ECC scheme is
+ unconstrained. For example,the commercial use of certain fields
+ F_p^k (i.e., fixation of certain primes p) is covered in Apple
+ patents.
+
+ * Updates:
+ * 3 Apr 98 REC Creation
+ *
+ * c. 1998 Perfectly Scientific, Inc.
+ * All Rights Reserved.
+ *
+ *
+ *************************************************************/
+
+/* include files */
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#ifdef _WIN32
+
+#include <process.h>
+
+#endif
+
+#include <string.h>
+#include "giants.h"
+#include "ellproj.h"
+#include "tools.h"
+
+/* global variables */
+
+static giant t0 = NULL, t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL,
+ t5 = NULL, t6 = NULL, t7 = NULL;
+
+/**************************************************************
+ *
+ * Maintenance functions
+ *
+ **************************************************************/
+
+point_proj
+new_point_proj(int shorts)
+{
+ point_proj pt;
+
+ if(t0 == NULL) init_ell_proj(shorts);
+ pt = (point_proj) malloc(sizeof(point_struct_proj));
+ pt->x = newgiant(shorts);
+ pt->y = newgiant(shorts);
+ pt->z = newgiant(shorts);
+ return(pt);
+}
+
+void
+free_point_proj(point_proj pt)
+{
+ free(pt->x); free(pt->y); free(pt->z);
+ free(pt);
+}
+
+void
+ptop_proj(point_proj pt1, point_proj pt2)
+{
+ gtog(pt1->x, pt2->x);
+ gtog(pt1->y, pt2->y);
+ gtog(pt1->z, pt2->z);
+}
+
+void
+init_ell_proj(int shorts)
+/* Called by new_point_proj(), to set up giant registers. */
+{
+ t0 = newgiant(shorts);
+ t1 = newgiant(shorts);
+ t2 = newgiant(shorts);
+ t3 = newgiant(shorts);
+ t4 = newgiant(shorts);
+ t5 = newgiant(shorts);
+ t6 = newgiant(shorts);
+ t7 = newgiant(shorts);
+}
+
+
+/**************************************************************
+ *
+ * Elliptic curve operations
+ *
+ **************************************************************/
+
+/* Begin projective-format functions for
+
+ y^2 = x^3 + a x + b.
+
+ These are useful in elliptic curve cryptography (ECC).
+ A point is kept as a triple {X,Y,Z}, with the true (x,y)
+ coordinates given by
+
+ {x,y} = {X/Z^2, Y/Z^3}
+
+ The function normalize_proj() performs the inverse conversion to get
+ the true (x,y) pair.
+ */
+
+void
+ell_double_proj(point_proj pt, giant a, giant p)
+/* pt := 2 pt on the curve. */
+{
+ giant x = pt->x, y = pt->y, z = pt->z;
+
+ if(isZero(y) || isZero(z)) {
+ itog(1,x); itog(1,y); itog(0,z);
+ return;
+ }
+ gtog(z,t1); squareg(t1); modg(p, t1);
+ squareg(t1); modg(p, t1);
+ mulg(a, t1); modg(p, t1); /* t1 := a z^4. */
+ gtog(x, t2); squareg(t2); smulg(3, t2); modg(p, t2); /* t2 := 3x^2. */
+ addg(t2, t1); modg(p, t1); /* t1 := slope m. */
+ mulg(y, z); addg(z,z); modg(p, z); /* z := 2 y z. */
+ gtog(y, t2); squareg(t2); modg(p, t2); /* t2 := y^2. */
+ gtog(t2, t3); squareg(t3); modg(p, t3); /* t3 := y^4. */
+ gshiftleft(3, t3); /* t3 := 8 y^4. */
+ mulg(x, t2); gshiftleft(2, t2); modg(p, t2); /* t2 := 4xy^2. */
+ gtog(t1, x); squareg(x); modg(p, x);
+ subg(t2, x); subg(t2, x); modg(p, x); /* x done. */
+ gtog(t1, y); subg(x, t2); mulg(t2, y); subg(t3, y);
+ modg(p, y);
+}
+/*
+elldouble[pt_] := Block[{x,y,z,m,y2,s},
+ x = pt[[1]]; y = pt[[2]]; z = pt[[3]];
+ If[(y==0) || (z==0), Return[{1,1,0}]];
+ m = Mod[3 x^2 + a Mod[Mod[z^2,p]^2,p],p];
+ z = Mod[2 y z, p];
+ y2 = Mod[y^2,p];
+ s = Mod[4 x y2,p];
+ x = Mod[m^2 - 2s,p];
+ y = Mod[m(s - x) - 8 y2^2,p];
+ Return[{x,y,z}];
+];
+*/
+
+void
+ell_add_proj(point_proj pt0, point_proj pt1, giant a, giant p)
+/* pt0 := pt0 + pt1 on the curve. */
+{
+ giant x0 = pt0->x, y0 = pt0->y, z0 = pt0->z,
+ x1 = pt1->x, y1 = pt1->y, z1 = pt1->z;
+
+ if(isZero(z0)) {
+ gtog(x1,x0); gtog(y1,y0); gtog(z1,z0);
+ return;
+ }
+ if(isZero(z1)) return;
+ gtog(x0, t1); gtog(y0,t2); gtog(z0, t3);
+ gtog(x1, t4); gtog(y1, t5);
+ if(!isone(z1)) {
+ gtog(z1, t6);
+ gtog(t6, t7); squareg(t7); modg(p, t7);
+ mulg(t7, t1); modg(p, t1);
+ mulg(t6, t7); modg(p, t7);
+ mulg(t7, t2); modg(p, t2);
+ }
+ gtog(t3, t7); squareg(t7); modg(p, t7);
+ mulg(t7, t4); modg(p, t4);
+ mulg(t3, t7); modg(p, t7);
+ mulg(t7, t5); modg(p, t5);
+ negg(t4); addg(t1, t4); modg(p, t4);
+ negg(t5); addg(t2, t5); modg(p, t5);
+ if(isZero(t4)) {
+ if(isZero(t5)) {
+ ell_double_proj(pt0, a, p);
+ } else {
+ itog(1, x0); itog(1, y0); itog(0, z0);
+ }
+ return;
+ }
+ addg(t1, t1); subg(t4, t1); modg(p, t1);
+ addg(t2, t2); subg(t5, t2); modg(p, t2);
+ if(!isone(z1)) {
+ mulg(t6, t3); modg(p, t3);
+ }
+ mulg(t4, t3); modg(p, t3);
+ gtog(t4, t7); squareg(t7); modg(p, t7);
+ mulg(t7, t4); modg(p, t4);
+ mulg(t1, t7); modg(p, t7);
+ gtog(t5, t1); squareg(t1); modg(p, t1);
+ subg(t7, t1); modg(p, t1);
+ subg(t1, t7); subg(t1, t7); modg(p, t7);
+ mulg(t7, t5); modg(p, t5);
+ mulg(t2, t4); modg(p, t4);
+ gtog(t5, t2); subg(t4,t2); modg(p, t2);
+ if(t2->n[0] & 1) { /* Test if t2 is odd. */
+ addg(p, t2);
+ }
+ gshiftright(1, t2);
+ gtog(t1, x0); gtog(t2, y0); gtog(t3, z0);
+}
+
+/*
+elladd[pt0_, pt1_] := Block[
+ {x0,y0,z0,x1,y1,z1,
+ t1,t2,t3,t4,t5,t6,t7},
+ x0 = pt0[[1]]; y0 = pt0[[2]]; z0 = pt0[[3]];
+ x1 = pt1[[1]]; y1 = pt1[[2]]; z1 = pt1[[3]];
+ If[z0 == 0, Return[pt1]];
+ If[z1 == 0, Return[pt0]];
+
+ t1 = x0;
+ t2 = y0;
+ t3 = z0;
+ t4 = x1;
+ t5 = y1;
+ If[(z1 != 1),
+ t6 = z1;
+ t7 = Mod[t6^2, p];
+ t1 = Mod[t1 t7, p];
+ t7 = Mod[t6 t7, p];
+ t2 = Mod[t2 t7, p];
+ ];
+ t7 = Mod[t3^2, p];
+ t4 = Mod[t4 t7, p];
+ t7 = Mod[t3 t7, p];
+ t5 = Mod[t5 t7, p];
+ t4 = Mod[t1-t4, p];
+ t5 = Mod[t2 - t5, p];
+ If[t4 == 0, If[t5 == 0,
+ Return[elldouble[pt0]],
+ Return[{1,1,0}]
+ ]
+ ];
+ t1 = Mod[2t1 - t4,p];
+ t2 = Mod[2t2 - t5, p];
+ If[z1 != 1, t3 = Mod[t3 t6, p]];
+ t3 = Mod[t3 t4, p];
+ t7 = Mod[t4^2, p];
+ t4 = Mod[t4 t7, p];
+ t7 = Mod[t1 t7, p];
+ t1 = Mod[t5^2, p];
+ t1 = Mod[t1-t7, p];
+ t7 = Mod[t7 - 2t1, p];
+ t5 = Mod[t5 t7, p];
+ t4 = Mod[t2 t4, p];
+ t2 = Mod[t5-t4, p];
+ If[EvenQ[t2], t2 = t2/2, t2 = (p+t2)/2];
+ Return[{t1, t2, t3}];
+];
+*/
+
+void
+ell_neg_proj(point_proj pt, giant p)
+/* pt := -pt on the curve. */
+{
+ negg(pt->y); modg(p, pt->y);
+}
+
+void
+ell_sub_proj(point_proj pt0, point_proj pt1, giant a, giant p)
+/* pt0 := pt0 - pt1 on the curve. */
+{
+ ell_neg_proj(pt1, p);
+ ell_add_proj(pt0, pt1,a,p);
+ ell_neg_proj(pt1,p);
+}
+
+void
+ell_mul_proj(point_proj pt0, point_proj pt1, giant k, giant a, giant p)
+/* General elliptic multiplication;
+ pt1 := k*pt0 on the curve,
+ with k an arbitrary integer.
+ */
+{
+ giant x = pt0->x, y = pt0->y, z = pt0->z,
+ xx = pt1->x, yy = pt1->y, zz = pt1->z;
+ int ksign, hlen, klen, b, hb, kb;
+
+ if(isZero(k)) {
+ itog(1, xx);
+ itog(1, yy);
+ itog(0, zz);
+ return;
+ }
+ ksign = k->sign;
+ if(ksign < 0) negg(k);
+ gtog(x,xx); gtog(y,yy); gtog(z,zz);
+ gtog(k, t0); addg(t0, t0); addg(k, t0); /* t0 := 3k. */
+ hlen = bitlen(t0);
+ klen = bitlen(k);
+ for(b = hlen-2; b > 0; b--) {
+ ell_double_proj(pt1,a,p);
+ hb = bitval(t0, b);
+ if(b < klen) kb = bitval(k, b); else kb = 0;
+ if((hb != 0) && (kb == 0))
+ ell_add_proj(pt1, pt0, a, p);
+ else if((hb == 0) && (kb !=0))
+ ell_sub_proj(pt1, pt0, a, p);
+ }
+ if(ksign < 0) {
+ ell_neg_proj(pt1, p);
+ k->sign = -k->sign;
+ }
+}
+
+/*
+elliptic[pt_, k_] := Block[{pt2, hh, kk, hb, kb, lenh, lenk},
+ If[k==0, Return[{1,1,0}]];
+ hh = Reverse[bitList[3k]];
+ kk = Reverse[bitList[k]];
+ pt2 = pt;
+ lenh = Length[hh];
+ lenk = Length[kk];
+ Do[
+ pt2 = elldouble[pt2];
+ hb = hh[[b]];
+ If[b <= lenk, kb = kk[[b]], kb = 0];
+ If[{hb,kb} == {1,0},
+ pt2 = elladd[pt2, pt],
+ If[{hb, kb} == {0,1},
+ pt2 = ellsub[pt2, pt]]
+ ]
+ ,{b, lenh-1, 2,-1}
+ ];
+ Return[pt2];
+];
+*/
+
+void
+normalize_proj(point_proj pt, giant p)
+/* Obtain actual x,y coords via normalization:
+ {x,y,z} := {x/z^2, y/z^3, 1}.
+ */
+
+{ giant x = pt->x, y = pt->y, z = pt->z;
+
+ if(isZero(z)) {
+ itog(1,x); itog(1,y);
+ return;
+ }
+ binvaux(p, z); gtog(z, t1);
+ squareg(z); modg(p, z);
+ mulg(z, x); modg(p, x);
+ mulg(t1, z); mulg(z, y); modg(p, y);
+ itog(1, z);
+}
+
+/*
+normalize[pt_] := Block[{z,z2,z3},
+ If[pt[[3]] == 0, Return[pt]];
+ z = ellinv[pt[[3]]];
+ z2 = Mod[z^2,p];
+ z3 = Mod[z z2,p];
+ Return[{Mod[pt[[1]] z2, p], Mod[pt[[2]] z3, p], 1}];
+ ];
+*/
+
+
+void
+find_point_proj(point_proj pt, giant seed, giant a, giant b, giant p)
+/* Starting with seed, finds a random (projective) point {x,y,1} on curve.
+ */
+{ giant x = pt->x, y = pt->y, z = pt->z;
+
+ modg(p, seed);
+ while(1) {
+ gtog(seed, x);
+ squareg(x); modg(p, x);
+ addg(a, x);
+ mulg(seed,x); addg(b, x);
+ modg(p, x); /* x := seed^3 + a seed + b. */
+ if(sqrtmod(p, x)) break; /* Test if cubic form has root. */
+ iaddg(1, seed);
+ }
+ gtog(x, y);
+ gtog(seed,x);
+ itog(1, z);
+}