quadraticsieve.java

factorization.zip

package javax.math.factorization;

/**
 * <p>Title: Factorization Library</p>
 * <p>Description: </p>
 * <p>Copyright: Copyright (c) 2004</p>
 * <p>Company: </p>
 * @author Vladimir Silva
 * @version 1.0
 */
import java.math.BigInteger;

/**
 * Sorry about all this awful spagethi code :-<
 * I got it from: Dario Alejandro Alpern (Buenos Aires - Argentina)
 * Last updated October 16th, 2003. See http://www.alpertron.com.ar/ECM.HTM
 * I am not to blame for its crappy design and implementation, but....
 * what the hell. it works!
 * I modified it though, to make it less horrible (extracted the pieces i needed)
 * If you are a brave soul,  try to port this stuff to an object oriented, well-designed
 * structure. I tried, but just gave up....
 */
public class QuadraticSieve implements Constants
{
      static long polynomialsSieved ;
      static long trialDivisions ;
      static long smoothsFound;
      static long totalPartials;
      static long partialsFound;

      static final long TestNbr2[] = new long[MP_INITIAL_SIZE];

      /**
       *
       */


  static BigInteger FactoringSIQS(BigInteger NbrToFactor) {

    BigInteger RealNbrToFactor = null;
    long FactorBase;
    int SieveLimit;
    int Expon;
    int NbrPrimes, NbrPrimes2;
    int s, F, F1, F2, F3, F4, X1, X2;
    long currentPrime;
    int NbrMod;
    long modsqrt[];
    long prime[];
    byte logar[];
    long ainv[];
    int Bainv2[][];
    int soln1[];
    int soln2[];
    long afact[];
    int aindex[];
    int amodq[];
    byte SieveArray[];
    long Power, Power2, Square, SqrRootMod, fact;
    int index;
    int index2;
    int PolynomialIndex;
    long D, E, Q, V, W, X, Y, Z, T1, V1, W1, Y1;
    long t1, t2, t3;
    double Temp, Prod;
    byte G, H, Q1, Q2;
    int NbrPolynomials;
    int PolynomialSetNbr = 0;
    int pp=0;
    double dp,fks,top;
    int i,j,k,bk,nk,kk,r;
    int K[]={0,1,2,3,5,6,7,10,11,13,14,15,17,0};
    long seed = 0;
    int span,min;
    int mask;
    int vectExpParity[];
    int matrixV[];
    int matrixB[][];
    int rowMatrixB[];
    long vectLeftHandSide[][];
    int nbrColumns, expParity;
    int matrixPartial[][];
    int matrixPartialHashIndex[];
    int prev;
    int nbrPartials = 0;
    int smallPrimeLowerLimit, smallPrimeUpperLimit;
    int firstLimit;
    int secondLimit;
    int thirdLimit;
    double dNumberToFactor;
    String SIQSInfoText;
    long Divid, Divisor, Rem;
    long startTime;
    long biR0=0, biR1=0, biR2=0, biR3=0, biR4=0, biR5=0;
    long biR6=0, biR7=0, biR8=0, biR9=0, biR10=0;
    boolean cond = false;
    int S1, S2, G0, G1, G2, G3;
    boolean polyadd;

    Temp = (double)NbrToFactor.bitLength();
    NbrPrimes = (int)Math.exp(Math.sqrt(Temp*Math.log(Temp))*0.255);
    SieveLimit = (int)Math.exp(8.5 + 0.010*Temp);
    if (SieveLimit > 30000) {SieveLimit = 30000;}
    s = NbrToFactor.bitLength()/28 + 1;
    prime = new long[NbrPrimes+3];
    modsqrt = new long[NbrPrimes+3];
    logar = new byte[NbrPrimes+3];
    ainv = new long[NbrPrimes+3];
    Bainv2 = new int[s][NbrPrimes+3];
    soln1 = new int[NbrPrimes+3];
    soln2 = new int[NbrPrimes+3];
    afact = new long[s];
    aindex = new int[s];
    amodq = new int[s];
    rowMatrixB = new int[200];
    NbrPolynomials = (1 << (s-1)) - 1;

     BigNumOps.BigNbrToBigInt(NbrToFactor);

    Factorizer.TestNbr[Factorizer.NumberLength] = 0;
    for (i=Factorizer.NumberLength; i>=0; i--) {
      TestNbr2[i] = Factorizer.TestNbr[i];
    }
    Factorizer.NumberLength++;
    matrixPartialHashIndex = new int[1024];
    for (i=matrixPartialHashIndex.length-1; i>=0; i--) {
      matrixPartialHashIndex[i] = -1;
    }

    Util.verbose("Factoring ("+NbrToFactor.toString().length()+" digits)");

    Util.verbose("\n\nSIQS parameters: " + NbrPrimes +
      " primes, sieve limit: " + SieveLimit + "\nSearching for Knuth-Schroeppel multiplier...");

    /************************/
    /* Compute startup data */
    /************************/

    top=-10.0e0;
    nk=0;
    bk=0;
    prime[0]=1;
    prime[1]=2;
    currentPrime = 3;
    while (true) { // search for best Knuth-Schroepel multiplier
      ShowSIQSStatus(0, 1, 0);     // Show time in applet
      kk=K[++nk];
      if (kk==0)
      { // finished
        kk=K[bk];
        BigNumOps.MultBigNbrByLong(TestNbr2, kk, Factorizer.TestNbr);
        if (Factorizer.TestNbr[Factorizer.NumberLength-1] != 0 || Factorizer.TestNbr[Factorizer.NumberLength-2] > Constants.Mi) {
          Factorizer.TestNbr[Factorizer.NumberLength] = 0;
          Factorizer.NumberLength++;
        }
        break;
      }
      BigNumOps.MultBigNbrByLong(TestNbr2, kk, Factorizer.TestNbr);
      RealNbrToFactor = NbrToFactor.multiply(BigInteger.valueOf(kk));
      fks=0.34657359;          //  (ln 2)/2
      r = RealNbrToFactor.and(BigInteger.valueOf(7)).intValue();
      if (r==1) fks*=(4.0e0);
      if (r==5) fks*=(2.0e0);
      fks-=Math.log((double)kk)/(2.0e0);
      j=2;
      currentPrime = 3;
      while (j<NbrPrimes) { // select small primes
        NbrMod = (int)BigNumOps.RemDivBigNbrByLong(Factorizer.TestNbr, currentPrime);
        if (BigNumOps.modPow(NbrMod, (currentPrime - 1)/2, currentPrime) == 1) {
          // use only if Jacobi symbol = 0 or 1
          j++;
          dp=(double)currentPrime;
          if (kk%currentPrime==0) {
            fks+=Math.log(dp)/dp;
          }
          else {
            fks+=2*Math.log(dp)/(dp-1.0e0);
          }
        }
        calculate_new_prime1:
          do {
          currentPrime += 2;
          for (Q=3; Q*Q<=currentPrime; Q+=2) {// Check if currentPrime is prime
            if (currentPrime%Q == 0) {
              continue calculate_new_prime1;
            }
          }
          break;                  // Prime found
        } while (true);
      }
      if (fks>top) { // find biggest fks
        top=fks;
        bk=nk;
      }
    }          // end while
    matrixPartial = new int[NbrPrimes*8][Factorizer.NumberLength+2];
    Factorizer.dN = (double)Factorizer.TestNbr[Factorizer.NumberLength-2];
    if (Factorizer.NumberLength > 1) {
      Factorizer.dN += (double)Factorizer.TestNbr[Factorizer.NumberLength-3]/Constants.dDosALa32;
    }
    if (Factorizer.NumberLength > 2) {
      Factorizer.dN += (double)Factorizer.TestNbr[Factorizer.NumberLength-4]/Constants.dDosALa64;
    }
    FactorBase = currentPrime;
    matrixB = new int[(int)((float)NbrPrimes*1.05 + 40)][];
    vectLeftHandSide = new long[matrixB.length][];
    vectExpParity = new int[matrixB.length];
    modsqrt[1] = NbrToFactor.testBit(0)?1:0;
    logar[1] = 1;
    j=2;
    currentPrime = 3;
    while (j<NbrPrimes) { // select small primes
      NbrMod = (int)BigNumOps.RemDivBigNbrByLong(Factorizer.TestNbr, currentPrime);
      if (BigNumOps.modPow(NbrMod, (currentPrime - 1)/2, currentPrime) == 1) {
        // use only if Jacobi symbol = 0 or 1
        prime[j]=currentPrime;
        NbrMod = (int)BigNumOps.RemDivBigNbrByLong(Factorizer.TestNbr, currentPrime);
        if (currentPrime%4 == 3) {
          SqrRootMod = BigNumOps.modPow(NbrMod, (currentPrime+1)/4, currentPrime);
        }
        else {
          if (currentPrime%8 == 5) {
            SqrRootMod = BigNumOps.modPow(NbrMod*2, (currentPrime-5)/8, currentPrime);
            SqrRootMod = ((((2*NbrMod*SqrRootMod % currentPrime)*SqrRootMod-1)
                  % currentPrime) * NbrMod % currentPrime) * SqrRootMod
              % currentPrime;
          }
          else {                // p = 1 (mod 8)
            Q = currentPrime - 1;
            E = 0;
            Power2 = 1;
            do {
              E++;
              Q /= 2;
              Power2 *= 2;
            } while ((Q & 1) == 0);        // E >= 3
            Power2 /= 2;
            X = 1;
            do {
              X++;
              Z = BigNumOps.modPow(X, Q, currentPrime);
            } while (BigNumOps.modPow(Z, Power2, currentPrime) == 1);
            Y = Z;
            X = BigNumOps.modPow(NbrMod, (Q-1)/2, currentPrime);
            V = NbrMod*X % currentPrime;
            W = V*X % currentPrime;
            while (W != 1) {
              T1=0; D=W;
              while (D != 1) {
                D = D*D % currentPrime;
                T1++;
              }
              D = BigNumOps.modPow(Y, 1 << (E-T1-1), currentPrime);
              Y1 = D*D % currentPrime;
              E = T1;
              V1 = V*D % currentPrime;
              W1 = W*Y1 % currentPrime;
              Y = Y1; V = V1; W = W1;
            }    // end while
            SqrRootMod = V;
          }      // end if
        }        // end if
        modsqrt[j] = SqrRootMod;
        logar[j] = (byte)(Math.round(Math.log((double)currentPrime)/Math.log(2)));
        j++;
      }          // end while
      calculate_new_prime2:
        do {
        currentPrime += 2;
        for (Q=3; Q*Q<=currentPrime; Q+=2) {// Check if currentPrime is prime
          if (currentPrime%Q == 0) {
            continue calculate_new_prime2;
          }
        }
        break;                  // Prime found
      } while (true);
    }                         // End while

    FactorBase = currentPrime;
    SieveArray = new byte[2*SieveLimit+1 > (int)FactorBase? 2*SieveLimit+5000: (int)FactorBase];
    dNumberToFactor = NbrToFactor.doubleValue();

    Util.verbose("Multiplier: "+kk + ", factor base: "+FactorBase);
//    lowerTextArea.setText(SIQSInfoText);

    firstLimit = 2;
    for (j=2; j<NbrPrimes; j++) {
      firstLimit *= (int)prime[j];
      if (firstLimit > 2*SieveLimit) {break;}
      }
    smallPrimeLowerLimit = j;
    smallPrimeUpperLimit = j+1;
    int logarsmall = logar[j+1];
    G = (byte)(Math.log(Math.sqrt(dNumberToFactor) * SieveLimit/(FactorBase*64)/prime[j+1])/Math.log(2));
    firstLimit = (int)(Math.log(dNumberToFactor)/3);
    for (secondLimit=firstLimit; secondLimit<NbrPrimes; secondLimit++) {
      if (prime[secondLimit]*2 > SieveLimit) {
        break;
        }
      }
    for (thirdLimit=secondLimit; thirdLimit<NbrPrimes; thirdLimit++) {
      if (prime[thirdLimit] > 2*SieveLimit) {
        break;
        }
      }
    NbrPrimes2 = NbrPrimes - 4;
    startTime=System.currentTimeMillis();    // Sieve start time in milliseconds.
    sieve_stage:
      do {
      //
      // Initialization stage for first polynomial
      //
      PolynomialIndex = 1;
      Prod = Math.sqrt(2*dNumberToFactor) / (double)SieveLimit;
      fact = (long)Math.pow(Prod,1/(float)s);
      for (i = 2; ; i++) {
        if (prime[i] > fact) {
          break;
        }
      }
      span = (int)(NbrPrimes / s / s / 2);
      if (NbrPrimes < 500) {span *= 2;}
      min = i - span/2;
      for (index = 0; index < s; index++) {
        do {
          seed = (1141592621*seed+321435) & Constants.MaxUInt;
          i = (int)(((seed*span) >> 32)+min);
          for (index2 = 0; index2 < index; index2++) {
            if (aindex[index2] == i || aindex[index2] == i+1) {
              break;
            }
          }
        } while (index2 < index);
        afact[index] = prime[i];
        aindex[index] = i;