1/************************************************************** 2 * 3 * fmodule.c 4 * 5 * Factoring utilities. 6 * 7 * Updates: 8 * 13 Apr 98 REC - creation 9 * 10 * c. 1998 Perfectly Scientific, Inc. 11 * All Rights Reserved. 12 * 13 * 14 *************************************************************/ 15 16/* include files */ 17 18#include <stdio.h> 19#include <math.h> 20#include <stdlib.h> 21#include <time.h> 22#ifdef _WIN32 23 24#include <process.h> 25 26#endif 27 28#include <string.h> 29#include "giants.h" 30#include "fmodule.h" 31#include "ellmont.h" 32 33#define NUM_PRIMES 6542 /* PrimePi[2^16]. */ 34#define GENERAL_MOD 0 35#define FERMAT_MOD 1 36#define MERSENNE_MOD (-1) 37#define D 100 /* A decimation parameter for stage-2 ECM factoring. */ 38 39/* compiler options */ 40 41#ifdef _WIN32 42#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */ 43#endif 44 45 46/* global variables */ 47 48extern int pr[NUM_PRIMES]; /* Primes array from tools.c. */ 49 50unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES]; 51int modmode = GENERAL_MOD; 52int curshorts = 0; 53static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL; 54static giant An = NULL, Ad = NULL; 55static point_mont pt1, pt2; 56point_mont pb[D+2]; 57giant xzb[D+2]; 58 59static int verbosity = 0; 60 61/************************************************************** 62 * 63 * Functions 64 * 65 **************************************************************/ 66 67int 68init_fmodule(int shorts) { 69 curshorts = shorts; 70 pb[0] = NULL; /* To guarantee proper ECM initialization. */ 71 t1 = newgiant(shorts); 72 t2 = newgiant(shorts); 73 t3 = newgiant(shorts); 74 t4 = newgiant(shorts); 75 An = newgiant(shorts); 76 Ad = newgiant(shorts); 77 pt1 = new_point_mont(shorts); 78 pt2 = new_point_mont(shorts); 79} 80 81void 82verbose(int state) 83/* Call verbose(1) for output during factoring processes, 84 call verbose(0) to silence all that. 85 */ 86{ 87 verbosity = state; 88} 89 90void 91dot(void) 92{ 93 printf("."); 94 fflush(stdout); 95} 96 97void 98set_mod_mode(int mode) 99/* Call this with mode := 1, 0, -1, for Fermat-mod, general mod, and Mersenne mod, 100 respectively; the point being that the special cases of 101 Fermat- and Mersenne-mod are much faster than 102 general mod. If all mods will be with respect to a number-to-be-factored, 103 of the form N = 2^m + 1, use Fermat mod; while if N = 2^m-1, use Mersenne mod. 104 */ 105{ 106 modmode = mode; 107} 108 109void 110special_modg( 111 giant N, 112 giant t 113) 114{ 115 switch (modmode) 116 { 117 case MERSENNE_MOD: 118 mersennemod(bitlen(N), t); 119 break; 120 case FERMAT_MOD: 121 fermatmod(bitlen(N)-1, t); 122 break; 123 default: 124 modg(N, t); 125 break; 126 } 127} 128 129unsigned short * 130prime_list() { 131 return(&factors[0]); 132} 133 134unsigned short * 135exponent_list() { 136 return(&exponents[0]); 137} 138 139int 140sieve(giant N, int sievelim) 141/* Returns number of N's prime factors < min(sievelim, 2^16), 142 with N reduced accordingly by said factors. 143 The n-th entry of factors[] becomes the n-th prime 144 factor of N, with corresponding exponent 145 becoming the n-th element of exponents[]. 146 */ 147{ int j, pcount, rem; 148 unsigned short pri; 149 150 pcount = 0; 151 exponents[0] = 0; 152 for (j=0; j < NUM_PRIMES; j++) 153 { 154 pri = pr[j]; 155 if(pri > sievelim) break; 156 do { 157 gtog(N, t1); 158 rem = idivg(pri, t1); 159 if(rem == 0) { 160 ++exponents[pcount]; 161 gtog(t1, N); 162 } 163 } while(rem == 0); 164 if(exponents[pcount] > 0) { 165 factors[pcount] = pr[j]; 166 ++pcount; 167 exponents[pcount] = 0; 168 } 169 } 170 return(pcount); 171} 172 173int 174pollard_rho(giant N, giant fact, int steps, int abort) 175/* Perform Pollard-rho x:= 3; loop(x:= x^2 + 2), a total of steps times. 176 Parameter fact will be a nontrivial factor found, in which case 177 N is also modified as: N := N/fact. 178 The function returns 0 if no nontrivial factor found, else returns 1. 179 The abort parameter, when set, causes the factorer to exit on the 180 first nontrivial factor found (the requisite GCD is checked 181 every 1000 steps). If abort := 0, the full number 182 of steps are always performed, then one solitary GCD is taken, 183 before exit. 184 */ 185{ 186 int j, found = 0; 187 188 itog(3, t1); 189 gtog(t1, t2); 190 itog(1, fact); 191 for(j=0; j < steps; j++) { 192 squareg(t1); iaddg(2, t1); special_modg(N, t1); 193 squareg(t2); iaddg(2, t2); special_modg(N, t2); 194 squareg(t2); iaddg(2, t2); special_modg(N, t2); 195 gtog(t2, t3); subg(t1,t3); mulg(t3, fact); special_modg(N, fact); 196 if(((j % 1000 == 999) && abort) || (j == steps-1)) { 197 if(verbosity) dot(); 198 gcdg(N, fact); 199 if(!isone(fact)) { 200 found = (gcompg(N, fact) == 0) ? 0 : 1; 201 break; 202 } 203 } 204 } 205 if(verbosity) { printf("\n"); fflush(stdout); } 206 if(found) { 207 divg(fact, N); 208 return(1); 209 } 210 itog(1, fact); 211 return(0); 212} 213 214int 215pollard_pminus1(giant N, giant fact, int lim, int abort) 216/* Perform Pollard-(p-1); where we test 217 218 GCD[N, 3^P - 1], 219 220 where P is an accumulation of primes <= min(lim, 2^16), 221 to appropriate powers. 222 Parameter fact will be a nontrivial factor found, in which case 223 N is also modified as: N := N/fact. 224 The function returns 0 if no nontrivial factor found, else returns 1. 225 The abort parameter, when set, causes the factorer to exit on the 226 first nontrivial factor found (the requisite GCD is checked 227 every 100 steps). If abort := 0, the full number 228 of steps are always performed, then one solitary GCD is taken, 229 before exit. 230 */ 231{ int cnt, j, k, pri, found = 0; 232 233 itog(3, fact); 234 for (j=0; j< NUM_PRIMES; j++) 235 { 236 pri = pr[j]; 237 if((pri > lim) || (j == NUM_PRIMES-1) || (abort && (j % 100 == 99))) { 238 if(verbosity) dot(); 239 gtog(fact, t1); 240 itog(1, t2); 241 subg(t2, t1); 242 special_modg(N, t1); 243 gcdg(N, t1); 244 if(!isone(t1)) { 245 found = (gcompg(N, t1) == 0) ? 0 : 1; 246 break; 247 } 248 if(pri > lim) break; 249 } 250 if(pri < 19) { cnt = 20-pri; /* Smaller primes get higher powers. */ 251 } else if(pri < 100) { 252 cnt = 2; 253 } else cnt = 1; 254 for (k=0; k< cnt; k++) 255 { 256 powermod(fact, pri, N); 257 } 258 } 259 if(verbosity) { printf("\n"); fflush(stdout); } 260 if(found) { 261 gtog(t1, fact); 262 divg(fact, N); 263 return(1); 264 } 265 itog(1, fact); 266 return(0); 267} 268 269int 270ecm(giant N, giant fact, int S, unsigned int B, unsigned int C) 271/* Perform elliptic curve method (ECM), with: 272 Brent seed parameter = S 273 Stage-one limit = B 274 Stage-two limit = C 275 This function: 276 returns 1 if nontrivial factor is found in stage 1 of ECM; 277 returns 2 if nontrivial factor is found in stage 2 of ECM; 278 returns 0 otherwise. 279 In the positive return, parameter fact is the factor and N := N/fact. 280 */ 281{ unsigned int pri, q; 282 int j, cnt, count, k; 283 284 if(verbosity) { 285 printf("Finding curve and point, B = %u, C = %u, seed = %d...", B, C, S); 286 fflush(stdout); 287 } 288 find_curve_point_brent(pt1, S, An, Ad, N); 289 if(verbosity) { 290 printf("\n"); 291 printf("Commencing stage 1 of ECM...\n"); 292 fflush(stdout); 293 } 294 295 q = pr[NUM_PRIMES-1]; 296 count = 0; 297 for(j=0; ; j++) { 298 if(j < NUM_PRIMES) { 299 pri = pr[j]; 300 } else { 301 q += 2; 302 if(primeq(q)) pri = q; 303 else continue; 304 } 305 if(verbosity) if((++count) % 100 == 0) dot(); 306 if(pri > B) break; 307 if(pri < 19) { cnt = 20-pri; 308 } else if(pri < 100) { 309 cnt = 2; 310 } else cnt = 1; 311 for(k = 0; k < cnt; k++) 312 ell_mul_int_brent(pt1, pri, An, Ad, N); 313 } 314 k = 19; 315 while (k<B) 316 { 317 if (primeq(k)) 318 { 319 ell_mul_int_brent(pt1, k, An, Ad, N); 320 if (k<100) 321 ell_mul_int_brent(pt1, k, An, Ad, N); 322 if (cnt++ %100==0) 323 dot(); 324 } 325 k += 2; 326 } 327 if(verbosity) { printf("\n"); fflush(stdout); } 328 329/* Next, test stage-1 attempt. */ 330 gtog(pt1->z, fact); 331 gcdg(N, fact); 332 if((!isone(fact)) && (gcompg(N, fact) != 0)) { 333 divg(fact, N); 334 return(1); 335 } 336 if(B >= C) { /* No stage 2 planned. */ 337 itog(1, fact); 338 return(0); 339 } 340 341/* Commence second stage of ECM. */ 342 if(verbosity) { 343 printf("\n"); 344 printf("Commencing stage 2 of ECM...\n"); 345 fflush(stdout); 346 } 347 if(pb[0] == NULL) { 348 for(k=0; k < D+2; k++) { 349 pb[k] = new_point_mont(curshorts); 350 xzb[k] = newgiant(curshorts); 351 352 } 353 } 354 k = ((int)B)/D; 355 ptop_mont(pt1, pb[0]); 356 ell_mul_int_brent(pb[0], k*D+1 , An, Ad, N); 357 ptop_mont(pt1, pb[D+1]); 358 ell_mul_int_brent(pb[D+1], (k+2)*D+1 , An, Ad, N); 359 360 for (j=1; j <= D; j++) 361 { 362 ptop_mont(pt1, pb[j]); 363 ell_mul_int_brent(pb[j], 2*j , An, Ad, N); 364 gtog(pb[j]->z, xzb[j]); 365 mulg(pb[j]->x, xzb[j]); 366 special_modg(N, xzb[j]); 367 } 368 itog(1, fact); 369 count = 0; 370 while (1) { 371 if(verbosity) if((++count) % 10 == 0) dot(); 372 gtog(pb[0]->z, xzb[0]); 373 mulg(pb[0]->x, xzb[0]); 374 special_modg(N, xzb[0]); 375 mulg(pb[0]->z, fact); 376 special_modg(N, fact); /* Accumulate. */ 377 for (j = 1; j < D; j++) { 378 if (!primeq(k*D+1+2*j)) continue; 379/* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */ 380 gtog(pb[0]->x, t1); 381 subg(pb[j]->x, t1); 382 gtog(pb[0]->z, t2); 383 addg(pb[j]->z, t2); 384 mulg(t1, t2); 385 special_modg(N, t1); 386 subg(xzb[0], t2); 387 addg(xzb[j], t2); 388 special_modg(N, t2); 389 mulg(t2, fact); 390 special_modg(N, fact); 391 } 392 k += 2; 393 if(k*D > C) 394 break; 395 ptop_mont(pb[D+1], pt2); 396 ell_odd_brent(pb[D], pb[D+1], pb[0], N); 397 ptop_mont(pt2, pb[0]); 398 } 399 if(verbosity) { printf("\n"); fflush(stdout); } 400 401 gcdg(N, fact); 402 if((!isone(fact)) && (gcompg(N, fact) != 0)) { 403 divg(fact, N); 404 return(2); /* Return value of 2 for stage-2 success! */ 405 } 406 itog(1, fact); 407 return(0); 408} 409 410 411