Security/libsecurity_cryptkit/lib/CurveParamDocs/fmodule.c

   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
  48 extern int pr[NUM_PRIMES];  /* Primes array from tools.c. */
  49
  50 unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES];
  51 int             modmode = GENERAL_MOD;
  52 int             curshorts = 0;
  53 static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL;
  54 static giant An = NULL, Ad = NULL;
  55 static point_mont pt1, pt2;
  56 point_mont pb[D+2];
  57 giant xzb[D+2];
  58
  59 static int verbosity = 0;
  60
  61 /**************************************************************
  62  *
  63  *      Functions
  64  *
  65  **************************************************************/
  66
  67 int
  68 init_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
  81 void
  82 verbose(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
  90 void
  91 dot(void)
  92 {
  93         printf(".");
  94         fflush(stdout);
  95 }
  96
  97 void
  98 set_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
 109 void
 110 special_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
 129 unsigned short *
 130 prime_list() {
 131         return(&factors[0]);
 132 }
 133
 134 unsigned short *
 135 exponent_list() {
 136         return(&exponents[0]);
 137 }
 138
 139 int
 140 sieve(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
 173 int
 174 pollard_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
 214 int
 215 pollard_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
 269 int
 270 ecm(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