primechecker.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;

/**
 * I am truly 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 PrimeChecker implements Constants
{

  // Prime checking routine
  // Return codes: 0 = Number is prime.
  //               1 = Number is composite.
  static int ProbabilisticPrimeTest(BigInteger N) {
    long Base,Q;
    int baseNbr, nbrBases, exp, index, j, k;
    long mask;

//    System.out.println("Starting Rabin probabilistic prime check routine.");
    BigNumOps.BigNbrToBigInt(N);


    exp = N.subtract(Constants.BigInt1).getLowestSetBit();

    MontgomeryParams.GetMontgomeryParams();

    Base = 1;
    nbrBases = N.bitLength()/2;
    for (baseNbr=0; baseNbr<nbrBases; baseNbr++) {
      if (Base<3) {
        Base++;
      }
      else {
        calculate_new_prime4:
          do {
          Base+=2;
          for (Q=3; Q*Q<=Base; Q+=2) {     /* Check if Base is prime */
            if (Base%Q == 0) {
              continue calculate_new_prime4;   /* Composite */
            }
          }
          break;                      /* Prime found */
        } while (true);
      }                 /* end if */

      Util.verbose("Rabin probabilistic prime check routine\n\nBase used: "+Base+" ("+baseNbr*100/nbrBases+"%)");

      //System.arraycopy(MontgomeryMultR1, 0, biN, 0, Factorizer.NumberLength);
      System.arraycopy(MontgomeryParams.getR1(), 0, biN, 0, Factorizer.NumberLength);
      index = Factorizer.NumberLength-1;

      mask = Constants.DosALa31;

      for (k=Factorizer.NumberLength*32; k>exp; k--) {
        // MontgomeryMult(biN, biN, biT);
        MontgomeryParams.MontgomeryMult(biN, biN, biT);

        if ((Factorizer.TestNbr[index] & mask) != 0) {
          BigNumOps.MultBigNbrByLongModN(biT, Base, biT);
        }
        System.arraycopy(biT, 0, biN, 0, Factorizer.NumberLength);
        mask /= 2;
        if (mask == 0) {
          index--;
          mask = Constants.DosALa31;
        }
      }
      for (j=0; j<Factorizer.NumberLength; j++) {
        //if (biN[j] != MontgomeryMultR1[j]) {break;}
        if (biN[j] != MontgomeryParams.getR1()[j]) {break;}
      }
      if (j==Factorizer.NumberLength) {continue;}  /* Probable prime, go to next base */
      for (k=0; k<exp; k++) {
        for (j=0; j<Factorizer.NumberLength; j++) {
          //if (biN[j] != MontgomeryMultR1[j]) {break;}
          if (biN[j] != MontgomeryParams.getR1()[j]) {break;}
        }

        if (j==Factorizer.NumberLength) {return 1;}  /* Composite number */

        //AddBigNbr(biN, MontgomeryMultR1, biT);
        BigNumOps.AddBigNbr(biN, MontgomeryParams.getR1(), biT);

        for (j=0; j<Factorizer.NumberLength; j++) {
          if (biT[j] != Factorizer.TestNbr[j]) {break;}
        }
        if (j==Factorizer.NumberLength) {break;}  /* Probable prime, go to next base */

        //MontgomeryMult(biN, biN, biT);
        MontgomeryParams.MontgomeryMult(biN, biN, biT);
        System.arraycopy(biT, 0, biN, 0, Factorizer.NumberLength);
      }
      if (k==exp) {
        return 1;                      /* Composite number */
      }
    }
    return 0;                        /* Indicate probable prime */
  }


  // Prime checking routine
  // Return codes: 0 = Number is prime.
  //               1 = Number is composite.
  static int AprtCle(BigInteger N) {
	int i, j, G, H, I, J, K, P, Q, T, U, W, X;
	int IV, InvX, LEVELnow, NP, PK, PL, PM, SW, VK, TestedQs, TestingQs;
	int QQ, T1, T3, U1, U3, V1, V3;
	int LengthN, LengthS;
	long Mask;
	double dS;

	String primalityString = "";

//	System.out.println("Starting Prime Check routine.");

	BigNumOps.BigNbrToBigInt(N);

	MontgomeryParams.GetMontgomeryParams();
        long[] MontgomeryMultR1 = MontgomeryParams.getR1();

//	if (Computing3Squares == false) {
	      Util.verbose("Testing primality of " + N  + "(" + N.toString().length() + " digits)" +
                                 "\nAPRT-CLE progress: ");
//	}

	j = PK = PL = PM = 0;
	for (I = 0; I < Factorizer.NumberLength; I++) {
	      biS[I] = 0;
	      for (J = 0; J < MAX_PW; J++) {
		    aiJX[J][I] = 0;
	      }
	}


    GetPrimes2Test:
      for (i=0; i< MAX_LEVEL; i++) {
      biS[0] = 2;
      for (I=1; I<Factorizer.NumberLength; I++) {
        biS[I] = 0;
      }
      for (j=0; j<aiNQ[i]; j++) {
        Q = aiQ[j];
        U = aiT[i]*Q;
        do {
          U /= Q;
          BigNumOps.MultBigNbrByLong(biS, Q, biS);
        } while (U % Q == 0);

        // Exit loop if S^2 > N.

        if (BigNumOps.CompareSquare(biS, Factorizer.TestNbr) > 0) {break GetPrimes2Test;}
      }      // End for j
    }        // End for i


    if (i == MAX_LEVEL) {    // too big
      return ProbabilisticPrimeTest(N);
    }
    LEVELnow = i;
    TestingQs = j;
    T = aiT[LEVELnow];
    NP = aiNP[LEVELnow];


    MainStart:
      do {
      for (i=0; i<NP; i++) {
        P = aiP[i];
        SW = TestedQs = 0;
        Q = W = (int)BigNumOps.BigNbrModLong(Factorizer.TestNbr, P*P);
        for (J=P-2; J>0; J--) {
          W = (W*Q)%(P*P);
        }
        if (P>2 && W!=1) {
          SW = 1;
        }
        do {
          for (j=TestedQs; j<=TestingQs; j++) {
            Q = aiQ[j]-1; G = aiG[j];
            K = 0;
            while (Q%P == 0) {
              K++;
              Q /= P;
            }
            Q = aiQ[j];
            if (K == 0) {continue;}

            //if (Computing3Squares == false) {
              Util.verbose(primalityString + "P = "+P+",  Q = "+Q+"  ("+(i*(TestingQs+1)+j)*100/(NP*(TestingQs+1))+"%)");
            //}

            PM = 1;
            for (I=1; I<K; I++) {
              PM = PM*P;
            }
            PL = (P-1)*PM;
            PK = P*PM;
            J = 1;
            for (I=1; I<Q; I++) {
              J = J*G % Q;
              aiIndx[J] = I;
            }
            J = 1;
            for (I=1; I<=Q-2; I++) {
              J = J*G % Q;
              aiF[I] = aiIndx[(Q+1-J) % Q];
            }
            for (I=0; I<PK; I++) {
              for (J=0; J<Factorizer.NumberLength; J++) {
                aiJ0[I][J] = aiJ1[I][J] = 0;
              }
            }
            if (P>2) {
              BigNumOps.JacobiSum(1,1,P,PK,PL,PM,Q);
            }
            else {
              if (K != 1) {
                BigNumOps.JacobiSum(1,1,P,PK,PL,PM,Q);
                for (I=0; I<PK; I++) {
                  for (J=0; J<Factorizer.NumberLength; J++) {
                    aiJW[I][J] = 0;
                  }
                }
                if (K != 2) {
                  for (I=0; I<PM; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJW[I][J] = aiJ0[I][J];
                    }
                  }
                  BigNumOps.JacobiSum(2,1,P,PK,PL,PM,Q);
                  for (I=0; I<PM; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJS[I][J] = aiJ0[I][J];
                    }
                  }
                  BigNumOps.JS_JW(PK, PL, PM, P);
                  BigNumOps.NormalizeJS(PK, PL, PM, P);
                  for (I=0; I<PM; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJ1[I][J] = aiJS[I][J];
                    }
                  }
                  BigNumOps.JacobiSum(3 << (K-3),1 << (K-3),P,PK,PL,PM,Q);
                  for (J=0; J<Factorizer.NumberLength; J++) {
                    for (I=0; I<PK; I++) {
                      aiJW[I][J] = 0;
                    }
                    for (I=0; I<PM; I++) {
                      aiJS[I][J] = aiJ0[I][J];
                    }
                  }
                  BigNumOps.JS_2(PK, PL, PM, P);
                  BigNumOps.NormalizeJS(PK, PL, PM, P);
                  for (I=0; I<PM; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJ2[I][J] = aiJS[I][J];
                    }
                  }
                }
              }
            }
            for (J=0; J<Factorizer.NumberLength; J++) {
              aiJ00[0][J] = aiJ01[0][J] = MontgomeryParams.MontgomeryMultR1[J];
              for (I=1; I<PK; I++) {
                aiJ00[I][J] = aiJ01[I][J] = 0;
              }
            }
            VK = (int)BigNumOps.BigNbrModLong(Factorizer.TestNbr,PK);
            for (I=1; I<PK; I++) {
              if (I%P != 0) {
                U1 = 1; U3 = I; V1 = 0; V3 = PK;
                while (V3 != 0) {
                  QQ = U3/V3;
                  T1 = U1 - V1*QQ;
                  T3 = U3 - V3*QQ;
                  U1 = V1; U3 = V3;
                  V1 = T1; V3 = T3;
                }
                aiInv[I] = (U1+PK)%PK;
              }
              else {
                aiInv[I] = 0;
              }
            }
            if (P != 2) {
              for (IV = 0; IV <= 1; IV++) {
                for (X = 1; X < PK; X++) {
                  for (I=0; I<PK; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJS[I][J] = aiJ0[I][J];
                    }
                  }
                  if (X%P == 0) {continue;}
                  if (IV == 0) {
                    BigNumOps.LongToBigNbr(X, biExp);
                  }
                  else {
                    BigNumOps.LongToBigNbr(VK * X / PK, biExp);
                    if (VK * X / PK == 0) {continue;}
                  }
                  BigNumOps.JS_E(PK, PL, PM, P);
                  for (I=0; I<PK; I++) {
                    for (J=0; J<Factorizer.NumberLength; J++) {
                      aiJW[I][J] = 0;
                    }
                  }
                  InvX = aiInv[X];
                  for (I=0; I<PK; I++) {
                    J = I*InvX % PK;
                    BigNumOps.AddBigNbrModN(aiJW[J],aiJS[I],aiJW[J]);
                  }
                  BigNumOps.NormalizeJW(PK, PL, PM, P);
                  if (IV == 0) {
                    for (I=0; I<PK; I++) {
                      for (J=0; J<Factorizer.NumberLength; J++) {
                        aiJS[I][J] = aiJ00[I][J];
                      }
                    }
                  }
                  else {
                    for (I=0; I<PK; I++) {
                      for (J=0; J<Factorizer.NumberLength; J++) {
                        aiJS[I][J] = aiJ01[I][J];
                      }
                    }
                  }
                  BigNumOps.JS_JW(PK, PL, PM, P);
                  if (IV == 0) {
                    for (I=0; I<PK; I++) {
                      for (J=0; J<Factorizer.NumberLength; J++) {
                        aiJ00[I][J] = aiJS[I][J];
                      }
                    }
                  }
                  else {
                    for (I=0; I<PK; I++) {
                      for (J=0; J<Factorizer.NumberLength; J++) {
                        aiJ01[I][J] = aiJS[I][J];
                      }
                    }
                  }
                }                // end for X
              }                  // end for IV
            }