-/**************************************************************
- *
- * 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);
-}