--- /dev/null
+/**************************************************************
+ *
+ * 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 <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#ifdef _WIN32
+
+#include <process.h>
+
+#endif
+
+#include <string.h>
+#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<D+2;j++)
+ {
+ xb[j] = newgiant(nsh);
+ zb[j] = newgiant(nsh);
+ xzb[j] = newgiant(nsh);
+ }
+}
+
+int
+bigprimeq(
+ giant x
+)
+{
+ itog(1, t1);
+ gtog(x, t2);
+ subg(t1, t2);
+ itog(5, t1);
+ powermodg(t1, t2, x);
+ if (isone(t1))
+ return(1);
+ return(0);
+}
+
+
+void
+dot(
+ void
+)
+{
+ printf(".");
+ fflush(stdout);
+}
+
+/**************************************************************
+ *
+ * Main Function
+ *
+ **************************************************************/
+
+main(
+ int argc,
+ char *argv[]
+)
+{
+ int j, k, C, nshorts, cnt, count,
+ limitbits = 0, pass, npr, rem;
+ long B;
+ int randmode = 0;
+
+ if (!strcmp(argv[argc-1], "-r"))
+ {
+ randmode = 1;
+ if (argc > 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;j<nshorts;j++)
+ {
+ ell_mul(xr, zr, 1<<16, An, Ad, N);
+ ell_mul(xr, zr, 3*3*3*3*3*3*3*3*3*3*3, An, Ad, N);
+ ell_mul(xr, zr, 5*5*5*5*5*5*5, An, Ad, N);
+ ell_mul(xr, zr, 7*7*7*7*7*7, An, Ad, N);
+ ell_mul(xr, zr, 11*11*11*11, An, Ad, N);
+ ell_mul(xr, zr, 13*13*13*13, An, Ad, N);
+ ell_mul(xr, zr, 17*17*17, An, Ad, N);
+ }
+ k = 19;
+ while (k<B)
+ {
+ if (isprime(k))
+ {
+ ell_mul(xr, zr, k, An, Ad, N);
+ if (k<100)
+ ell_mul(xr, zr, k, An, Ad, N);
+ if (cnt++ %100==0)
+ dot();
+ }
+ k += 2;
+ }
+ count = 0;
+
+ gtog(zr, gg);
+ 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;
+ }
+ 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;
+}
+
projects / apple / security.git / blobdiff