+++ /dev/null
-/**************************************************************
- *
- * 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);
-}
-
-