home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / lenstra.c < prev    next >
Text File  |  1989-04-19  |  9KB  |  253 lines
  1. /*
  2.  *   Program to factor big numbers using Lenstras elliptic curve method.
  3.  *   Works when for some prime divisor p of n, p+1+d has only
  4.  *   small factors, where d depends on the particular curve used.
  5.  *   See "Speeding the Pollard and Elliptic Curve Methods"
  6.  *   by Peter Montgomery, Math. Comp. Vol. 48 Jan. 1987 pp243-264
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include "miracl.h"
  11.  
  12. #define LIMIT1 2000        /* must be int, and > MULT/2 */
  13. #define LIMIT2 100000L     /* may be long */
  14. #define MULT   2310        /* must be int, product of small primes 2.3... */
  15. #define NCURVES 20         /* no. of curves to try */
  16.  
  17. big ak,t,s1,d1,s2,d2;
  18.  
  19. void duplication(sum,diff,n,x,z)
  20. big sum,diff,n,x,z;
  21. { /* double a point on the curve P(x,z)=2.P(x1,z1) */
  22.     mad(sum,sum,sum,n,n,t);   
  23.     mad(diff,diff,diff,n,n,z); 
  24.     mad(t,z,z,n,n,x);           /* x = sum^2.diff^2 */
  25.     subtract(t,z,t);            /* t = sum^2-diff^2 */
  26.     mad(ak,t,z,n,n,z);          /* z = ak*t +diff^2 */
  27.     mad(z,t,t,n,n,z);           /* z = z.t          */
  28. }
  29.  
  30. void addition(xd,zd,sm1,df1,sm2,df2,n,x,z)
  31. big xd,zd,sm1,df1,sm2,df2,n,x,z;
  32. { /* add two points on the curve P(x,z)=P(x1,z1)+P(x2,z2) *
  33.    * given their difference P(xd,zd)                      */
  34.         mad(df2,sm1,sm1,n,n,x);
  35.         mad(df1,sm2,sm2,n,n,z);
  36.         add(z,x,t);
  37.         subtract(z,x,z);
  38.         mad(t,t,t,n,n,x);
  39.         mad(zd,x,x,n,n,x);     /* x = zd.[df1.sm2+sm1.df2]^2 */
  40.         mad(z,z,z,n,n,z);
  41.         mad(xd,z,z,n,n,z);     /* z = xd.[df1.sm2-sm1.df2]^2 */
  42. }
  43.  
  44. void ellipse(x,z,r,n,x1,z1,x2,z2)
  45. big x,z,n,x1,z1,x2,z2;
  46. int r;
  47. { /* calculate point r.P(x,y) on curve */
  48.     int k,rr;
  49.     k=1;
  50.     rr=r;
  51.     copy(x,x1);            
  52.     copy(z,z1);
  53.     add(x1,z1,s1);
  54.     subtract(x1,z1,d1);
  55.     duplication(s1,d1,n,x2,z2);  /* generate 2.P */
  56.     while ((rr/=2)>1) k*=2;
  57.     while (k>0)
  58.     { /* use binary method */
  59.         add(x1,z1,s1);         /* form sums and differences */
  60.         subtract(x1,z1,d1);    /* x+z and x-z for P1 and P2 */
  61.         add(x2,z2,s2);
  62.         subtract(x2,z2,d2);
  63.         if ((r&k)==0)
  64.         { /* double P(x1,z1) mP to 2mP */
  65.             addition(x,z,s1,d1,s2,d2,n,x2,z2);
  66.             duplication(s1,d1,n,x1,z1);
  67.         }
  68.         else
  69.         { /* double P(x2,z2) (m+1)P to (2m+2)P */
  70.             addition(x,z,s1,d1,s2,d2,n,x1,z1);
  71.             duplication(s2,d2,n,x2,z2);
  72.         }    
  73.         k/=2;
  74.     }
  75. }
  76.  
  77. main()
  78. {  /*  factoring program using Lenstras Elliptic Curve method */
  79.     int phase,m,k,nc,iv,pos,btch;
  80.     int *primes;
  81.     long i,p,pa;
  82.     big q,x,z,a,x1,z1,x2,z2,xt,zt,n,fvw;
  83.     static big fu[1+MULT/2];
  84.     static bool cp[1+MULT/2];
  85.     mirsys(50,MAXBASE);
  86.     q=mirvar(0);
  87.     x=mirvar(0);
  88.     z=mirvar(0);
  89.     a=mirvar(0);
  90.     x1=mirvar(0);
  91.     z1=mirvar(0);
  92.     x2=mirvar(0);
  93.     z2=mirvar(0);
  94.     n=mirvar(0);
  95.     t=mirvar(0);
  96.     s1=mirvar(0);
  97.     d1=mirvar(0);
  98.     s2=mirvar(0);
  99.     d2=mirvar(0);
  100.     ak=mirvar(0);
  101.     xt=mirvar(0);
  102.     zt=mirvar(0);
  103.     fvw=mirvar(0);
  104.     primes=gprime(LIMIT1);
  105.     for (m=1;m<=MULT/2;m+=2)
  106.         if (igcd(MULT,m)==1)
  107.         {
  108.             fu[m]=mirvar(0);
  109.             cp[m]=TRUE;
  110.         }
  111.         else cp[m]=FALSE;
  112.     printf("input number to be factored\n");
  113.     cinnum(n,stdin);
  114.     if (prime(n))
  115.     {
  116.         printf("this number is prime!\n");
  117.         exit(0);
  118.     }
  119.     for (nc=1,k=3;k<100;k++)
  120.     { /* try a new curve */
  121.         convert(2,x);         /* generating an elliptic curve */
  122.         convert(1,z);
  123.         convert((2*(k*k-1)),t);
  124.         if (xgcd(t,n,t,t,t)!=1) continue;
  125.         convert((4-k*k),a);
  126.         mad(a,t,t,n,n,t);
  127.         if (size(t)<0) add(t,n,t);
  128.         copy(t,a);        
  129.         if (xgcd(t,n,t,t,t)!=1) continue;
  130.         nc++;
  131.         add(a,t,a);
  132.         divide(a,n,n);
  133.         convert(4,ak);
  134.         xgcd(ak,n,ak,ak,ak);
  135.         incr(a,2,t);
  136.         mad(ak,t,t,n,n,ak);    /* ak = (a+2)/4 */
  137.         phase=1;
  138.         p=0;
  139.         i=0;
  140.         btch=50;
  141.         printf("phase 1 - trying all primes less than %d\n",LIMIT1);
  142.         printf("prime= %8ld",p);
  143.         forever
  144.         { /* main loop */
  145.             if (phase==1)
  146.             {
  147.                 p=primes[i];
  148.                 if (primes[i+1]==0)
  149.                 { /* now change gear */
  150.                     phase=2;
  151.                     printf("\nphase 2 - trying last prime less than %ld\n",
  152.                             LIMIT2);
  153.                     printf("prime= %8ld",p);
  154.                     copy(x,xt);
  155.                     copy(z,zt);
  156.                     add(x,z,s2);
  157.                     subtract(x,z,d2);                    /*   P = (s2,d2) */
  158.                     duplication(s2,d2,n,x,z);
  159.                     add(x,z,s1);
  160.                     subtract(x,z,d1);                    /* 2.P = (s1,d1) */
  161.                     xgcd(z1,n,fu[1],fu[1],fu[1]);        /* fu[1] = x1/z1 */
  162.                     mad(x1,fu[1],x1,n,n,fu[1]);
  163.                     addition(x1,z1,s1,d1,s2,d2,n,x2,z2); /* 3.P = (x2,z2) */
  164.                     for (m=5;m<=MULT/2;m+=2)
  165.                     { /* calculate m.P = (x,y) and store fu[m] = x/y */
  166.                         add(x2,z2,s2);
  167.                         subtract(x2,z2,d2);
  168.                         addition(x1,z1,s2,d2,s1,d1,n,x,z);
  169.                         copy(x2,x1);
  170.                         copy(z2,z1);
  171.                         copy(x,x2);
  172.                         copy(z,z2);
  173.                         if (!cp[m]) continue;
  174.                         copy(z2,fu[m]);
  175.                         xgcd(fu[m],n,fu[m],fu[m],fu[m]);
  176.                         mad(x2,fu[m],x2,n,n,fu[m]);
  177.                     }
  178.                     ellipse(xt,zt,MULT,n,x,z,x2,z2);
  179.                     add(x,z,xt);
  180.                     subtract(x,z,zt);              /* MULT.P = (xt,zt) */
  181.                     iv=p/MULT;
  182.                     if (p%MULT>MULT/2) iv++,p=2*(long)iv*MULT-p;
  183.                     ellipse(x,z,iv,n,x1,z1,x2,z2); /* (x1,z1) = iv.MULT.P */
  184.                     xgcd(z1,n,fvw,fvw,fvw);
  185.                     mad(x1,fvw,x1,n,n,fvw);        /* fvw = x1/z1 */
  186.                     subtract(fvw,fu[p%MULT],q);
  187.                     btch*=10;
  188.                     i++;
  189.                     continue;
  190.                 }
  191.                 pa=p;
  192.                 while ((LIMIT1/p) > pa) pa*=p;
  193.                 convert((int)pa,t);
  194.                 ellipse(x,z,(int)pa,n,x1,z1,x2,z2);
  195.                 copy(x1,x);
  196.                 copy(z1,z);
  197.                 copy(z,q);
  198.             }
  199.             else
  200.             { /* phase 2 - looking for last large prime factor of (p+1+d) */
  201.                 p+=2;
  202.                 pos=p%MULT;
  203.                 if (pos>MULT/2)
  204.                 { /* increment giant step */
  205.                     iv++;
  206.                     p=(long)iv*MULT+1;
  207.                     pos=1;
  208.                     xgcd(z2,n,fvw,fvw,fvw);
  209.                     mad(x2,fvw,x2,n,n,fvw);
  210.                     add(x2,z2,s2);
  211.                     subtract(x2,z2,d2);
  212.                     addition(x1,z1,s2,d2,xt,zt,n,x,z);
  213.                     copy(x2,x1);
  214.                     copy(z2,z1);
  215.                     copy(x,x2);
  216.                     copy(z,z2);
  217.                 }
  218.                 if (!cp[pos]) continue;
  219.                 subtract(fvw,fu[pos],t);
  220.                 mad(q,t,t,n,n,q);         /* batch gcds */
  221.             }
  222.             if (i++%btch==0)
  223.             { /* try for a solution */
  224.                 printf("\b\b\b\b\b\b\b\b%8ld",p);
  225.                 gcd(q,n,t);
  226.                 if (size(t)==1)
  227.                 {
  228.                     if (p>LIMIT2) break;
  229.                     else continue;
  230.                 }
  231.                 if (compare(t,n)==0)
  232.                 {
  233.                     printf("\ndegenerate case");
  234.                     break;
  235.                 }
  236.                 printf("\nfactors are\n");
  237.                 if (prime(t)) printf("prime factor     ");
  238.                 else          printf("composite factor ");
  239.                 cotnum(t,stdout);
  240.                 divide(n,t,n);
  241.                 if (prime(n)) printf("prime factor     ");
  242.                 else          printf("composite factor ");
  243.                 cotnum(n,stdout);
  244.                 exit(0);
  245.             }
  246.         }
  247.         if (nc>NCURVES) break;
  248.         printf("\ntrying a different curve %d\n",nc);
  249.     } 
  250.     printf("\nfailed to factor\n");
  251. }
  252.  
  253.