montgomeryparams.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
 */

/**
 * 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 MontgomeryParams
{
      static final long MontgomeryMultR1[] = new long[Constants.MP_INITIAL_SIZE];
      static final long MontgomeryMultR2[] = new long[Constants.MP_INITIAL_SIZE];
      static  final long montyAfterInv[] = new long[Constants.MP_INITIAL_SIZE];
      static long montyMultN;

      /**
       *
       */
      static void GetMontgomeryParams() {
            int _length = Factorizer.NumberLength;
            long[] n= Factorizer.TestNbr;

            int N, x, j;

            Factorizer.dN = (double) n[_length - 1]; //TestNbr[_length - 1];
            if (_length > 1) {
                  Factorizer.dN += (double) n[_length - 2] / Constants.dDosALa32; // MAX_WORD;
            }
            if (_length > 2) {
                  Factorizer.dN += (double) n[_length - 3] / Constants.dDosALa64; // DTWO_POWER64;
            }

            x = N = (int) n[0]; // 2 least significant bits of inverse correct.
            x = x * (2 - N * x); // 4 least significant bits of inverse correct.
            x = x * (2 - N * x); // 8 least significant bits of inverse correct.
            x = x * (2 - N * x); // 16 least significant bits of inverse correct.
            x = x * (2 - N * x); // 32 least significant bits of inverse correct.
            montyMultN = ( -x) & Constants.MaxUInt;
            j = _length;
            MontgomeryMultR1[j] = 1;
            do {
                  MontgomeryMultR1[--j] = 0;
            }
            while (j > 0);

//            _ops = new BigNumOps(_mpx); //, dN);
            BigNumOps.AdjustModN(MontgomeryMultR1);
            BigNumOps.MultBigNbrModN(MontgomeryMultR1, MontgomeryMultR1, MontgomeryMultR2);
            MontgomeryMult(MontgomeryMultR2, MontgomeryMultR2, montyAfterInv);
            BigNumOps.AddBigNbrModN(MontgomeryMultR1, MontgomeryMultR1, MontgomeryMultR2);
      }


      /**
       *
       */
      static void MontgomeryMult( long Nbr1[], long Nbr2[], long Prod[]) {
	    int i, j;
	    long MaxUInt = Constants.MaxUInt;
	    long Pr, Nbr, ProdB0;
	    long Prod0, Prod1, Prod2, Prod3, Prod4, Prod5, Prod6, Prod7;
	    long Prod8, Prod9, Prod10;
	    long TestNbr0, TestNbr1, TestNbr2, TestNbr3, TestNbr4, TestNbr5,
		  TestNbr6, TestNbr7;
	    long TestNbr8, TestNbr9, TestNbr10;
	    long Nbr2_0, Nbr2_1, Nbr2_2, Nbr2_3, Nbr2_4, Nbr2_5, Nbr2_6, Nbr2_7;
	    long Nbr2_8, Nbr2_9, Nbr2_10;
	    long New, MontDig, Prd;

	    long TestNbr[] = Factorizer.TestNbr;
	    int _length = Factorizer.NumberLength;
	    int NumberLength1 = _length - 1;
//	    String EcmInfo;


	    TestNbr0 = TestNbr[0];
	    TestNbr1 = TestNbr[1];
	    Nbr2_0 = Nbr2[0];
	    Nbr2_1 = Nbr2[1];
	    switch (_length) {
		  case 2:
			Prod0 = Prod1 = 0;
			i = 0;
			do {
			      ProdB0 = (Nbr = Nbr1[i]) * Nbr2_0 + Prod0;
			      MontDig = ( (int) ProdB0 * montyMultN) &
				    MaxUInt;
			      Prd = (ProdB0 >>> 32) + Nbr * Nbr2_1 + Prod1;
			      Prod0 = (Pr =
				       ( (MontDig * TestNbr0 + MaxUInt) >>> 32) +
				       MontDig * TestNbr1 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prod1 = (Prd >>> 32) + (Pr >>> 32);
			      i++;
			}
			while (i < 2);
			if (Prod1 > TestNbr1 || Prod1 == TestNbr1 &&
			    (Prod0 >= TestNbr0)) {
			      Prod0 = (Pr = Prod0 - TestNbr0) & MaxUInt;
			      Prod1 = ( (Pr >> 32) + Prod1 - TestNbr1) &
				    MaxUInt;
			}
			Prod[0] = Prod0;
			Prod[1] = Prod1;
			break;

		  case 3:
			Prod0 = Prod1 = Prod2 = 0;
			TestNbr2 = TestNbr[2];
			Nbr2_2 = Nbr2[2];
			i = 0;
			do {
			      ProdB0 = (Nbr = Nbr1[i]) * Nbr2_0 + Prod0;
			      MontDig = ( (int) ProdB0 * montyMultN) &
				    MaxUInt;
			      Prd = (ProdB0 >>> 32) + Nbr * Nbr2_1 + Prod1;
			      Prod0 = (Pr =
				       ( (MontDig * TestNbr0 + MaxUInt) >>> 32) +
				       MontDig * TestNbr1 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_2 + Prod2;
			      Prod1 = (Pr = (Pr >>> 32) + MontDig * TestNbr2 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prod2 = (Prd >>> 32) + (Pr >>> 32);
			      i++;
			}
			while (i < 3);
			if (Prod2 > TestNbr2 || Prod2 == TestNbr2 &&
			    (Prod1 > TestNbr1 || Prod1 == TestNbr1 &&
			     (Prod0 >= TestNbr0))) {
			      Prod0 = (Pr = Prod0 - TestNbr0) & MaxUInt;
			      Prod1 = (Pr = (Pr >> 32) + Prod1 - TestNbr1) &
				    MaxUInt;
			      Prod2 = ( (Pr >> 32) + Prod2 - TestNbr2) &
				    MaxUInt;
			}
			Prod[0] = Prod0;
			Prod[1] = Prod1;
			Prod[2] = Prod2;
			break;

		  case 4:
			Prod0 = Prod1 = Prod2 = Prod3 = 0;
			TestNbr2 = TestNbr[2];
			TestNbr3 = TestNbr[3];
			Nbr2_2 = Nbr2[2];
			Nbr2_3 = Nbr2[3];
			i = 0;
			do {
			      ProdB0 = (Nbr = Nbr1[i]) * Nbr2_0 + Prod0;
			      MontDig = ( (int) ProdB0 * montyMultN) &
				    MaxUInt;
			      Prd = (ProdB0 >>> 32) + Nbr * Nbr2_1 + Prod1;
			      Prod0 = (Pr =
				       ( (MontDig * TestNbr0 + MaxUInt) >>> 32) +
				       MontDig * TestNbr1 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_2 + Prod2;
			      Prod1 = (Pr = (Pr >>> 32) + MontDig * TestNbr2 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_3 + Prod3;
			      Prod2 = (Pr = (Pr >>> 32) + MontDig * TestNbr3 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prod3 = (Prd >>> 32) + (Pr >>> 32);
			      i++;
			}
			while (i < 4);
			if (Prod3 > TestNbr3 || Prod3 == TestNbr3 &&
			    (Prod2 > TestNbr2 || Prod2 == TestNbr2 &&
			     (Prod1 > TestNbr1 || Prod1 == TestNbr1 &&
			      (Prod0 >= TestNbr0)))) {
			      Prod0 = (Pr = Prod0 - TestNbr0) & MaxUInt;
			      Prod1 = (Pr = (Pr >> 32) + Prod1 - TestNbr1) &
				    MaxUInt;
			      Prod2 = (Pr = (Pr >> 32) + Prod2 - TestNbr2) &
				    MaxUInt;
			      Prod3 = ( (Pr >> 32) + Prod3 - TestNbr3) &
				    MaxUInt;
			}
			Prod[0] = Prod0;
			Prod[1] = Prod1;
			Prod[2] = Prod2;
			Prod[3] = Prod3;
			break;

		  case 5:
			Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = 0;
			TestNbr2 = TestNbr[2];
			TestNbr3 = TestNbr[3];
			TestNbr4 = TestNbr[4];
			Nbr2_2 = Nbr2[2];
			Nbr2_3 = Nbr2[3];
			Nbr2_4 = Nbr2[4];
			i = 0;
			do {
			      ProdB0 = (Nbr = Nbr1[i]) * Nbr2_0 + Prod0;
			      MontDig = ( (int) ProdB0 * montyMultN) &
				    MaxUInt;
			      Prd = (ProdB0 >>> 32) + Nbr * Nbr2_1 + Prod1;
			      Prod0 = (Pr =
				       ( (MontDig * TestNbr0 + MaxUInt) >>> 32) +
				       MontDig * TestNbr1 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_2 + Prod2;
			      Prod1 = (Pr = (Pr >>> 32) + MontDig * TestNbr2 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_3 + Prod3;
			      Prod2 = (Pr = (Pr >>> 32) + MontDig * TestNbr3 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_4 + Prod4;
			      Prod3 = (Pr = (Pr >>> 32) + MontDig * TestNbr4 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prod4 = (Prd >>> 32) + (Pr >>> 32);
			      i++;
			}
			while (i < 5);
			if (Prod4 > TestNbr4 || Prod4 == TestNbr4 &&
			    (Prod3 > TestNbr3 || Prod3 == TestNbr3 &&
			     (Prod2 > TestNbr2 || Prod2 == TestNbr2 &&
			      (Prod1 > TestNbr1 || Prod1 == TestNbr1 &&
			       (Prod0 >= TestNbr0))))) {
			      Prod0 = (Pr = Prod0 - TestNbr0) & MaxUInt;
			      Prod1 = (Pr = (Pr >> 32) + Prod1 - TestNbr1) &
				    MaxUInt;
			      Prod2 = (Pr = (Pr >> 32) + Prod2 - TestNbr2) &
				    MaxUInt;
			      Prod3 = (Pr = (Pr >> 32) + Prod3 - TestNbr3) &
				    MaxUInt;
			      Prod4 = ( (Pr >> 32) + Prod4 - TestNbr4) &
				    MaxUInt;
			}
			Prod[0] = Prod0;
			Prod[1] = Prod1;
			Prod[2] = Prod2;
			Prod[3] = Prod3;
			Prod[4] = Prod4;
			break;

		  case 6:
			Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = 0;
			TestNbr2 = TestNbr[2];
			TestNbr3 = TestNbr[3];
			TestNbr4 = TestNbr[4];
			TestNbr5 = TestNbr[5];
			Nbr2_2 = Nbr2[2];
			Nbr2_3 = Nbr2[3];
			Nbr2_4 = Nbr2[4];
			Nbr2_5 = Nbr2[5];
			i = 0;
			do {
			      ProdB0 = (Nbr = Nbr1[i]) * Nbr2_0 + Prod0;
			      MontDig = ( (int) ProdB0 * montyMultN) &
				    MaxUInt;
			      Prd = (ProdB0 >>> 32) + Nbr * Nbr2_1 + Prod1;
			      Prod0 = (Pr =
				       ( (MontDig * TestNbr0 + MaxUInt) >>> 32) +
				       MontDig * TestNbr1 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_2 + Prod2;
			      Prod1 = (Pr = (Pr >>> 32) + MontDig * TestNbr2 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_3 + Prod3;
			      Prod2 = (Pr = (Pr >>> 32) + MontDig * TestNbr3 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_4 + Prod4;
			      Prod3 = (Pr = (Pr >>> 32) + MontDig * TestNbr4 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prd = (Prd >>> 32) + Nbr * Nbr2_5 + Prod5;
			      Prod4 = (Pr = (Pr >>> 32) + MontDig * TestNbr5 +
				       (Prd & MaxUInt)) & MaxUInt;
			      Prod5 = (Prd >>> 32) + (Pr >>> 32);
			      i++;
			}
			while (i < 6);
			if (Prod5 > TestNbr5 || Prod5 == TestNbr5 &&
			    (Prod4 > TestNbr4 || Prod4 == TestNbr4 &&
			     (Prod3 > TestNbr3 || Prod3 == TestNbr3 &&
			      (Prod2 > TestNbr2 || Prod2 == TestNbr2 &&
			       (Prod1 > TestNbr1 || Prod1 == TestNbr1 &&
				(Prod0 >= TestNbr0)))))) {
			      Prod0 = (Pr = Prod0 - TestNbr0) & MaxUInt;
			      Prod1 = (Pr = (Pr >> 32) + Prod1 - TestNbr1) &
				    MaxUInt;
			      Prod2 = (Pr = (Pr >> 32) + Prod2 - TestNbr2) &
				    MaxUInt;
			      Prod3 = (Pr = (Pr >> 32) + Prod3 - TestNbr3) &