+/**************************************************************
+ *
+ * fmodule.c
+ *
+ * Factoring utilities.
+ *
+ * Updates:
+ * 13 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 "fmodule.h"
+#include "ellmont.h"
+
+#define NUM_PRIMES 6542 /* PrimePi[2^16]. */
+#define GENERAL_MOD 0
+#define FERMAT_MOD 1
+#define MERSENNE_MOD (-1)
+#define D 100 /* A decimation parameter for stage-2 ECM factoring. */
+
+/* compiler options */
+
+#ifdef _WIN32
+#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
+#endif
+
+
+/* global variables */
+
+extern int pr[NUM_PRIMES]; /* Primes array from tools.c. */
+
+unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES];
+int modmode = GENERAL_MOD;
+int curshorts = 0;
+static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL;
+static giant An = NULL, Ad = NULL;
+static point_mont pt1, pt2;
+point_mont pb[D+2];
+giant xzb[D+2];
+
+static int verbosity = 0;
+
+/**************************************************************
+ *
+ * Functions
+ *
+ **************************************************************/
+
+int
+init_fmodule(int shorts) {
+ curshorts = shorts;
+ pb[0] = NULL; /* To guarantee proper ECM initialization. */
+ t1 = newgiant(shorts);
+ t2 = newgiant(shorts);
+ t3 = newgiant(shorts);
+ t4 = newgiant(shorts);
+ An = newgiant(shorts);
+ Ad = newgiant(shorts);
+ pt1 = new_point_mont(shorts);
+ pt2 = new_point_mont(shorts);
+}
+
+void
+verbose(int state)
+/* Call verbose(1) for output during factoring processes,
+ call verbose(0) to silence all that.
+ */
+{
+ verbosity = state;
+}
+
+void
+dot(void)
+{
+ printf(".");
+ fflush(stdout);
+}
+
+void
+set_mod_mode(int mode)
+/* Call this with mode := 1, 0, -1, for Fermat-mod, general mod, and Mersenne mod,
+ respectively; the point being that the special cases of
+ Fermat- and Mersenne-mod are much faster than
+ general mod. If all mods will be with respect to a number-to-be-factored,
+ of the form N = 2^m + 1, use Fermat mod; while if N = 2^m-1, use Mersenne mod.
+ */
+{
+ modmode = mode;
+}
+
+void
+special_modg(
+ giant N,
+ giant t
+)
+{
+ switch (modmode)
+ {
+ case MERSENNE_MOD:
+ mersennemod(bitlen(N), t);
+ break;
+ case FERMAT_MOD:
+ fermatmod(bitlen(N)-1, t);
+ break;
+ default:
+ modg(N, t);
+ break;
+ }
+}
+
+unsigned short *
+prime_list() {
+ return(&factors[0]);
+}
+
+unsigned short *
+exponent_list() {
+ return(&exponents[0]);
+}
+
+int
+sieve(giant N, int sievelim)
+/* Returns number of N's prime factors < min(sievelim, 2^16),
+ with N reduced accordingly by said factors.
+ The n-th entry of factors[] becomes the n-th prime
+ factor of N, with corresponding exponent
+ becoming the n-th element of exponents[].
+ */
+{ int j, pcount, rem;
+ unsigned short pri;
+
+ pcount = 0;
+ exponents[0] = 0;
+ for (j=0; j < NUM_PRIMES; j++)
+ {
+ pri = pr[j];
+ if(pri > sievelim) break;
+ do {
+ gtog(N, t1);
+ rem = idivg(pri, t1);
+ if(rem == 0) {
+ ++exponents[pcount];
+ gtog(t1, N);
+ }
+ } while(rem == 0);
+ if(exponents[pcount] > 0) {
+ factors[pcount] = pr[j];
+ ++pcount;
+ exponents[pcount] = 0;
+ }
+ }
+ return(pcount);
+}
+
+int
+pollard_rho(giant N, giant fact, int steps, int abort)
+/* Perform Pollard-rho x:= 3; loop(x:= x^2 + 2), a total of steps times.
+ Parameter fact will be a nontrivial factor found, in which case
+ N is also modified as: N := N/fact.
+ The function returns 0 if no nontrivial factor found, else returns 1.
+ The abort parameter, when set, causes the factorer to exit on the
+ first nontrivial factor found (the requisite GCD is checked
+ every 1000 steps). If abort := 0, the full number
+ of steps are always performed, then one solitary GCD is taken,
+ before exit.
+ */
+{
+ int j, found = 0;
+
+ itog(3, t1);
+ gtog(t1, t2);
+ itog(1, fact);
+ for(j=0; j < steps; j++) {
+ squareg(t1); iaddg(2, t1); special_modg(N, t1);
+ squareg(t2); iaddg(2, t2); special_modg(N, t2);
+ squareg(t2); iaddg(2, t2); special_modg(N, t2);
+ gtog(t2, t3); subg(t1,t3); mulg(t3, fact); special_modg(N, fact);
+ if(((j % 1000 == 999) && abort) || (j == steps-1)) {
+ if(verbosity) dot();
+ gcdg(N, fact);
+ if(!isone(fact)) {
+ found = (gcompg(N, fact) == 0) ? 0 : 1;
+ break;
+ }
+ }
+ }
+ if(verbosity) { printf("\n"); fflush(stdout); }
+ if(found) {
+ divg(fact, N);
+ return(1);
+ }
+ itog(1, fact);
+ return(0);
+}
+
+int
+pollard_pminus1(giant N, giant fact, int lim, int abort)
+/* Perform Pollard-(p-1); where we test
+
+ GCD[N, 3^P - 1],
+
+ where P is an accumulation of primes <= min(lim, 2^16),
+ to appropriate powers.
+ Parameter fact will be a nontrivial factor found, in which case
+ N is also modified as: N := N/fact.
+ The function returns 0 if no nontrivial factor found, else returns 1.
+ The abort parameter, when set, causes the factorer to exit on the
+ first nontrivial factor found (the requisite GCD is checked
+ every 100 steps). If abort := 0, the full number
+ of steps are always performed, then one solitary GCD is taken,
+ before exit.
+ */
+{ int cnt, j, k, pri, found = 0;
+
+ itog(3, fact);
+ for (j=0; j< NUM_PRIMES; j++)
+ {
+ pri = pr[j];
+ if((pri > lim) || (j == NUM_PRIMES-1) || (abort && (j % 100 == 99))) {
+ if(verbosity) dot();
+ gtog(fact, t1);
+ itog(1, t2);
+ subg(t2, t1);
+ special_modg(N, t1);
+ gcdg(N, t1);
+ if(!isone(t1)) {
+ found = (gcompg(N, t1) == 0) ? 0 : 1;
+ break;
+ }
+ if(pri > lim) break;
+ }
+ if(pri < 19) { cnt = 20-pri; /* Smaller primes get higher powers. */
+ } else if(pri < 100) {
+ cnt = 2;
+ } else cnt = 1;
+ for (k=0; k< cnt; k++)
+ {
+ powermod(fact, pri, N);
+ }
+ }
+ if(verbosity) { printf("\n"); fflush(stdout); }
+ if(found) {
+ gtog(t1, fact);
+ divg(fact, N);
+ return(1);
+ }
+ itog(1, fact);
+ return(0);
+}
+
+int
+ecm(giant N, giant fact, int S, unsigned int B, unsigned int C)
+/* Perform elliptic curve method (ECM), with:
+ Brent seed parameter = S
+ Stage-one limit = B
+ Stage-two limit = C
+ This function:
+ returns 1 if nontrivial factor is found in stage 1 of ECM;
+ returns 2 if nontrivial factor is found in stage 2 of ECM;
+ returns 0 otherwise.
+ In the positive return, parameter fact is the factor and N := N/fact.
+ */
+{ unsigned int pri, q;
+ int j, cnt, count, k;
+
+ if(verbosity) {
+ printf("Finding curve and point, B = %u, C = %u, seed = %d...", B, C, S);
+ fflush(stdout);
+ }
+ find_curve_point_brent(pt1, S, An, Ad, N);
+ if(verbosity) {
+ printf("\n");
+ printf("Commencing stage 1 of ECM...\n");
+ fflush(stdout);
+ }
+
+ q = pr[NUM_PRIMES-1];
+ count = 0;
+ for(j=0; ; j++) {
+ if(j < NUM_PRIMES) {
+ pri = pr[j];
+ } else {
+ q += 2;
+ if(primeq(q)) pri = q;
+ else continue;
+ }
+ if(verbosity) if((++count) % 100 == 0) dot();
+ if(pri > B) break;
+ if(pri < 19) { cnt = 20-pri;
+ } else if(pri < 100) {
+ cnt = 2;
+ } else cnt = 1;
+ for(k = 0; k < cnt; k++)
+ ell_mul_int_brent(pt1, pri, An, Ad, N);
+ }
+ k = 19;
+ while (k<B)
+ {
+ if (primeq(k))
+ {
+ ell_mul_int_brent(pt1, k, An, Ad, N);
+ if (k<100)
+ ell_mul_int_brent(pt1, k, An, Ad, N);
+ if (cnt++ %100==0)
+ dot();
+ }
+ k += 2;
+ }
+ if(verbosity) { printf("\n"); fflush(stdout); }
+
+/* Next, test stage-1 attempt. */
+ gtog(pt1->z, fact);
+ gcdg(N, fact);
+ if((!isone(fact)) && (gcompg(N, fact) != 0)) {
+ divg(fact, N);
+ return(1);
+ }
+ if(B >= C) { /* No stage 2 planned. */
+ itog(1, fact);
+ return(0);
+ }
+
+/* Commence second stage of ECM. */
+ if(verbosity) {
+ printf("\n");
+ printf("Commencing stage 2 of ECM...\n");
+ fflush(stdout);
+ }
+ if(pb[0] == NULL) {
+ for(k=0; k < D+2; k++) {
+ pb[k] = new_point_mont(curshorts);
+ xzb[k] = newgiant(curshorts);
+
+ }
+ }
+ k = ((int)B)/D;
+ ptop_mont(pt1, pb[0]);
+ ell_mul_int_brent(pb[0], k*D+1 , An, Ad, N);
+ ptop_mont(pt1, pb[D+1]);
+ ell_mul_int_brent(pb[D+1], (k+2)*D+1 , An, Ad, N);
+
+ for (j=1; j <= D; j++)
+ {
+ ptop_mont(pt1, pb[j]);
+ ell_mul_int_brent(pb[j], 2*j , An, Ad, N);
+ gtog(pb[j]->z, xzb[j]);
+ mulg(pb[j]->x, xzb[j]);
+ special_modg(N, xzb[j]);
+ }
+ itog(1, fact);
+ count = 0;
+ while (1) {
+ if(verbosity) if((++count) % 10 == 0) dot();
+ gtog(pb[0]->z, xzb[0]);
+ mulg(pb[0]->x, xzb[0]);
+ special_modg(N, xzb[0]);
+ mulg(pb[0]->z, fact);
+ special_modg(N, fact); /* Accumulate. */
+ for (j = 1; j < D; j++) {
+ if (!primeq(k*D+1+2*j)) continue;
+/* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */
+ gtog(pb[0]->x, t1);
+ subg(pb[j]->x, t1);
+ gtog(pb[0]->z, t2);
+ addg(pb[j]->z, t2);
+ mulg(t1, t2);
+ special_modg(N, t1);
+ subg(xzb[0], t2);
+ addg(xzb[j], t2);
+ special_modg(N, t2);
+ mulg(t2, fact);
+ special_modg(N, fact);
+ }
+ k += 2;
+ if(k*D > C)
+ break;
+ ptop_mont(pb[D+1], pt2);
+ ell_odd_brent(pb[D], pb[D+1], pb[0], N);
+ ptop_mont(pt2, pb[0]);
+ }
+ if(verbosity) { printf("\n"); fflush(stdout); }
+
+ gcdg(N, fact);
+ if((!isone(fact)) && (gcompg(N, fact) != 0)) {
+ divg(fact, N);
+ return(2); /* Return value of 2 for stage-2 success! */
+ }
+ itog(1, fact);
+ return(0);
+}
+
+
projects / apple / security.git / blobdiff