X-Git-Url: https://git.saurik.com/apple/security.git/blobdiff_plain/5dd5f9ec28f304ca377c42fd7f711d6cf12b90e1..5c19dc3ae3bd8e40a9c028b0deddd50ff337692c:/OSX/libsecurity_cryptkit/lib/CurveParamDocs/factor.c?ds=inline diff --git a/OSX/libsecurity_cryptkit/lib/CurveParamDocs/factor.c b/OSX/libsecurity_cryptkit/lib/CurveParamDocs/factor.c new file mode 100644 index 00000000..a0d0186f --- /dev/null +++ b/OSX/libsecurity_cryptkit/lib/CurveParamDocs/factor.c @@ -0,0 +1,844 @@ +/************************************************************** + * + * factor.c + * + * General purpose factoring program + * + * Updates: + * 18 May 97 REC - invoked new, fast divide functions + * 26 Apr 97 RDW - fixed tabs and unix EOL + * 20 Apr 97 RDW - conversion to TC4.5 + * + * c. 1997 Perfectly Scientific, Inc. + * All Rights Reserved. + * + * + *************************************************************/ + +/* include files */ + +#include +#include +#include +#include +#ifdef _WIN32 + +#include + +#endif + +#include +#include "giants.h" + + +/* definitions */ + +#define D 100 +#define NUM_PRIMES 6542 /* PrimePi[2^16]. */ + + +/* compiler options */ + +#ifdef _WIN32 +#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */ +#endif + + +/* global variables */ + +extern giant scratch2; +int pr[NUM_PRIMES]; +giant xr = NULL, xs = NULL, zs = NULL, zr = NULL, xorg = NULL, + zorg = NULL, t1 = NULL, t2 = NULL, t3 = NULL, N = NULL, + gg = NULL, An = NULL, Ad = NULL; +giant xb[D+2], zb[D+2], xzb[D+2]; +int modmode = 0, Q, modcount = 0; + + +/* function prototypes */ + +void ell_even(giant x1, giant z1, giant x2, giant z2, giant An, + giant Ad, giant N); +void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xor, + giant zor, giant N); +void ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N); +int least_2(int n); +void dot(void); +int psi_rand(); + + +/************************************************************** + * + * Functions + * + **************************************************************/ + + +int +psi_rand( + void +) +{ + unsigned short hi; + unsigned short low; + time_t tp; + int result; + + time(&tp); + low = (unsigned short)rand(); + hi = (unsigned short)rand(); + result = ((hi << 16) | low) ^ ((int)tp); + + return (result & 0x7fffffff); +} + + +void +set_random_seed( + void +) +{ + /* Start the random number generator at a new position. */ + time_t tp; + + time(&tp); + srand((int)tp + (int)getpid()); +} + + +int +isprime( + int odd +) +{ + int j; + int p; + + for (j=1; ; j++) + { + p = pr[j]; + if (p*p > odd) + return(1); + if (odd % p == 0) + return(0); + } +} + + +int +primeq( + int p +) +{ + register int j=3; + + if ((p&1)==0) + return ((p==2)?1:0); + if (j>=p) + return (1); + while ((p%j)!=0) + { + j+=2; + if (j*j>p) + return(1); + } + return(0); +} + + +void +s_modg( + giant N, + giant t +) +{ + ++modcount; + switch (modmode) + { + case 0: + modg(N, t); + break; + case -1: + mersennemod(Q, t); + break; + case 1: + fermatmod(Q, t); + break; + } +} + + +void +reset_mod( + giant x, + giant N +) +/* Perform a divide (by the discovered factor) and switch back + to non-Fermat-non-Mersenne (i.e. normal) mod mode. */ +{ + divg(x, N); + modmode = 0; +} + +void +ell_even( + giant x1, + giant z1, + giant x2, + giant z2, + giant An, + giant Ad, + giant N +) +{ + gtog(x1, t1); + addg(z1, t1); + squareg(t1); + s_modg(N, t1); + gtog(x1, t2); + subg(z1, t2); + squareg(t2); + s_modg(N, t2); + gtog(t1, t3); + subg(t2, t3); + gtog(t2, x2); + mulg(t1, x2); + gshiftleft(2, x2); + s_modg(N, x2); + mulg(Ad, x2); + s_modg(N, x2); + mulg(Ad, t2); + gshiftleft(2, t2); + s_modg(N, t2); + gtog(t3, t1); + mulg(An, t1); + s_modg(N, t1); + addg(t1, t2); + mulg(t3, t2); + s_modg(N, t2); + gtog(t2,z2); +} + + +void +ell_odd( + giant x1, + giant z1, + giant x2, + giant z2, + giant xor, + giant zor, + giant N +) +{ + gtog(x1, t1); + subg(z1, t1); + gtog(x2, t2); + addg(z2, t2); + mulg(t1, t2); + s_modg(N, t2); + gtog(x1, t1); + addg(z1, t1); + gtog(x2, t3); + subg(z2, t3); + mulg(t3, t1); + s_modg(N, t1); + gtog(t2, x2); + addg(t1, x2); + squareg(x2); + s_modg(N, x2); + gtog(t2, z2); + subg(t1, z2); + squareg(z2); + s_modg(N, z2); + mulg(zor, x2); + s_modg(N, x2); + mulg(xor, z2); + s_modg(N, z2); +} + + +void +ell_mul( + giant xx, + giant zz, + int n, + giant An, + giant Ad, + giant N +) +{ + unsigned int c = (unsigned int)0x80000000; + + if (n==1) + return; + if (n==2) + { + ell_even(xx, zz, xx, zz, An, Ad, N); + return; + } + gtog(xx, xorg); + gtog(zz, zorg); + ell_even(xx, zz, xs, zs, An, Ad, N); + + while((c&n) == 0) + { + c >>= 1; + } + + c>>=1; + do + { + if (c&n) + { + ell_odd(xs, zs, xx, zz, xorg, zorg, N); + ell_even(xs, zs, xs, zs, An, Ad, N); + } + else + { + ell_odd(xx, zz, xs, zs, xorg, zorg, N); + ell_even(xx, zz, xx, zz, An, Ad, N); + } + c >>= 1; + } while(c); +} + + + +/* From R. P. Brent, priv. comm. 1996: +Let s > 5 be a pseudo-random seed (called $\sigma$ in the Tech. Report), + + u/v = (s^2 - 5)/(4s) + +Then starting point is (x_1, y_1) where + + x_1 = (u/v)^3 +and + a = (v-u)^3(3u+v)/(4u^3 v) - 2 +*/ + +void +choose12( + giant x, + giant z, + int k, + giant An, + giant Ad, + giant N +) +{ + itog(k, zs); + gtog(zs, xs); + squareg(xs); + itog(5, t2); + subg(t2, xs); + s_modg(N, xs); + addg(zs, zs); + addg(zs, zs); + s_modg(N, zs); + gtog(xs, x); + squareg(x); + s_modg(N, x); + mulg(xs, x); + s_modg(N, x); + gtog(zs, z); + squareg(z); + s_modg(N, z); + mulg(zs, z); + s_modg(N, z); + + /* Now for A. */ + gtog(zs, t2); + subg(xs, t2); + gtog(t2, t3); + squareg(t2); + s_modg(N, t2); + mulg(t3, t2); + s_modg(N, t2); /* (v-u)^3. */ + gtog(xs, t3); + addg(t3, t3); + addg(xs, t3); + addg(zs, t3); + s_modg(N, t3); + mulg(t3, t2); + s_modg(N, t2); /* (v-u)^3 (3u+v). */ + gtog(zs, t3); + mulg(xs, t3); + s_modg(N, t3); + squareg(xs); + s_modg(N, xs); + mulg(xs, t3); + s_modg(N, t3); + addg(t3, t3); + addg(t3, t3); + s_modg(N, t3); + gtog(t3, Ad); + gtog(t2, An); /* An/Ad is now A + 2. */ +} + + +void +ensure( + int q +) +{ + int nsh, j; + + N = newgiant(INFINITY); + if(!q) + { + gread(N,stdin); + q = bitlen(N) + 1; + } + nsh = q/4; /* Allowing (easily) enough space per giant, + since N is about 2^q, which is q bits, or + q/16 shorts. But squaring, etc. is allowed, + so we need at least q/8, and we choose q/4 + to be conservative. */ + if (!xr) + xr = newgiant(nsh); + if (!zr) + zr = newgiant(nsh); + if (!xs) + xs = newgiant(nsh); + if (!zs) + zs = newgiant(nsh); + if (!xorg) + xorg = newgiant(nsh); + if (!zorg) + zorg = newgiant(nsh); + if (!t1) + t1 = newgiant(nsh); + if (!t2) + t2 = newgiant(nsh); + if (!t3) + t3 = newgiant(nsh); + if (!gg) + gg = newgiant(nsh); + if (!An) + An = newgiant(nsh); + if (!Ad) + Ad = newgiant(nsh); + for (j=0;j 4) + /* This segment only takes effect in random mode. */ + limitbits = atoi(argv[argc-2]); + } + else + { + randmode = 0; + } + + modmode = 0; + if (argc > 2) + { + modmode = atoi(argv[1]); + Q = atoi(argv[2]); + } + if (modmode==0) + Q = 0; + ensure(Q); + if (modmode) + { + itog(1, N); + gshiftleft(Q, N); + itog(modmode, t1); + addg(t1, N); + } + pr[0] = 2; + for (k=0, npr=1;; k++) + { + if (primeq(3+2*k)) + { + pr[npr++] = 3+2*k; + if (npr >= NUM_PRIMES) + break; + } + } + + if (randmode == 0) + { + printf("Sieving...\n"); + fflush(stdout); + for (j=0; j < NUM_PRIMES; j++) + { + gtog(N, t1); + rem = idivg(pr[j], t1); + if (rem == 0) + { + printf("%d ", pr[j]); + gtog(t1, N); + if (isone(N)) + { + printf("\n"); + exit(0); + } + else + { + printf("* "); + fflush(stdout); + } + --j; + } + } + + if (bigprimeq(N)) + { + gout(N); + exit(0); + } + + printf("\n"); + printf("Commencing Pollard rho...\n"); + fflush(stdout); + itog(1, gg); + itog(3, t1); itog(3, t2); + + for (j=0; j < 15000; j++) + { + if((j%100) == 0) + { + dot(); + gcdg(N, gg); + if (!isone(gg)) + break; + } + squareg(t1); + iaddg(2, t1); + s_modg(N, t1); + squareg(t2); + iaddg(2, t2); + s_modg(N, t2); + squareg(t2); + iaddg(2, t2); + s_modg(N, t2); + gtog(t2, t3); + subg(t1, t3); + t3->sign = abs(t3->sign); + mulg(t3, gg); + s_modg(N, gg); + } + gcdg(N, gg); + + if ((gcompg(N,gg) != 0) && (!isone(gg))) + { + fprintf(stdout,"\n"); + gout(gg); + reset_mod(gg, N); + if (isone(N)) + { + printf("\n"); + exit(0); + } + else + { + printf("* "); + fflush(stdout); + } + if (bigprimeq(N)) + { + gout(N); + exit(0); + } + } + + printf("\n"); + printf("Commencing Pollard (p-1)...\n"); + fflush(stdout); + itog(1, gg); + itog(3, t1); + for (j=0; j< NUM_PRIMES; j++) + { + cnt = (int)(8*log(2.0)/log(1.0*pr[j])); + if (cnt < 2) + cnt = 1; + for (k=0; k< cnt; k++) + { + powermod(t1, pr[j], N); + } + itog(1, t2); + subg(t1, t2); + mulg(t2, gg); + s_modg(N, gg); + + if (j % 100 == 0) + { + dot(); + gcdg(N, gg); + if (!isone(gg)) + break; + } + } + gcdg(N, gg); + if ((gcompg(N,gg) != 0) && (!isone(gg))) + { + fprintf(stdout,"\n"); + gout(gg); + reset_mod(gg, N); + if (isone(N)) + { + printf("\n"); + exit(0); + } + else + { + printf("* "); + fflush(stdout); + } + if (bigprimeq(N)) + { + gout(N); + exit(0); + } + } + } /* This is the end of (randmode == 0) */ + + printf("\n"); + printf("Commencing ECM...\n"); + fflush(stdout); + + if (randmode) + set_random_seed(); + pass = 0; + while (++pass) + { + if (randmode == 0) + { + if (pass <= 3) + { + B = 1000; + } + else if (pass <= 10) + { + B = 10000; + } + else if (pass <= 100) + { + B = 100000L; + } else + { + B = 1000000L; + } + } + else + { + B = 2000000L; + } + C = 50*((int)B); + + /* Next, choose curve with order divisible by 16 and choose + * a point (xr/zr) on said curve. + */ + + /* Order-div-12 case. + * cnt = 8020345; Brent's parameter for stage one discovery + * of 27-digit factor of F_13. + */ + + cnt = psi_rand(); /* cnt = 8020345; */ + choose12(xr, zr, cnt, An, Ad, N); + printf("Choosing curve %d, with s = %d, B = %d, C = %d:\n", pass,cnt, B, C); fflush(stdout); + cnt = 0; + nshorts = 1; + count = 0; + for (j=0;jlimitbits)) + { + fprintf(stdout,"\n"); + gwriteln(gg, stdout); + fflush(stdout); + reset_mod(gg, N); + if (isone(N)) + { + printf("\n"); + exit(0); + } + else + { + printf("* "); + fflush(stdout); + } + if (bigprimeq(N)) + { + gout(N); + exit(0); + } + continue; + } + else + { + printf("\n"); + fflush(stdout); + } + + /* Continue; Invoke, to test Stage 1 only. */ + k = ((int)B)/D; + gtog(xr, xb[0]); + gtog(zr, zb[0]); + ell_mul(xb[0], zb[0], k*D+1 , An, Ad, N); + gtog(xr, xb[D+1]); + gtog(zr, zb[D+1]); + ell_mul(xb[D+1], zb[D+1], (k+2)*D+1 , An, Ad, N); + + for (j=1; j <= D; j++) + { + gtog(xr, xb[j]); + gtog(zr, zb[j]); + ell_mul(xb[j], zb[j], 2*j , An, Ad, N); + gtog(zb[j], xzb[j]); + mulg(xb[j], xzb[j]); + s_modg(N, xzb[j]); + } + modcount = 0; + printf("\nCommencing second stage, curve %d...\n",pass); fflush(stdout); + count = 0; + itog(1, gg); + + while (1) + { + gtog(zb[0], xzb[0]); + mulg(xb[0], xzb[0]); + s_modg(N, xzb[0]); + mulg(zb[0], gg); + s_modg(N,gg); /* Accumulate. */ + for (j = 1; j < D; j++) + { + if (!isprime(k*D+1+ 2*j)) + continue; + + /* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */ + gtog(xb[0], t1); + subg(xb[j], t1); + gtog(zb[0], t2); + addg(zb[j], t2); + mulg(t1, t2); + s_modg(N, t1); + subg(xzb[0], t2); + addg(xzb[j], t2); + s_modg(N, t2); + --modcount; + mulg(t2, gg); + s_modg(N, gg); + if((++count)%1000==0) + dot(); + } + + k += 2; + if(k*D > C) + break; + gtog(xb[D+1], xs); + gtog(zb[D+1], zs); + ell_odd(xb[D], zb[D], xb[D+1], zb[D+1], xb[0], zb[0], N); + gtog(xs, xb[0]); + gtog(zs, zb[0]); + } + + gcdg(N, gg); + if((!isone(gg))&&(bitlen(gg)>limitbits)) + { + fprintf(stdout,"\n"); + gwriteln(gg, stdout); + fflush(stdout); + reset_mod(gg, N); + if (isone(N)) + { + printf("\n"); + exit(0); + } + else + { + printf("* "); + fflush(stdout); + } + if (bigprimeq(N)) + { + gout(N); + exit(0); + } + continue; + } + + printf("\n"); + fflush(stdout); + } + + return 0; +} +