home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_02 / qsieve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  12.9 KB  |  465 lines
  1. /*
  2.  *   Program to factor big numbers using Pomerance-Silverman-Montgomery  
  3.  *   multiple polynomial quadratic sieve.
  4.  *   See "The Multiple Polynomial Quadratic Sieve", R.D. Silverman,
  5.  *   Math. Comp. Vol. 48, 177, Jan. 1987, pp329-339
  6.  *
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include <math.h> 
  11. #include "miracl.h"
  12.  
  13. #define MEM 240        /* Maximum size of factor base   */
  14. #define MLF 100        /* Max. number of large factors  */
  15. #define PAK 1+MEM/(2*BTS)
  16. #define SSIZE 10000    /* Maximum sieve size            */
  17. #define LOG log        /* if log not available use any other log */
  18.  
  19. big NN,TT,DD,RR,VV,PP,XX,YY,DG,IG,AA,BB;
  20. big x[MEM+1],y[MEM+1],z[MLF+1],w[MLF+1];
  21. int ep[MEM+1],r1[MEM+1],r2[MEM+1],rp[MEM+1],pr[MLF+1];
  22. int b[MEM+1],e[MEM+1],hash[2*MLF+1];
  23. unsigned int EE[MEM+1][PAK],G[MLF+1][PAK];
  24. unsigned char logp[MEM+1],sieve[SSIZE+1];
  25. int mm,jj,nbts,nlp,lp,hmod,hmod2;
  26. bool partial;
  27.  
  28. extern char *malloc();
  29.  
  30. int knuth(ep,N,R,D,T)
  31. int ep[];
  32. big N,R,D,T;
  33. { /* initialization */
  34.     double fks,dp,top;
  35.     bool found;
  36.     int *primes;
  37.     int i,j,mm,bk,nk,m,d,kk,r,p;
  38.     static int K[]={0,1,2,3,5,6,7,10,11,13,14,15,17,0};
  39.     printf("input number to be factored N= \n");
  40.     d=cinnum(N,stdin);
  41.     if (prime(N))
  42.     {
  43.         printf("this number is prime!\n");
  44.         exit(0);
  45.     }
  46.     if (d<10) m=d;
  47.     else m=25;
  48.     if (d>20) m=30*(d-20);
  49.     if (m>MEM) m=MEM;
  50.     primes=gprime((-2)*(m+5));
  51.     top=(-10.0e0);
  52.     found=FALSE;
  53.     nk=0;
  54.     bk=0;
  55.     ep[0]=1;
  56.     ep[1]=2;
  57.     do
  58.     { /* search for best Knuth-Schroepel multiplier */
  59.         mm=m;
  60.         kk=K[++nk];
  61.         if (kk==0)
  62.         { /* finished */
  63.             kk=K[bk];
  64.             found=TRUE;
  65.         }
  66.         premult(N,kk,D);
  67.         if (root(D,2,R))
  68.         {
  69.             printf("%dN is a perfect square!\n",kk);
  70.             continue;
  71.         }
  72.         fks=LOG(2.0e0)/(2.0e0);
  73.         r=subdiv(D,8,T);
  74.         if (r==1) fks*=(4.0e0);
  75.         if (r==5) fks*=(2.0e0);
  76.         fks-=LOG((double)kk)/(2.0e0);
  77.         i=0;
  78.         j=1;
  79.         while (j<m)
  80.         { /* select small primes */
  81.             p=primes[++i];
  82.             if (p==0)
  83.             {
  84.                 mm=j;
  85.                 break;
  86.             }
  87.             r=subdiv(D,p,T);
  88.             if (spmd(r,(p-1)/2,p)<=1)
  89.             { /* use only if Jacobi symbol = 0 or 1 */
  90.                 ep[++j]=p;
  91.                 dp=(double)p;
  92.                 if (kk%p==0) fks+=LOG(dp)/dp;
  93.                 else         fks+=2*LOG(dp)/(dp-1.0e0);
  94.             }
  95.         }
  96.         if (fks>top)
  97.         { /* find biggest fks */
  98.             top=fks;
  99.             bk=nk;
  100.         }
  101.     } while (!found);
  102.     printf("using multiplier k= %d\n",kk);
  103.     printf("and %d small primes as factor base\n",mm);
  104.     free(primes);
  105.     return mm;
  106. }
  107.  
  108. bool factored(lptr,T)
  109. big T;
  110. long lptr;
  111. { /* factor quadratic residue */
  112.     bool facted;
  113.     int i,j,r,st;  
  114.     partial=FALSE;
  115.     facted=FALSE;
  116.     for (j=1;j<=mm;j++)
  117.     { /* now attempt complete factorisation of T */
  118.         r=lptr%ep[j];
  119.         if (r<0) r+=ep[j];
  120.         if (r!=r1[j] && r!=r2[j]) continue;
  121.         while (subdiv(T,ep[j],XX)==0)
  122.         { /* cast out ep[j] */
  123.             e[j]++;
  124.             copy(XX,T);
  125.         }
  126.         st=size(T);
  127.         if (st==1)
  128.         {
  129.            facted=TRUE;
  130.            break;
  131.         }
  132.         if (size(XX)<=ep[j])
  133.         { /* st is prime < ep[mm]^2 */
  134.             if (st>=TOOBIG || (st/ep[mm])>=(1+MLF/50)) break;
  135.             if (st<=ep[mm])
  136.                 for (i=j;i<=mm;i++)
  137.                 if (st==ep[i])
  138.                 {
  139.                     e[i]++;
  140.                     facted=TRUE;
  141.                     break;
  142.                 }
  143.             if (facted) break;
  144.             lp=st;  /* factored with large prime */
  145.             partial=TRUE;
  146.             facted=TRUE;
  147.             break;
  148.         }
  149.     }
  150.     return facted;
  151. }
  152.  
  153. bool gotcha()
  154. { /* use new factorisation */
  155.     int r,j,i,k,n,rb,had,hp;
  156.     bool found;
  157.     found=TRUE;
  158.     if (partial)
  159.     { /* check partial factorisation for usefulness */
  160.         had=lp%hmod;
  161.         forever
  162.         { /* hash search for matching large prime */
  163.             hp=hash[had];
  164.             if (hp<0)
  165.             { /* failed to find match */
  166.                 found=FALSE;
  167.                 break;
  168.             }
  169.             if (pr[hp]==lp) break; /* hash hit! */
  170.             had=(had+(hmod2-lp%hmod2))%hmod;
  171.         }
  172.         if (!found && nlp>=MLF) return FALSE;
  173.     }
  174.     copy(PP,XX);
  175.     convert(1,YY);
  176.     for (k=1;k<=mm;k++)
  177.     { /* build up square part in YY  *
  178.        * reducing e[k] to 0s and 1s */
  179.         if (e[k]<2) continue;
  180.         r=e[k]/2;
  181.         e[k]%=2;
  182.         convert(ep[k],TT);
  183.         power(TT,r,TT,TT);
  184.         multiply(TT,YY,YY);
  185.     }
  186. /* debug only - print out details of each factorization 
  187.     cotnum(XX,stdout);
  188.     cotnum(YY,stdout);
  189.     if (e[0]==1) printf("-1");
  190.     else printf("1");
  191.     for (k=1;k<=mm;k++)
  192.     {
  193.         if (e[k]==0) continue;
  194.         printf(".%d",ep[k]);
  195.     }
  196.     if (partial) printf(".%d\n",lp);
  197.     else printf("\n");
  198. */
  199.     if (partial)
  200.     { /* factored with large prime */
  201.         if (!found)
  202.         { /* store new partial factorization */
  203.             hash[had]=nlp;
  204.             pr[nlp]=lp;
  205.             copy(XX,z[nlp]);
  206.             copy(YY,w[nlp]);
  207.             for (n=0,rb=0,j=0;j<=mm;j++)
  208.             {
  209.                 G[nlp][n]|=((e[j]&1)<<rb);
  210.                 if (++rb==nbts) n++,rb=0;
  211.             }
  212.             nlp++;
  213.         }
  214.         if (found)
  215.         { /* match found so use as factorization */
  216.             printf("\b\b\b\b\b*");
  217.             mad(XX,z[hp],XX,NN,NN,XX);
  218.             mad(YY,w[hp],YY,NN,NN,YY);
  219.             for (n=0,rb=0,j=0;j<=mm;j++)
  220.             {
  221.                 e[j]+=((G[hp][n]>>rb)&1);
  222.                 if (e[j]==2)
  223.                 {
  224.                     premult(YY,ep[j],YY);
  225.                     divide(YY,NN,NN);
  226.                     e[j]=0;
  227.                 }
  228.                 if (++rb==nbts) n++,rb=0;
  229.             }
  230.             premult(YY,lp,YY);
  231.             divide(YY,NN,NN);
  232.         }
  233.     }
  234.     else printf("\b\b\b\b\b ");
  235.     if (found)
  236.     {
  237.         for (k=mm;k>=0;k--)
  238.         { /* use new factorization in search for solution */
  239.             if (e[k]%2==0) continue;
  240.             if (b[k]<0)
  241.             { /* no solution this time */
  242.                 found=FALSE;
  243.                 break;
  244.             }
  245.             i=b[k];
  246.             mad(XX,x[i],XX,NN,NN,XX);    /* This is very inefficient -  */
  247.             mad(YY,y[i],YY,NN,NN,YY);    /* There must be a better way! */
  248.             for (n=0,rb=0,j=0;j<=mm;j++)
  249.             { /* Gaussian elimination */
  250.                 e[j]+=((EE[i][n]>>rb)&1);
  251.                 if (++rb==nbts) n++,rb=0;
  252.             }
  253.         }
  254.         for (j=0;j<=mm;j++)
  255.         { /* update YY */
  256.             if (e[j]<2) continue;
  257.             convert(ep[j],TT);
  258.             power(TT,e[j]/2,NN,TT);
  259.             mad(YY,TT,YY,NN,NN,YY);
  260.         }
  261.         if (!found)
  262.         { /* store details in EE, x and y for later */
  263.             b[k]=jj;
  264.             copy(XX,x[jj]);
  265.             copy(YY,y[jj]);
  266.             for (n=0,rb=0,j=0;j<=mm;j++)
  267.             {
  268.                 EE[jj][n]|=((e[j]&1)<<rb);
  269.                 if (++rb==nbts) n++,rb=0;
  270.             }
  271.             jj++;
  272.             printf("%4d",jj);
  273.         }
  274.     }
  275.     if (found)
  276.     { /* check for false alarm */
  277.         printf("\ntrying...\n");
  278.         add(XX,YY,TT);
  279.         if (compare(XX,YY)==0 || compare(TT,NN)==0) found=FALSE;
  280.         if (!found) printf("working... %4d",jj);
  281.     }
  282.     return found;
  283. }
  284.  
  285. void initv()
  286. { /* initialize big numbers and arrays */
  287.     int i,j;
  288.     NN=mirvar(0); 
  289.     TT=mirvar(0);
  290.     DD=mirvar(0);
  291.     RR=mirvar(0);
  292.     VV=mirvar(0);
  293.     PP=mirvar(0);
  294.     XX=mirvar(0);
  295.     YY=mirvar(0);
  296.     DG=mirvar(0);
  297.     IG=mirvar(0);
  298.     AA=mirvar(0);
  299.     BB=mirvar(0);
  300.     for (i=0;i<=MEM;i++)
  301.     {
  302.         x[i]=mirvar(0);
  303.         y[i]=mirvar(0);
  304.     }
  305.     for (i=0;i<=MLF;i++)
  306.     {
  307.         z[i]=mirvar(0);
  308.         w[i]=mirvar(0);
  309.     }
  310.  
  311.     nbts=8*sizeof(int);
  312.     mm=knuth(ep,NN,RR,DD,TT);
  313.  
  314.     for (i=0;i<=mm;i++)
  315.     { 
  316.         b[i]=(-1);
  317.         for (j=0;j<=mm/nbts;j++)
  318.             EE[i][j]=0;
  319.     }
  320.     for (i=0;i<=MLF;i++)
  321.         for (j=0;j<=mm/nbts;j++)
  322.             G[i][j]=0;
  323. }
  324.  
  325. main()
  326. { /* factoring via quadratic sieve */
  327.     int i,j,S,r,s1,s2,s,NS,logm,a,ptr,threshold,isz;
  328.     unsigned int  *isp;
  329.     long M,la,lptr;
  330.     mirsys(-30,MAXBASE);
  331.     initv();
  332.     isz=SSIZE/sizeof(int);
  333.     isp=(unsigned int *)sieve;  /* equivalence sieve to int array */
  334.  
  335.     hmod=2*MLF+1;               /* set up hash table */
  336.     convert(hmod,TT);
  337.     while (!prime(TT)) decr(TT,2,TT);
  338.     hmod=size(TT);
  339.     hmod2=hmod-2;
  340.     for (i=0;i<hmod;i++) hash[i]=(-1);
  341.  
  342.     M=50*(long)mm;
  343.     NS=M/SSIZE;
  344.     if (M%SSIZE!=0) NS++;
  345.     M=SSIZE*(long)NS;
  346.     logm=0;
  347.     la=M;
  348.     while ((la/=2)>0) logm++;   /* logm = log(M) */
  349.  
  350.     rp[0]=logp[0]=0;
  351.     for (i=1;i<=mm;i++)
  352.     { /* find root mod each prime, and approx log of each prime */
  353.         r=subdiv(DD,ep[i],TT);
  354.         rp[i]=sqrmp(r,ep[i]);
  355.         logp[i]=0;
  356.         r=ep[i];
  357.         while((r/=2)>0) logp[i]++;
  358.     }
  359.     r=subdiv(DD,8,TT);    /* take special care of 2 */
  360.     if (r==5) logp[1]++;
  361.     if (r==1) logp[1]+=2;
  362.  
  363.     threshold=logm+logb2(RR)-2*logp[mm];
  364.  
  365.     jj=0;
  366.     nlp=0;
  367.     premult(DD,2,DG);
  368.     root(DG,2,DG);
  369.     lconv(M,TT);
  370.     divide(DG,TT,DG);
  371.     root(DG,2,DG);
  372.     if (subdiv(DG,2,TT)==0) incr(DG,1,DG);
  373.     if (subdiv(DG,4,TT)==1) incr(DG,2,DG);
  374.     printf("working...    0");
  375.     forever
  376.     { /* try a new polynomial */
  377.         r=NTRY;
  378.         NTRY=1;         /* speed up search for prime */
  379.         do
  380.         { /* looking for suitable prime DG = 3 mod 4 */
  381.             do incr(DG,4,DG); while(!prime(DG));
  382.             decr(DG,1,TT);
  383.             subdiv(TT,2,TT);
  384.             powmod(DD,TT,DG,TT);  /* check D is quad residue */
  385.         } while (size(TT)!=1);
  386.         NTRY=r;
  387.         incr(DG,1,TT);
  388.         subdiv(TT,4,TT);
  389.         powmod(DD,TT,DG,BB);
  390.         negate(DD,TT);
  391.         mad(BB,BB,TT,DG,TT,TT);
  392.         negate(TT,TT);
  393.         premult(BB,2,AA);
  394.         xgcd(AA,DG,AA,AA,AA);
  395.         mad(AA,TT,TT,DG,DG,AA);
  396.         multiply(AA,DG,TT);
  397.         add(BB,TT,BB);        /* BB^2 = DD mod DG^2 */
  398.         multiply(DG,DG,AA);   /* AA = DG*DG         */
  399.         xgcd(DG,DD,IG,IG,IG); /* IG = 1/DG mod DD  */
  400.  
  401.         r1[0]=r2[0]=0;
  402.         for (i=1;i<=mm;i++) 
  403.         { /* find roots of quadratic mod each prime */
  404.             s=subdiv(BB,ep[i],TT);
  405.             r=subdiv(AA,ep[i],TT);
  406.             r=inverse(r,ep[i]);     /* r = 1/AA mod p */
  407.             s1=(-s+rp[i]);
  408.             s2=(-s-rp[i]);
  409.             r1[i]=smul(s1,r,ep[i]);
  410.             r2[i]=smul(s2,r,ep[i]);
  411.         }
  412.  
  413.         for (ptr=(-NS);ptr<NS;ptr++)
  414.         { /* sieve over next period */
  415.             la=(long)ptr*SSIZE;
  416.             for (i=0;i<isz;i++) isp[i]=0;
  417.             for (i=1;i<=mm;i++)
  418.             { /* sieving with each prime */
  419.                 r=la%ep[i];
  420.                 s1=(r1[i]-r)%ep[i];
  421.                 if (s1<0) s1+=ep[i];
  422.                 s2=(r2[i]-r)%ep[i];
  423.                 if (s2<0) s2+=ep[i];
  424.                 for (j=s1;j<SSIZE;j+=ep[i]) sieve[j]+=logp[i];
  425.                 if (s1==s2) continue;
  426.                 for (j=s2;j<SSIZE;j+=ep[i]) sieve[j]+=logp[i];
  427.             }
  428.  
  429.             for (a=0;a<SSIZE;a++)
  430.             { /* main loop - look for fully factored residues */
  431.                 if (sieve[a]<threshold) continue;
  432.                 lptr=la+a;
  433.                 lconv(lptr,TT);
  434.                 S=0;
  435.                 multiply(AA,TT,TT);           /* TT = AAx + BB      */
  436.                 add(TT,BB,TT);
  437.                 mad(TT,IG,TT,DD,DD,PP);       /* PP = (AAx + BB)/G  */
  438.                 if (size(PP)<0) add(PP,DD,PP);
  439.                 mad(PP,PP,PP,DD,DD,VV);       /* VV = PP^2 mod kN  */
  440.                 absol(TT,TT);
  441.                 if (compare(TT,RR)<0) S=1;    /* check for -ve VV */
  442.                 if (S==1) subtract(DD,VV,VV);
  443.                 copy(VV,TT);
  444.                 e[0]=S;
  445.                 for (i=1;i<=mm;i++) e[i]=0;
  446.                 if (!factored(lptr,TT)) continue;
  447.                 if (gotcha())
  448.                 { /* factors found! */
  449.                     gcd(TT,NN,PP);
  450.                     printf("factors are\n");
  451.                     if (prime(PP)) printf("prime factor     ");
  452.                     else          printf("composite factor ");
  453.                     cotnum(PP,stdout);
  454.                     divide(NN,PP,NN);
  455.                     if (prime(NN)) printf("prime factor     ");
  456.                     else          printf("composite factor ");
  457.                     cotnum(NN,stdout);
  458.                     exit(0);
  459.                 }
  460.             }
  461.         }
  462.     }
  463. }
  464.  
  465.