X-Git-Url: https://git.saurik.com/apple/security.git/blobdiff_plain/5dd5f9ec28f304ca377c42fd7f711d6cf12b90e1..5c19dc3ae3bd8e40a9c028b0deddd50ff337692c:/OSX/libsecurity_cryptkit/lib/ellipticProj.c

diff --git a/OSX/libsecurity_cryptkit/lib/ellipticProj.c b/OSX/libsecurity_cryptkit/lib/ellipticProj.c
new file mode 100644
index 00000000..bdb24316
--- /dev/null
+++ b/OSX/libsecurity_cryptkit/lib/ellipticProj.c
@@ -0,0 +1,565 @@
+/* Copyright (c) 1998,2011,2014 Apple Inc.  All Rights Reserved.
+ *
+ * NOTICE: USE OF THE MATERIALS ACCOMPANYING THIS NOTICE IS SUBJECT
+ * TO THE TERMS OF THE SIGNED "FAST ELLIPTIC ENCRYPTION (FEE) REFERENCE
+ * SOURCE CODE EVALUATION AGREEMENT" BETWEEN APPLE, INC. AND THE
+ * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE,
+ * INC.  ANY USE OF THESE MATERIALS NOT PERMITTED BY SUCH AGREEMENT WILL
+ * EXPOSE YOU TO LIABILITY.
+ ***************************************************************************
+ *
+ * ellipticProj.c - elliptic projective algebra routines.
+ *
+ * Revision History
+ * ----------------
+ * 1 Sep 1998 at Apple
+ *	Created.
+ *
+ **************************************************************
+
+   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 normalizeProj() 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
+
+      ellMulProj()
+
+   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 ellMulProj(), then afterwards
+   calls normalizeProj(),
+
+   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
+
+		perspective," Springer-Verlag, manuscript
+
+   Solinas J 1998, IEEE P1363 Annex A (draft standard)
+
+***********************************************************/
+
+#include "ckconfig.h"
+#if CRYPTKIT_ELL_PROJ_ENABLE
+
+#include "ellipticProj.h"
+#include "falloc.h"
+#include "curveParams.h"
+#include "elliptic.h"
+#include "feeDebug.h"
+
+/*
+ * convert REC-style smulg to generic imulg
+ */
+#define smulg(s, g) imulg((unsigned)s, g)
+
+pointProj newPointProj(unsigned numDigits)
+{
+	pointProj pt;
+
+	pt = (pointProj) fmalloc(sizeof(pointProjStruct));
+	pt->x = newGiant(numDigits);
+	pt->y = newGiant(numDigits);
+	pt->z = newGiant(numDigits);
+	return(pt);
+}
+
+void freePointProj(pointProj pt)
+{
+	clearGiant(pt->x); freeGiant(pt->x);
+	clearGiant(pt->y); freeGiant(pt->y);
+	clearGiant(pt->z); freeGiant(pt->z);
+	ffree(pt);
+}
+
+void ptopProj(pointProj pt1, pointProj pt2)
+{
+	gtog(pt1->x, pt2->x);
+	gtog(pt1->y, pt2->y);
+	gtog(pt1->z, pt2->z);
+}
+
+/**************************************************************
+ *
+ *	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 normalizeProj() performs the inverse conversion to get
+   the true (x,y) pair.
+ */
+
+void ellDoubleProj(pointProj pt, curveParams *cp)
+/* pt := 2 pt on the curve. */
+{
+	giant x = pt->x, y = pt->y, z = pt->z;
+	giant t1;
+	giant t2;
+	giant t3;
+
+	if(isZero(y) || isZero(z)) {
+		int_to_giant(1,x); int_to_giant(1,y); int_to_giant(0,z);
+		return;
+	}
+	t1 = borrowGiant(cp->maxDigits);
+	t2 = borrowGiant(cp->maxDigits);
+	t3 = borrowGiant(cp->maxDigits);
+
+	if((cp->a->sign >= 0) || (cp->a->n[0] != 3)) { /* Path prior to Apr2001. */
+		gtog(z,t1); gsquare(t1); feemod(cp, t1);
+		gsquare(t1); feemod(cp, t1);
+		mulg(cp->a, t1); feemod(cp, t1);			/* t1 := a z^4. */
+		gtog(x, t2); gsquare(t2); feemod(cp, t2);
+			smulg(3, t2);                           /* t2 := 3x^2. */
+		addg(t2, t1); feemod(cp, t1);	 			/* t1 := slope m. */
+	} else { /* New optimization for a = -3 (post Apr 2001). */
+		gtog(z, t1); gsquare(t1); feemod(cp, t1);   /* t1 := z^2. */
+		gtog(x, t2); subg(t1, t2);                  /* t2 := x-z^2. */
+		addg(x, t1); smulg(3, t1);                  /* t1 := 3(x+z^2). */
+		mulg(t2, t1); feemod(cp, t1);               /* t1 := slope m. */
+	}
+	mulg(y, z); addg(z,z); feemod(cp, z);	  	/* z := 2 y z. */
+	gtog(y, t2); gsquare(t2); feemod(cp, t2); 	/* t2 := y^2. */
+	gtog(t2, t3); gsquare(t3); feemod(cp, t3);	/* t3 := y^4. */
+	gshiftleft(3, t3);  				/* t3 := 8 y^4. */
+	mulg(x, t2); gshiftleft(2, t2); feemod(cp, t2);	/* t2 := 4xy^2. */
+	gtog(t1, x); gsquare(x); feemod(cp, x);
+	subg(t2, x); subg(t2, x); feemod(cp, x);	/* x done. */
+	gtog(t1, y); subg(x, t2); mulg(t2, y); subg(t3, y);
+	feemod(cp, y);
+	returnGiant(t1);
+	returnGiant(t2);
+	returnGiant(t3);
+}
+
+void ellAddProj(pointProj pt0, pointProj pt1, curveParams *cp)
+/* 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;
+	giant t1;
+	giant t2;
+	giant t3;
+	giant t4;
+	giant t5;
+	giant t6;
+	giant t7;
+
+	if(isZero(z0)) {
+		gtog(x1,x0); gtog(y1,y0); gtog(z1,z0);
+		return;
+	}
+	if(isZero(z1)) return;
+
+	t1 = borrowGiant(cp->maxDigits);
+	t2 = borrowGiant(cp->maxDigits);
+	t3 = borrowGiant(cp->maxDigits);
+	t4 = borrowGiant(cp->maxDigits);
+	t5 = borrowGiant(cp->maxDigits);
+	t6 = borrowGiant(cp->maxDigits);
+	t7 = borrowGiant(cp->maxDigits);
+
+	gtog(x0, t1); gtog(y0,t2); gtog(z0, t3);
+	gtog(x1, t4); gtog(y1, t5);
+	if(!isone(z1)) {
+		gtog(z1, t6);
+		gtog(t6, t7); gsquare(t7); feemod(cp, t7);
+		mulg(t7, t1); feemod(cp, t1);
+		mulg(t6, t7); feemod(cp, t7);
+		mulg(t7, t2); feemod(cp, t2);
+	}
+	gtog(t3, t7); gsquare(t7); feemod(cp, t7);
+	mulg(t7, t4); feemod(cp, t4);
+	mulg(t3, t7); feemod(cp, t7);
+	mulg(t7, t5); feemod(cp, t5);
+	negg(t4); addg(t1, t4); feemod(cp, t4);
+	negg(t5); addg(t2, t5); feemod(cp, t5);
+	if(isZero(t4)) {
+		if(isZero(t5)) {
+			ellDoubleProj(pt0, cp);
+	    	} else {
+			int_to_giant(1, x0); int_to_giant(1, y0);
+			int_to_giant(0, z0);
+		}
+		goto out;
+	}
+	addg(t1, t1); subg(t4, t1); feemod(cp, t1);
+	addg(t2, t2); subg(t5, t2); feemod(cp, t2);
+	if(!isone(z1)) {
+		mulg(t6, t3); feemod(cp, t3);
+	}
+	mulg(t4, t3); feemod(cp, t3);
+	gtog(t4, t7); gsquare(t7); feemod(cp, t7);
+	mulg(t7, t4); feemod(cp, t4);
+	mulg(t1, t7); feemod(cp, t7);
+	gtog(t5, t1); gsquare(t1); feemod(cp, t1);
+	subg(t7, t1); feemod(cp, t1);
+	subg(t1, t7); subg(t1, t7); feemod(cp, t7);
+	mulg(t7, t5); feemod(cp, t5);
+	mulg(t2, t4); feemod(cp, t4);
+	gtog(t5, t2); subg(t4,t2); feemod(cp, t2);
+	if(t2->n[0] & 1) { /* Test if t2 is odd. */
+		addg(cp->basePrime, t2);
+	}
+	gshiftright(1, t2);
+	gtog(t1, x0); gtog(t2, y0); gtog(t3, z0);
+out:
+	returnGiant(t1);
+	returnGiant(t2);
+	returnGiant(t3);
+	returnGiant(t4);
+	returnGiant(t5);
+	returnGiant(t6);
+	returnGiant(t7);
+}
+
+
+void ellNegProj(pointProj pt, curveParams *cp)
+/* pt := -pt on the curve. */
+{
+	negg(pt->y); feemod(cp, pt->y);
+}
+
+void ellSubProj(pointProj pt0, pointProj pt1, curveParams *cp)
+/* pt0 := pt0 - pt1 on the curve. */
+{
+	ellNegProj(pt1, cp);
+	ellAddProj(pt0, pt1,cp);
+	ellNegProj(pt1, cp);
+}
+
+/*
+ * Simple projective multiply.
+ *
+ * pt := pt * k, result normalized.
+ */
+void ellMulProjSimple(pointProj pt0, giant k, curveParams *cp)
+{
+	pointProjStruct pt1;	// local, giants borrowed
+
+	CKASSERT(isone(pt0->z));
+	CKASSERT(cp->curveType == FCT_Weierstrass);
+
+	/* ellMulProj assumes constant pt0, can't pass as src and dst */
+	pt1.x = borrowGiant(cp->maxDigits);
+	pt1.y = borrowGiant(cp->maxDigits);
+	pt1.z = borrowGiant(cp->maxDigits);
+	ellMulProj(pt0, &pt1, k, cp);
+	normalizeProj(&pt1, cp);
+	CKASSERT(isone(pt1.z));
+
+	ptopProj(&pt1, pt0);
+	returnGiant(pt1.x);
+	returnGiant(pt1.y);
+	returnGiant(pt1.z);
+}
+
+void ellMulProj(pointProj pt0, pointProj pt1, giant k, curveParams *cp)
+/* 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;
+    	giant t0;
+
+	CKASSERT(cp->curveType == FCT_Weierstrass);
+	if(isZero(k)) {
+		int_to_giant(1, xx);
+		int_to_giant(1, yy);
+		int_to_giant(0, zz);
+		return;
+	}
+	t0 = borrowGiant(cp->maxDigits);
+    	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--) {
+		ellDoubleProj(pt1,cp);
+		hb = bitval(t0, b);
+		if(b < klen) kb = bitval(k, b); else kb = 0;
+		if((hb != 0) && (kb == 0))
+			ellAddProj(pt1, pt0, cp);
+		else if((hb == 0) && (kb !=0))
+			ellSubProj(pt1, pt0, cp);
+	}
+	if(ksign < 0) {
+		ellNegProj(pt1, cp);
+		k->sign = -k->sign;
+	}
+	returnGiant(t0);
+}
+
+void normalizeProj(pointProj pt, curveParams *cp)
+/* 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;
+	giant t1;
+
+	CKASSERT(cp->curveType == FCT_Weierstrass);
+	if(isZero(z)) {
+		int_to_giant(1,x); int_to_giant(1,y);
+		return;
+	}
+	t1 = borrowGiant(cp->maxDigits);
+	binvg_cp(cp, z);		// was binvaux(p, z);
+		gtog(z, t1);
+	gsquare(z); feemod(cp, z);
+	mulg(z, x); feemod(cp, x);
+	mulg(t1, z); mulg(z, y); feemod(cp, y);
+	int_to_giant(1, z);
+	returnGiant(t1);
+}
+
+static int
+jacobi_symbol(giant a, curveParams *cp)
+/* Standard Jacobi symbol (a/cp->basePrime).
+  basePrime must be odd, positive. */
+{
+	int t = 1, u;
+	giant t5 = borrowGiant(cp->maxDigits);
+	giant t6 = borrowGiant(cp->maxDigits);
+	giant t7 = borrowGiant(cp->maxDigits);
+	int rtn;
+
+	gtog(a, t5); feemod(cp, t5);
+	gtog(cp->basePrime, t6);
+	while(!isZero(t5)) {
+	    u = (t6->n[0]) & 7;
+		while((t5->n[0] & 1) == 0) {
+			gshiftright(1, t5);
+			if((u==3) || (u==5)) t = -t;
+		}
+		gtog(t5, t7); gtog(t6, t5); gtog(t7, t6);
+		u = (t6->n[0]) & 3;
+		if(((t5->n[0] & 3) == 3) && ((u & 3) == 3)) t = -t;
+		modg(t6, t5);
+	}
+	if(isone(t6)) {
+		rtn = t;
+	}
+	else {
+		rtn = 0;
+	}
+	returnGiant(t5);
+	returnGiant(t6);
+	returnGiant(t7);
+
+	return rtn;
+}
+
+static void
+powFp2(giant a, giant b, giant w2, giant n, curveParams *cp)
+/* Perform powering in the field F_p^2:
+   a + b w := (a + b w)^n (mod p), where parameter w2 is a quadratic
+   nonresidue (formally equal to w^2).
+ */
+{
+	int j;
+	giant t6;
+	giant t7;
+	giant t8;
+	giant t9;
+
+	if(isZero(n)) {
+		int_to_giant(1,a);
+		int_to_giant(0,b);
+		return;
+	}
+    	t6 = borrowGiant(cp->maxDigits);
+    	t7 = borrowGiant(cp->maxDigits);
+    	t8 = borrowGiant(cp->maxDigits);
+    	t9 = borrowGiant(cp->maxDigits);
+	gtog(a, t8); gtog(b, t9);
+	for(j = bitlen(n)-2; j >= 0; j--) {
+		gtog(b, t6);
+		mulg(a, b); addg(b,b); feemod(cp, b);  /* b := 2 a b. */
+		gsquare(t6); feemod(cp, t6);
+		mulg(w2, t6); feemod(cp, t6);
+		gsquare(a); addg(t6, a); feemod(cp, a);
+						/* a := a^2 + b^2 w2. */
+		if(bitval(n, j)) {
+			gtog(b, t6); mulg(t8, b); feemod(cp, b);
+			gtog(a, t7); mulg(t9, a); addg(a, b); feemod(cp, b);
+			mulg(t9, t6); feemod(cp, t6);
+			mulg(w2, t6); feemod(cp, t6);
+			mulg(t8, a); addg(t6, a); feemod(cp, a);
+		}
+	}
+	returnGiant(t6);
+	returnGiant(t7);
+	returnGiant(t8);
+	returnGiant(t9);
+	return;
+}
+
+static void
+powermodg(
+	giant		x,
+	giant		n,
+	curveParams	*cp
+)
+/* x becomes x^n (mod basePrime). */
+{
+	int 		len, pos;
+	giant		scratch2 = borrowGiant(cp->maxDigits);
+
+	gtog(x, scratch2);
+	int_to_giant(1, x);
+	len = bitlen(n);
+	pos = 0;
+	while (1)
+	{
+		if (bitval(n, pos++))
+		{
+			mulg(scratch2, x);
+			feemod(cp, x);
+		}
+		if (pos>=len)
+			break;
+		gsquare(scratch2);
+		feemod(cp, scratch2);
+	}
+	returnGiant(scratch2);
+}
+
+static int sqrtmod(giant x, curveParams *cp)
+/* If Sqrt[x] (mod p) exists, function returns 1, else 0.
+   In either case x is modified, but if 1 is returned,
+   x:= Sqrt[x] (mod p).
+ */
+{
+	int rtn;
+	giant t0 = borrowGiant(cp->maxDigits);
+	giant t1 = borrowGiant(cp->maxDigits);
+	giant t2 = borrowGiant(cp->maxDigits);
+	giant t3 = borrowGiant(cp->maxDigits);
+	giant t4 = borrowGiant(cp->maxDigits);
+
+	giant p = cp->basePrime;
+
+    	feemod(cp, x);			/* Justify the argument. */
+    	gtog(x, t0);  /* Store x for eventual validity check on square root. */
+    	if((p->n[0] & 3) == 3) {  /* The case p = 3 (mod 4). */
+		gtog(p, t1);
+		iaddg(1, t1); gshiftright(2, t1);
+		powermodg(x, t1, cp);
+		goto resolve;
+    	}
+	/* Next, handle case p = 5 (mod 8). */
+    	if((p->n[0] & 7) == 5) {
+		gtog(p, t1); int_to_giant(1, t2);
+		subg(t2, t1); gshiftright(2, t1);
+		gtog(x, t2);
+		powermodg(t2, t1, cp);  /* t2 := x^((p-1)/4) % p. */
+		iaddg(1, t1);
+		gshiftright(1, t1); /* t1 := (p+3)/8. */
+		if(isone(t2)) {
+			powermodg(x, t1, cp);  /* x^((p+3)/8) is root. */
+			goto resolve;
+		} else {
+			int_to_giant(1, t2); subg(t2, t1);
+				/* t1 := (p-5)/8. */
+			gshiftleft(2,x);
+			powermodg(x, t1, cp);
+			mulg(t0, x); addg(x, x); feemod(cp, x);
+				/* 2x (4x)^((p-5)/8. */
+			goto resolve;
+		}
+	}
+
+	/* Next, handle tougher case: p = 1 (mod 8). */
+	int_to_giant(2, t1);
+	while(1) {  /* Find appropriate nonresidue. */
+		gtog(t1, t2);
+		gsquare(t2); subg(x, t2); feemod(cp, t2);
+		if(jacobi_symbol(t2, cp) == -1) break;
+		iaddg(1, t1);
+	}  /* t2 is now w^2 in F_p^2. */
+   	int_to_giant(1, t3);
+   	gtog(p, t4); iaddg(1, t4); gshiftright(1, t4);
+	powFp2(t1, t3, t2, t4, cp);
+	gtog(t1, x);
+
+resolve:
+   	gtog(x,t1); gsquare(t1); feemod(cp, t1);
+    	if(gcompg(t0, t1) == 0) {
+		rtn = 1; 	/* Success. */
+	}
+	else {
+		rtn = 0;	/* no square root */
+	}
+	returnGiant(t0);
+	returnGiant(t1);
+	returnGiant(t2);
+	returnGiant(t3);
+	returnGiant(t4);
+	return rtn;
+}
+
+
+void findPointProj(pointProj pt, giant seed, curveParams *cp)
+/* Starting with seed, finds a random (projective) point {x,y,1} on curve.
+ */
+{
+	giant x = pt->x, y = pt->y, z = pt->z;
+
+	CKASSERT(cp->curveType == FCT_Weierstrass);
+	feemod(cp, seed);
+    	while(1) {
+		gtog(seed, x);
+		gsquare(x); feemod(cp, x);	// x := seed^2
+		addg(cp->a, x);			// x := seed^2 + a
+		mulg(seed,x); 			// x := seed^3 + a*seed
+		addg(cp->b, x);
+		feemod(cp, x);			// x := seed^3 + a seed + b.
+		/* test cubic form for having root. */
+		if(sqrtmod(x, cp)) break;
+		iaddg(1, seed);
+	}
+	gtog(x, y);
+    	gtog(seed,x);
+	int_to_giant(1, z);
+}
+
+#endif	/* CRYPTKIT_ELL_PROJ_ENABLE */