]> git.saurik.com Git - apple/security.git/blobdiff - Security/libsecurity_cryptkit/lib/CurveParamDocs/ellproj.c
Security-57031.1.35.tar.gz
[apple/security.git] / Security / libsecurity_cryptkit / lib / CurveParamDocs / ellproj.c
diff --git a/Security/libsecurity_cryptkit/lib/CurveParamDocs/ellproj.c b/Security/libsecurity_cryptkit/lib/CurveParamDocs/ellproj.c
new file mode 100644 (file)
index 0000000..b98b257
--- /dev/null
@@ -0,0 +1,448 @@
+/**************************************************************
+ *
+ * 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);
+}