// <XMP>
// Elliptic Curve Method (ECM) Prime Factorization
//
// Written by Dario Alejandro Alpern (Buenos Aires - Argentina)
// Last updated October 16th, 2003. See http://www.alpertron.com.ar/ECM.HTM
//
// Based in Yuji Kida's implementation for UBASIC interpreter
//
// No part of this code can be used for commercial purposes without
// the written consent from the author. Otherwise it can be used freely
// except that you have to write somewhere in the code this header.
// 

import java.math.BigInteger;
import java.util.Random;
import java.util.Date;

public class ecm31 {
	static BigInteger zero = new BigInteger("0");
	static BigInteger one = new BigInteger("1");
	static BigInteger two = new BigInteger("2"); 

	static BigInteger n = new BigInteger("0");
	static BigInteger startCurve = new BigInteger("0");
	static String s;

	static int olds1 = 0;
	static int s1 = 0;


  static final int NLen = 1200;


  int NumberLength;          /* Length of multiple precision nbrs */

  final long CalcAuxGcdU[] = new long[NLen];/* Used inside GCD calculations in */
  final long CalcAuxGcdV[] = new long[NLen];/* multiple precision numbers */
  final long CalcAuxGcdT[] = new long[NLen];


  long TestNbr[] = new long[NLen];


  final long GcdAccumulated[] = new long[NLen];



  long [] fieldAA, fieldTX, fieldTZ, fieldUX, fieldUZ;
  long [] fieldAux1, fieldAux2, fieldAux3, fieldAux4;



  static final long DosALa32 = (long)1 << 32;
  static final long DosALa31 = (long)1 << 31;

  static final double dDosALa32 = (double)DosALa32;
  static final double dDosALa64 = dDosALa32 * dDosALa32;

  static final long MaxUInt = DosALa32 - 1;
  static final long Mi = 1000000000;
  double dN;
  final long BigNbr1[] = new long[NLen];
  final int SmallPrime[] = new int[670];    /* Primes < 5000 */
  final long MontgomeryMultR1[] = new long[NLen];
  final long MontgomeryMultR2[] = new long[NLen];
  final long MontgomeryMultAfterInv[] = new long[NLen];
  long MontgomeryMultN;


  int indexM,maxIndexM;


	public boolean factorize(BigInteger n, BigInteger startCurve, BigInteger finishCurve) {
		BigInteger factor = new BigInteger("1");
		boolean prime = false;

		if (n.isProbablePrime(1000)) {
			prime = true;
//			System.out.println('\n' + n.toString());
			return prime;
		}


// bugfix JGW 2006-06-30
		this.startCurve = startCurve;


		factor = fnECM(n, this.startCurve, finishCurve);



		if (factor.compareTo(one) > 0 && factor.compareTo(n) < 0) {
//			factorize(factor, startCurve, finishCurve);
			factorize(n.divide(factor), this.startCurve, finishCurve);
		}
		else {
			prime = false;
//			System.out.println('\n' + "composite" + '\t' + n);
		}

		return prime;

	}

	public BigInteger evalRand(char op, BigInteger oldn) {
		BigInteger n = new BigInteger("1");


		switch (op) {
		case 'r':
		case 'R':
			Random r = new Random();

			n = new BigInteger(oldn.intValue(), r);
			break;

		default:
			n = oldn;
			break;
		}
		return n;
	}

	public BigInteger evalFact(BigInteger oldn, char op) {
		BigInteger n = new BigInteger("1");
		BigInteger i = new BigInteger("1");
		BigInteger j = new BigInteger("1");
		boolean prime = true;


		switch (op) {
		case '!':
			for (i = one; i.compareTo(oldn) <= 0; i = i.add(one))
				n = n.multiply(i);
			break;

		case '#':
			for (i = one; i.compareTo(oldn) <= 0; i = i.add(one)) {
				prime = true;
				for (j = two; (prime == true) && (j.multiply(j).compareTo(i) <= 0); j = j.add(one))
					if (i.remainder(j).compareTo(zero) == 0)
						prime = false;
				if (prime == true)
					n = n.multiply(i);
			}
			break;

		default:
			n = oldn;
			break;
		}
		return n;
	}

	public BigInteger evalPower(BigInteger oldn, BigInteger n1, char op) {
		BigInteger n = new BigInteger("0");

		switch (op) {
		case '^':
			n = oldn.pow(n1.intValue());
			break;
		default:
			n = n1;
			break;
		}
		return n;
	}

	public BigInteger evalProduct(BigInteger oldn, BigInteger n1, char op) {
		BigInteger n = new BigInteger("0");

		switch (op) {
		case '*':
			n = oldn.multiply(n1);
			break;
		case '/':
			n = oldn.divide(n1);
			break;
		case '%':
			n = oldn.remainder(n1);
			break;
		default:
			n = n1;
			break;
		}
		return n;
	}

	public BigInteger evalSum(BigInteger oldn, BigInteger n1, char op) {
		BigInteger n = new BigInteger("0");

		switch (op) {
		case '+':
			n = oldn.add(n1);
			break;
		case '-':
			n = oldn.subtract(n1);
			break;
		default:
			n = n1;
			break;
		}
		return n;
	}

	public BigInteger eval(String s) {
		BigInteger oldn0 = new BigInteger("0");
		BigInteger oldn1 = new BigInteger("0");
		BigInteger oldn2 = new BigInteger("0");
		BigInteger n = new BigInteger("0");

		char oldop0 = 0;
		char oldop1 = 0;
		char oldop2 = 0;
		char op = 0;

		while (s1 < s.length()) {
			switch (s.charAt(s1)) {
				case '(':
				case '[':
				case '{':
					s1++;
					n = eval(s);
					break;

				case '0':
				case '1':
				case '2':
				case '3':
				case '4':
				case '5':
				case '6':
				case '7':
				case '8':
				case '9':
					n = readNum(s);
					break;

				default:
					break;
			}
			if (s1 < s.length()) {
				switch (s.charAt(s1)) {

				case ')':
				case ']':
				case '}':
				case '!':
				case '#':
				case 'r':
				case 'R':
				case '^':
				case '*':
				case '/':
				case '%':
				case '+':
				case '-':
					op = s.charAt(s1);
					s1++;
					break;

				default:
					break;
				}
			}
			else
				op = 0;


			switch (op) {
				case 0:
				case ')':
				case ']':
				case '}':
					n = evalPower(oldn2, n, oldop2);
					n = evalProduct(oldn1, n, oldop1);
					n = evalSum(oldn0, n, oldop0);
					return n;

				case '!':
				case '#':
					n = evalFact(n, op);
					break;

				case 'r':
				case 'R':
					n = readNum(s);
					n = evalRand(op, n);
					break;

				case '^':
					n = evalPower(oldn2, n, oldop2);
					oldn2 = n;
					oldop2 = op;
					break;

				case '*':
				case '/':
				case '%':
					n = evalPower(oldn2, n, oldop2);
					oldop2 = 0;
					n = evalProduct(oldn1, n, oldop1);
					oldn1 = n;
					oldop1 = op;
					break;

				case '+':
				case '-':
					n = evalPower(oldn2, n, oldop2);
					oldop2 = 0;
					n = evalProduct(oldn1, n, oldop1);
					oldop1 = 0;
					n = evalSum(oldn0, n, oldop0);
					oldn0 = n;
					oldop0 = op;
					break;
				default:
					break;
			}
		}
		return n;
	}

	public BigInteger readNum(String s) {
		olds1 = s1;
		while (s1 < s.length() && Character.isDigit(s.charAt(s1)))
			s1++;
		BigInteger n = new BigInteger(s.substring(olds1, s1));

		return n;
	}





  BigInteger fnECM(BigInteger N, BigInteger startCurve, BigInteger finishCurve) {

// added JGW
  int EC;                    /* Elliptic Curve number */


    int I,J,Pass,Qaux;
    long L1,L2,LS,P,Q,K,IP,Paux=1;
    long [] A0 = new long[NLen];
    long [] A02 = new long[NLen];
    long [] A03 = new long[NLen];
    long [] AA = new long[NLen];
    long [] DX = new long[NLen];
    long [] DZ = new long[NLen];
    long [] GD = new long[NLen];
    long [] M = new long[NLen];
    long [] N1 = new long[NLen];
    long [] N2 = new long[NLen];
    long [] TX = new long[NLen]; fieldTX = TX;
    long [] TZ = new long[NLen]; fieldTZ = TZ;
    long [] UX = new long[NLen]; fieldUX = UX;
    long [] UZ = new long[NLen]; fieldUZ = UZ;
    long [] W0 = new long[NLen];
    long [] W1 = new long[NLen];
    long [] W2 = new long[NLen];
    long [] W3 = new long[NLen];
    long [] W4 = new long[NLen];
    long [] WX = new long[NLen];
    long [] WZ = new long[NLen];
    long [] X = new long[NLen];
    long [] Z = new long[NLen];
    long [] Aux1 = new long[NLen]; fieldAux1 = Aux1;
    long [] Aux2 = new long[NLen]; fieldAux2 = Aux2;
    long [] Aux3 = new long[NLen]; fieldAux3 = Aux3;
    long [] Aux4 = new long[NLen]; fieldAux4 = Aux4;
    long [] Xaux = new long[NLen];
    long [] Zaux = new long[NLen];
    long [][] root = new long[480][NLen];
    byte [] Result;
    byte [] sieve = new byte[23100];
    byte [] sieve2310 = new byte[2310];
    int [] sieveidx = new int[480];
    String UpperLine;
    String LowerLine;
    String primalityString;
    int Prob, i, j, u;

    fieldAA = AA;
    BigNbrToBigInt(N);
    GetMontgomeryParms();
    for (I=0; I<NumberLength; I++) {
      M[I] = DX[I] = DZ[I] = W3[I] = W4[I] = GD[I] = 0;
    }

// added JGW
	EC = startCurve.intValue()-1;
	if (EC<0) EC=0;
	

// added JGW
    BigNbr1[0] = 1;
    for (i=1; i<NLen; i++) {
      BigNbr1[i] = 0;
    }



    SmallPrime[0] = 2;
    P = 3;
    indexM = 1;
    for (indexM = 1; indexM < SmallPrime.length; indexM++) {
      SmallPrime[indexM] = (int)P;    /* Store prime */
      calculate_new_prime1:
        do {
        P+=2;
        for (Q=3; Q*Q<=P; Q+=2) {     /* Check if P is prime */
          if (P%Q == 0) {
            continue calculate_new_prime1;   /* Composite */
          }
        }
        break;                      /* Prime found */
      } while (true);
    }

    do {
      new_curve:
        do {
          EC++;

// added JGW
	System.out.print("curve " + EC + '\r');
	
	
// bugfix JGW 2006-06-30
	this.startCurve = this.startCurve.add(one);
	
	
// added JGW 06-06-27
	if (EC > finishCurve.intValue()) return N;
	

        L1 = 2000;
        L2 = 200000;
        LS = 45;
        Paux = EC;
        if (EC > 25) {
          if (EC < 326) {
            L1 = 50000;
            L2 = 5000000;
            LS = 224;
            Paux = EC - 24;
          }
          else {
            if (EC < 2000) {
              L1 = 1000000;
              L2 = 100000000;
              LS = 1001;
              Paux = EC - 299;
            }
            else {
              L1 = 11000000;
              L2 = 1100000000;
              LS = 3316;
              Paux = EC - 1900;
            }
          }
        }

        LongToBigNbr(2*(EC+1), Aux1);
        LongToBigNbr(3*(EC+1)*(EC+1)-1, Aux2);
        ModInvBigNbr(Aux2, Aux2, TestNbr);
        MultBigNbrModN(Aux1, Aux2, Aux3);
        MultBigNbrModN(Aux3,MontgomeryMultR1,A0);
        MontgomeryMult(A0, A0, A02);
        MontgomeryMult(A02, A0, A03);
        SubtractBigNbrModN(A03, A0, Aux1);
        MultBigNbrByLongModN(A02, 9, Aux2);
        SubtractBigNbrModN(Aux2, MontgomeryMultR1, Aux2);
        MontgomeryMult(Aux1, Aux2, Aux3);
        if (BigNbrIsZero(Aux3)) {
          continue;
        }
        MultBigNbrByLongModN(A0, 4, Z);
        MultBigNbrByLongModN(A02, 6, Aux1);
        SubtractBigNbrModN(MontgomeryMultR1, Aux1, Aux1);
        MontgomeryMult(A02, A02, Aux2);
        MultBigNbrByLongModN(Aux2, 3, Aux2);
        SubtractBigNbrModN(Aux1, Aux2, Aux1);
        MultBigNbrByLongModN(A03, 4, Aux2);
        ModInvBigNbr(Aux2, Aux2, TestNbr);
        MontgomeryMult(Aux2, MontgomeryMultAfterInv, Aux3);
        MontgomeryMult(Aux1, Aux3, A0);
        AddBigNbrModN(A0, MontgomeryMultR2, Aux1);
        LongToBigNbr(4, Aux2);
        ModInvBigNbr(Aux2, Aux3, TestNbr);
        MultBigNbrModN(Aux3,MontgomeryMultR1,Aux2);
        MontgomeryMult(Aux1, Aux2, AA);
        MultBigNbrByLongModN(A02, 3, Aux1);
        AddBigNbrModN(Aux1, MontgomeryMultR1, X);
        /**************/
        /* First step */
        /**************/
        System.arraycopy(X,0,Xaux,0,NumberLength);
        System.arraycopy(Z,0,Zaux,0,NumberLength);
        System.arraycopy(MontgomeryMultR1,0,GcdAccumulated,0,NumberLength);
        for (Pass=0; Pass<2; Pass++) {
          /* For powers of 2 */

          for (I=1; I<=L1; I<<=1) {
            duplicate(X, Z, X, Z);
          }
          for (I=3; I<=L1; I*=3) {
            duplicate(W1, W2, X, Z);
            add3(X, Z, X, Z, W1, W2, X, Z);
          }

          if (Pass == 0) {
            MontgomeryMult(GcdAccumulated, Z, Aux1);
            System.arraycopy(Aux1,0,GcdAccumulated,0,NumberLength);
          }
          else {
            GcdBigNbr(Z, TestNbr, GD);
            if (BigNbrAreEqual(GD, BigNbr1) == false) {
              break new_curve;
            }
          }

          /* for powers of odd primes */

          indexM=1;
          do {

            P = SmallPrime[indexM];
            for (IP=P; IP<=L1; IP*=P) {
              prac((int)P, X, Z, W1, W2, W3, W4);
            }
            indexM++;
            if (Pass == 0) {
              MontgomeryMult(GcdAccumulated, Z, Aux1);
              System.arraycopy(Aux1,0,GcdAccumulated,0,NumberLength);
            }
            else {
              GcdBigNbr(Z, TestNbr, GD);
              if (BigNbrAreEqual(GD, BigNbr1) == false) {
                break new_curve;
              }
            }
          } while (SmallPrime[indexM-1]<=LS);
          P+=2;

         /* Initialize sieve2310[n]: 1 if gcd(P+2n,2310) > 1, 0 otherwise */
          u = (int)P;
          for (i=0; i<2310; i++) {
            sieve2310[i] = (u%3 == 0 || u%5 == 0 || u%7 == 0 || u%11 == 0?
                                         (byte)1:(byte)0);
            u += 2;
          }
          do {
             /* Generate sieve */
            GenerateSieve((int)P, sieve, sieve2310, SmallPrime);

            /* Walk through sieve */

            for (i=0; i<23100; i++) {
              if (sieve[i] != 0) continue;   /* Do not process composites */
              if (P+2*i > L1) break;

              prac((int)(P+2*i), X, Z, W1, W2, W3, W4);
              if (Pass == 0) {
                MontgomeryMult(GcdAccumulated, Z, Aux1);
                System.arraycopy(Aux1,0,GcdAccumulated,0,NumberLength);
              }
              else {
                GcdBigNbr(Z, TestNbr, GD);
                if (BigNbrAreEqual(GD, BigNbr1) == false) {
                  break new_curve;
                }
              }
            }
            P += 46200;
          } while (P<L1);
          if (Pass == 0) {
            if (BigNbrIsZero(GcdAccumulated)) { // If GcdAccumulated is
              System.arraycopy(Xaux,0,X,0,NumberLength);
              System.arraycopy(Zaux,0,Z,0,NumberLength);
              continue;                        // multiple of TestNbr, continue.
            }
            GcdBigNbr(GcdAccumulated, TestNbr, GD);
            if (BigNbrAreEqual(GD, BigNbr1) == false) {
              break new_curve;
            }
            break;
          }
        }                     /* end for Pass */

        /******************************************************/
        /* Second step (using improved standard continuation) */
        /******************************************************/

        j = 0;
        for (u=1; u<2310; u+=2) {
          if (u%3 == 0 || u%5 == 0 || u%7 == 0 || u%11 == 0) {
            sieve2310[u/2] = (byte)1;
            }
          else {
            sieve2310[(sieveidx[j++] = u/2)] = (byte)0;
            }
          }
        System.arraycopy(sieve2310,0,sieve2310,1155,1155);
        System.arraycopy(X,0,Xaux,0,NumberLength);    // (X:Z) -> Q (output
        System.arraycopy(Z,0,Zaux,0,NumberLength);    //         from step 1)
        for (Pass=0; Pass<2; Pass++) {
          System.arraycopy(MontgomeryMultR1,0,GcdAccumulated,0,NumberLength);
          System.arraycopy(X,0,UX,0,NumberLength);
          System.arraycopy(Z,0,UZ,0,NumberLength);    // (UX:UZ) -> Q
          ModInvBigNbr(Z, Aux2, TestNbr);
          MontgomeryMult(Aux2, MontgomeryMultAfterInv, Aux1);
          MontgomeryMult(Aux1, X, root[0]);           // root[0] <- X/Z (Q)
          J = 0;
          AddBigNbrModN(X, Z, Aux1);
          MontgomeryMult(Aux1, Aux1, W1);
          SubtractBigNbrModN(X, Z, Aux1);
          MontgomeryMult(Aux1, Aux1, W2);
          MontgomeryMult(W1, W2, TX);
          SubtractBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, AA, Aux2);
          AddBigNbrModN(Aux2, W2, Aux3);
          MontgomeryMult(Aux1, Aux3, TZ);             // (TX:TZ) -> 2Q
          SubtractBigNbrModN(X, Z, Aux1);
          AddBigNbrModN(TX, TZ, Aux2);                
          MontgomeryMult(Aux1, Aux2, W1);
          AddBigNbrModN(X, Z, Aux1);
          SubtractBigNbrModN(TX, TZ, Aux2);
          MontgomeryMult(Aux1, Aux2, W2);
          AddBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, Aux1, Aux2);
          MontgomeryMult(Aux2, UZ, X);
          SubtractBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, Aux1, Aux2);
          MontgomeryMult(Aux2, UX, Z);                // (X:Z) -> 3Q
          for (I=5; I<2310; I+=2) {
            System.arraycopy(X,0,WX,0,NumberLength);
            System.arraycopy(Z,0,WZ,0,NumberLength);
            SubtractBigNbrModN(X, Z, Aux1);
            AddBigNbrModN(TX, TZ, Aux2);
            MontgomeryMult(Aux1, Aux2, W1);
            AddBigNbrModN(X, Z, Aux1);
            SubtractBigNbrModN(TX, TZ, Aux2);
            MontgomeryMult(Aux1, Aux2, W2);
            AddBigNbrModN(W1, W2, Aux1);
            MontgomeryMult(Aux1, Aux1, Aux2);
            MontgomeryMult(Aux2, UZ, X);
            SubtractBigNbrModN(W1, W2, Aux1);
            MontgomeryMult(Aux1, Aux1, Aux2);
            MontgomeryMult(Aux2, UX, Z);           // (X:Z) -> 5Q, 7Q, ...
            if (Pass == 0) {
              MontgomeryMult(GcdAccumulated, Aux1, Aux2);
              System.arraycopy(Aux2,0,GcdAccumulated,0,NumberLength);
            }
            else {
              GcdBigNbr(Aux1, TestNbr, GD);
              if (BigNbrAreEqual(GD, BigNbr1) == false) {
                break new_curve;
              }
            }
            if (I == 1155) {
              System.arraycopy(X,0,DX,0,NumberLength);
              System.arraycopy(Z,0,DZ,0,NumberLength);  // (DX:DZ) -> 1155Q
            }
            if (I%3 != 0 && I%5 != 0 && I%7 != 0 && I%11 != 0) {
              J++;
              ModInvBigNbr(Z, Aux2, TestNbr);
              MontgomeryMult(Aux2, MontgomeryMultAfterInv, Aux1);
              MontgomeryMult(Aux1, X, root[J]);         // root[J] <- X/Z
            }
            System.arraycopy(WX,0,UX,0,NumberLength);   // (UX:UZ) <-
            System.arraycopy(WZ,0,UZ,0,NumberLength);   // Previous (X:Z)
          }                         /* end for I */
          AddBigNbrModN(DX, DZ, Aux1);
          MontgomeryMult(Aux1, Aux1, W1);
          SubtractBigNbrModN(DX, DZ, Aux1);
          MontgomeryMult(Aux1, Aux1, W2);
          MontgomeryMult(W1, W2, X);
          SubtractBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, AA, Aux2);
          AddBigNbrModN(Aux2, W2, Aux3);
          MontgomeryMult(Aux1, Aux3, Z);
          System.arraycopy(X,0,UX,0,NumberLength);
          System.arraycopy(Z,0,UZ,0,NumberLength);    // (UX:UZ) -> 2310Q
          AddBigNbrModN(X, Z, Aux1);
          MontgomeryMult(Aux1, Aux1, W1);
          SubtractBigNbrModN(X, Z, Aux1);
          MontgomeryMult(Aux1, Aux1, W2);
          MontgomeryMult(W1, W2, TX);
          SubtractBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, AA, Aux2);
          AddBigNbrModN(Aux2, W2, Aux3);
          MontgomeryMult(Aux1, Aux3, TZ);             // (TX:TZ) -> 2*2310Q
          SubtractBigNbrModN(X, Z, Aux1);
          AddBigNbrModN(TX, TZ, Aux2);
          MontgomeryMult(Aux1, Aux2, W1);
          AddBigNbrModN(X, Z, Aux1);
          SubtractBigNbrModN(TX, TZ, Aux2);
          MontgomeryMult(Aux1, Aux2, W2);
          AddBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, Aux1, Aux2);
          MontgomeryMult(Aux2, UZ, X);
          SubtractBigNbrModN(W1, W2, Aux1);
          MontgomeryMult(Aux1, Aux1, Aux2);
          MontgomeryMult(Aux2, UX, Z);                // (X:Z) -> 3*2310Q
          Qaux = (int)(L1/4620);
          maxIndexM = (int)(L2/4620);
          for (indexM=0; indexM<=maxIndexM; indexM++) {
            if (indexM >= Qaux) {            // If inside step 2 range... 
              if (indexM == 0) {
                ModInvBigNbr(UZ, Aux2, TestNbr);
                MontgomeryMult(Aux2, MontgomeryMultAfterInv, Aux3);
                MontgomeryMult(UX, Aux3, Aux1);         // Aux1 <- X/Z (2310Q)
                }
              else {
                ModInvBigNbr(Z, Aux2, TestNbr);
                MontgomeryMult(Aux2, MontgomeryMultAfterInv, Aux3);
                MontgomeryMult(X, Aux3, Aux1);          // Aux1 <- X/Z (3,5,*
                }                                       //              2310Q)

               /* Generate sieve */
              if (indexM % 10 == 0 || indexM == Qaux) {
                GenerateSieve(indexM/10*46200+1, sieve, sieve2310, SmallPrime);
              }
              /* Walk through sieve */
              J = 1155 + (indexM%10) * 2310;
              for (i=0; i<480; i++) {
                j = sieveidx[i];                        // 0 < J < 1155
                if (sieve[J+j] != 0 && sieve[J-1-j] != 0) {
                  continue;   // Do not process if both are composite numbers.
                  }
                SubtractBigNbrModN(Aux1, root[i], M);
                MontgomeryMult(GcdAccumulated, M, Aux2);
                System.arraycopy(Aux2,0,GcdAccumulated,0,NumberLength);
              }
              if (Pass != 0) {
                GcdBigNbr(GcdAccumulated, TestNbr, GD);
                if (BigNbrAreEqual(GD, BigNbr1) == false) {
                  break new_curve;
                }
              }
            }
            if (indexM != 0) {      // Update (X:Z)
              System.arraycopy(X,0,WX,0,NumberLength);
              System.arraycopy(Z,0,WZ,0,NumberLength);
              SubtractBigNbrModN(X, Z, Aux1);
              AddBigNbrModN(TX, TZ, Aux2);
              MontgomeryMult(Aux1, Aux2, W1);
              AddBigNbrModN(X, Z, Aux1);
              SubtractBigNbrModN(TX, TZ, Aux2);
              MontgomeryMult(Aux1, Aux2, W2);
              AddBigNbrModN(W1, W2, Aux1);
              MontgomeryMult(Aux1, Aux1, Aux2);
              MontgomeryMult(Aux2, UZ, X);
              SubtractBigNbrModN(W1, W2, Aux1);
              MontgomeryMult(Aux1, Aux1, Aux2);
              MontgomeryMult(Aux2, UX, Z);
              System.arraycopy(WX,0,UX,0,NumberLength);
              System.arraycopy(WZ,0,UZ,0,NumberLength);
            }
          }                         // end for Q
          if (Pass == 0) {
            if (BigNbrIsZero(GcdAccumulated)) { // If GcdAccumulated is
              System.arraycopy(Xaux,0,X,0,NumberLength);
              System.arraycopy(Zaux,0,Z,0,NumberLength);
              continue;                        // multiple of TestNbr, continue.
            }
            GcdBigNbr(GcdAccumulated, TestNbr, GD);
            if (BigNbrAreEqual(GD, TestNbr) == true) {
              break;
            }
            if (BigNbrAreEqual(GD, BigNbr1) == false) {
              break new_curve;
            }
            break;
          }
        }                        /* end for Pass */
      } while (true);            /* end curve calculation */
    } while (BigNbrAreEqual(GD, TestNbr) == true);

    return BigIntToBigNbr(GD);
  }




  /**********************/
  /* Auxiliary routines */
  /**********************/

static void GenerateSieve(int initial, byte [] sieve, byte [] sieve2310,
                          int [] SmallPrime) {
  int i, j, Q;
  for (i=0; i<23100; i+=2310) {
    System.arraycopy(sieve2310,0,sieve,i,2310);
    }
  j = 5; Q = 13;                /* Point to prime 13 */
  do {
    if (initial>Q*Q) {
      for (i=(int)(((long)initial*((Q-1)/2))%Q); i<23100; i+=Q) {
        sieve[i] = 1;           /* Composite */
        }
      }
    else {
      i=Q*Q-initial;
      if (i<46200) {
        for (i=i/2; i<23100; i+=Q) {
          sieve[i] = 1;         /* Composite */
          }
        }
      else {
        break;
        }
      }
    Q = SmallPrime[++j];
    } while (Q < 5000);
  }


  void MontgomeryMult(long Nbr1[], long Nbr2[], long Prod[]) {
    int i,j;
    long MaxUInt = this.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;
    int MontgomeryMultN = (int)this.MontgomeryMultN;
    long TestNbr[] = this.TestNbr;
    int NumberLength = this.NumberLength;
    int NumberLength1 = NumberLength-1;
    String EcmInfo;

    TestNbr0 = TestNbr[0];
    TestNbr1 = TestNbr[1];
    Nbr2_0 = Nbr2[0];
    Nbr2_1 = Nbr2[1];
    switch (NumberLength) {
      case 2:
        Prod0 = Prod1 = 0;
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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 * MontgomeryMultN) & 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 * MontgomeryMultN) & 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 * MontgomeryMultN) & 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 * MontgomeryMultN) & 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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = ((Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      break;

      case 7:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      i = 0;
      do {

        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prod6 = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<7);
      if (Prod6 > TestNbr6 || Prod6 == TestNbr6 &&
        (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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = (Pr = (Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
        Prod6 = ((Pr>>32) + Prod6 - TestNbr6) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      break;

      case 8:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = Prod7 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      TestNbr7 = TestNbr[7];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      Nbr2_7 = Nbr2[7];
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_7 + Prod7;
        Prod6 = (Pr = (Pr >>> 32) + MontDig*TestNbr7 + (Prd & MaxUInt)) & MaxUInt;
        Prod7 = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<8);
      if (Prod7 > TestNbr7 || Prod7 == TestNbr7 &&
        (Prod6 > TestNbr6 || Prod6 == TestNbr6 &&
          (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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = (Pr = (Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
        Prod6 = (Pr = (Pr>>32) + Prod6 - TestNbr6) & MaxUInt;
        Prod7 = ((Pr>>32) + Prod7 - TestNbr7) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      Prod[7] = Prod7;
      break;

      case 9:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = Prod7 =
        Prod8 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      TestNbr7 = TestNbr[7];
      TestNbr8 = TestNbr[8];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      Nbr2_7 = Nbr2[7];
      Nbr2_8 = Nbr2[8];
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_7 + Prod7;
        Prod6 = (Pr = (Pr >>> 32) + MontDig*TestNbr7 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_8 + Prod8;
        Prod7 = (Pr = (Pr >>> 32) + MontDig*TestNbr8 + (Prd & MaxUInt)) & MaxUInt;
        Prod8 = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<9);
      if (Prod8 > TestNbr8 || Prod8 == TestNbr8 &&
        (Prod7 > TestNbr7 || Prod7 == TestNbr7 &&
          (Prod6 > TestNbr6 || Prod6 == TestNbr6 &&
            (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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = (Pr = (Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
        Prod6 = (Pr = (Pr>>32) + Prod6 - TestNbr6) & MaxUInt;
        Prod7 = (Pr = (Pr>>32) + Prod7 - TestNbr7) & MaxUInt;
        Prod8 = ((Pr>>32) + Prod8 - TestNbr8) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      Prod[7] = Prod7;
      Prod[8] = Prod8;
      break;

      case 10:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = Prod7 =
        Prod8 = Prod9 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      TestNbr7 = TestNbr[7];
      TestNbr8 = TestNbr[8];
      TestNbr9 = TestNbr[9];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      Nbr2_7 = Nbr2[7];
      Nbr2_8 = Nbr2[8];
      Nbr2_9 = Nbr2[9];
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_7 + Prod7;
        Prod6 = (Pr = (Pr >>> 32) + MontDig*TestNbr7 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_8 + Prod8;
        Prod7 = (Pr = (Pr >>> 32) + MontDig*TestNbr8 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_9 + Prod9;
        Prod8 = (Pr = (Pr >>> 32) + MontDig*TestNbr9 + (Prd & MaxUInt)) & MaxUInt;
        Prod9 = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<10);
      if (Prod9 > TestNbr9 || Prod9 == TestNbr9 &&
        (Prod8 > TestNbr8 || Prod8 == TestNbr8 &&
          (Prod7 > TestNbr7 || Prod7 == TestNbr7 &&
            (Prod6 > TestNbr6 || Prod6 == TestNbr6 &&
              (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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = (Pr = (Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
        Prod6 = (Pr = (Pr>>32) + Prod6 - TestNbr6) & MaxUInt;
        Prod7 = (Pr = (Pr>>32) + Prod7 - TestNbr7) & MaxUInt;
        Prod8 = (Pr = (Pr>>32) + Prod8 - TestNbr8) & MaxUInt;
        Prod9 = ((Pr>>32) + Prod9 - TestNbr9) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      Prod[7] = Prod7;
      Prod[8] = Prod8;
      Prod[9] = Prod9;
      break;

      case 11:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = Prod7 =
        Prod8 = Prod9 = Prod10 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      TestNbr7 = TestNbr[7];
      TestNbr8 = TestNbr[8];
      TestNbr9 = TestNbr[9];
      TestNbr10 = TestNbr[10];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      Nbr2_7 = Nbr2[7];
      Nbr2_8 = Nbr2[8];
      Nbr2_9 = Nbr2[9];
      Nbr2_10 = Nbr2[10];
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_7 + Prod7;
        Prod6 = (Pr = (Pr >>> 32) + MontDig*TestNbr7 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_8 + Prod8;
        Prod7 = (Pr = (Pr >>> 32) + MontDig*TestNbr8 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_9 + Prod9;
        Prod8 = (Pr = (Pr >>> 32) + MontDig*TestNbr9 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_10 + Prod10;
        Prod9 = (Pr = (Pr >>> 32) + MontDig*TestNbr10 + (Prd & MaxUInt)) & MaxUInt;
        Prod10 = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<11);
      if (Prod10 > TestNbr10 || Prod10 == TestNbr10 &&
        (Prod9 > TestNbr9 || Prod9 == TestNbr9 &&
          (Prod8 > TestNbr8 || Prod8 == TestNbr8 &&
            (Prod7 > TestNbr7 || Prod7 == TestNbr7 &&
              (Prod6 > TestNbr6 || Prod6 == TestNbr6 &&
                (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) & MaxUInt;
        Prod4 = (Pr = (Pr>>32) + Prod4 - TestNbr4) & MaxUInt;
        Prod5 = (Pr = (Pr>>32) + Prod5 - TestNbr5) & MaxUInt;
        Prod6 = (Pr = (Pr>>32) + Prod6 - TestNbr6) & MaxUInt;
        Prod7 = (Pr = (Pr>>32) + Prod7 - TestNbr7) & MaxUInt;
        Prod8 = (Pr = (Pr>>32) + Prod8 - TestNbr8) & MaxUInt;
        Prod9 = (Pr = (Pr>>32) + Prod9 - TestNbr9) & MaxUInt;
        Prod10 = ((Pr>>32) + Prod10 - TestNbr10) & MaxUInt;
      }
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      Prod[7] = Prod7;
      Prod[8] = Prod8;
      Prod[9] = Prod9;
      Prod[10] = Prod10;
      break;

      default:
        Prod0 = Prod1 = Prod2 = Prod3 = Prod4 = Prod5 = Prod6 = Prod7 =
        Prod8 = Prod9 = Prod10 = 0;
      TestNbr2 = TestNbr[2];
      TestNbr3 = TestNbr[3];
      TestNbr4 = TestNbr[4];
      TestNbr5 = TestNbr[5];
      TestNbr6 = TestNbr[6];
      TestNbr7 = TestNbr[7];
      TestNbr8 = TestNbr[8];
      TestNbr9 = TestNbr[9];
      TestNbr10 = TestNbr[10];
      Nbr2_2 = Nbr2[2];
      Nbr2_3 = Nbr2[3];
      Nbr2_4 = Nbr2[4];
      Nbr2_5 = Nbr2[5];
      Nbr2_6 = Nbr2[6];
      Nbr2_7 = Nbr2[7];
      Nbr2_8 = Nbr2[8];
      Nbr2_9 = Nbr2[9];
      Nbr2_10 = Nbr2[10];
      for (j=11; j<NumberLength; j++) {
        Prod[j] = 0;
      }
      i = 0;
      do {
        ProdB0 = (Nbr = Nbr1[i])*Nbr2_0 + Prod0;
        MontDig = ((int)ProdB0 * MontgomeryMultN) & 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;
        Prd = (Prd >>> 32) + Nbr*Nbr2_6 + Prod6;
        Prod5 = (Pr = (Pr >>> 32) + MontDig*TestNbr6 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_7 + Prod7;
        Prod6 = (Pr = (Pr >>> 32) + MontDig*TestNbr7 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_8 + Prod8;
        Prod7 = (Pr = (Pr >>> 32) + MontDig*TestNbr8 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_9 + Prod9;
        Prod8 = (Pr = (Pr >>> 32) + MontDig*TestNbr9 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2_10 + Prod10;
        Prod9 = (Pr = (Pr >>> 32) + MontDig*TestNbr10 + (Prd & MaxUInt)) & MaxUInt;
        Prd = (Prd >>> 32) + Nbr*Nbr2[11] + Prod[11];
        Prod10 = (Pr = (Pr >>> 32) + MontDig*TestNbr[11] + (Prd & MaxUInt)) & MaxUInt;
        for (j=12; j < NumberLength; j++) {
//          Prd = (Prd >>> 32) + Nbr*Nbr2[j] + Prod[j];
//          Prod[j-1] = (Pr = (Pr >>> 32) + MontDig*TestNbr[j] + (Prd & MaxUInt)) & MaxUInt;
          Prd = (Prd >>> 32) + Nbr*Nbr2[j] + Prod[j];
          Prod[j-1] = (Pr = (Pr >>> 32) + MontDig*TestNbr[j] + (Prd & 0xffffffffL)) & 0xffffffffL;
        }
        Prod[j-1] = (Prd >>> 32) + (Pr >>> 32);
        i++;
      } while (i<NumberLength);
      Prod[0] = Prod0;
      Prod[1] = Prod1;
      Prod[2] = Prod2;
      Prod[3] = Prod3;
      Prod[4] = Prod4;
      Prod[5] = Prod5;
      Prod[6] = Prod6;
      Prod[7] = Prod7;
      Prod[8] = Prod8;
      Prod[9] = Prod9;
      Prod[10] = Prod10;
      for (j=NumberLength-1; j>=0; j--) {
        if (Prod[j] != TestNbr[j]) {
          break;
        }
      }
      if (j<0 || j>=0 && Prod[j] >= TestNbr[j]) {
        Pr = 0;
        for (j=0; j<NumberLength; j++) {
          Prod[j] = (Pr = (Pr>>32) + Prod[j] - TestNbr[j]) & MaxUInt;
        }
      }
    }                      /* end switch */
  }




  void AddBigNbrModN(long Nbr1[], long Nbr2[], long Sum[]) {
    int NumberLength = this.NumberLength;
    long Cy=0;
    int i;

    for(i=0; i<NumberLength; i++) {
      Cy += Nbr1[i]+Nbr2[i]-TestNbr[i];
      Sum[i] = Cy & MaxUInt;
      Cy >>= 32;
    }
    if (Cy < 0) {
      Cy = 0;
      for(i=0; i<NumberLength; i++) {
        Cy += Sum[i]+TestNbr[i];
        Sum[i] = Cy & MaxUInt;
        Cy >>= 32;
      }
    }
  }


  void SubtractBigNbrModN(long Nbr1[], long Nbr2[], long Diff[]) {
    int NumberLength = this.NumberLength;
    long Cy=0;
    int i;

    for(i=0; i<NumberLength; i++) {
      Cy += Nbr1[i]-Nbr2[i];
      Diff[i] = Cy & MaxUInt;
      Cy >>= 32;
    }
    if (Cy < 0) {
      Cy = 0;
      for(i=0; i<NumberLength; i++) {
        Cy += Diff[i]+TestNbr[i];
        Diff[i] = Cy & MaxUInt;
        Cy >>= 32;
      }
    }
  }




  BigInteger BigIntToBigNbr(long [] GD) {
    byte [] Result;
    int I;
    int NL;

    NL = NumberLength*4;
    Result = new byte[NL];
    for (I=0;I<NumberLength;I++) {
      Result[NL-1-4*I] = (byte)(GD[I] & 0xFF);
      Result[NL-2-4*I] = (byte)(GD[I]/0x100 & 0xFF);
      Result[NL-3-4*I] = (byte)(GD[I]/0x10000 & 0xFF);
      Result[NL-4-4*I] = (byte)(GD[I]/0x1000000 & 0xFF);
    }
    return (new BigInteger(Result));
  }


  boolean BigNbrAreEqual(long Nbr1[], long Nbr2[]) {
    for(int i=0; i<NumberLength; i++) {
      if (Nbr1[i] != Nbr2[i]) {
        return false;
      }
    }
    return true;
  }

  // Gcd calculation:
  // Step 1: Set k<-0, and then repeatedly set k<-k+1, u<-u/2, v<-v/2
  //         zero or more times until u and v are not both even.
  // Step 2: If u is odd, set t<-(-v) and go to step 4. Otherwise set t<-u.
  // Step 3: Set t<-t/2
  // Step 4: If t is even, go back to step 3.
  // Step 5: If t>0, set u<-t, otherwise set v<-(-t).
  // Step 6: Set t<-u-v. If t!=0, go back to step 3.
  // Step 7: The GCD is u*2^k.
  void GcdBigNbr(long Nbr1[], long Nbr2[], long Gcd[]) {
    int i,k;
    int NumberLength = this.NumberLength;

    System.arraycopy(Nbr1,0,CalcAuxGcdU,0,NumberLength);
    System.arraycopy(Nbr2,0,CalcAuxGcdV,0,NumberLength);
    for (i=0; i<NumberLength; i++) {
      if (CalcAuxGcdU[i] != 0) {break;}
    }
    if (i == NumberLength) {
      System.arraycopy(CalcAuxGcdV,0,Gcd,0,NumberLength);
      return;
    }
    for (i=0; i<NumberLength; i++) {
      if (CalcAuxGcdV[i] != 0) {break;}
    }
    if (i == NumberLength) {
      System.arraycopy(CalcAuxGcdU,0,Gcd,0,NumberLength);
      return;
    }
    if (CalcAuxGcdU[NumberLength-1] >= DosALa31) {
      ChSignBigNbr(CalcAuxGcdU);
    }
    if (CalcAuxGcdV[NumberLength-1] >= DosALa31) {
      ChSignBigNbr(CalcAuxGcdV);
    }
    k=0;
    while ((CalcAuxGcdU[0] & 1)==0 && (CalcAuxGcdV[0] & 1)==0) {   // Step 1
      k++;
      DivBigNbrByLong(CalcAuxGcdU,2,CalcAuxGcdU);
      DivBigNbrByLong(CalcAuxGcdV,2,CalcAuxGcdV);
    }
    if ((CalcAuxGcdU[0] & 1)==1) {                                 // Step 2
      System.arraycopy(CalcAuxGcdV,0,CalcAuxGcdT,0,NumberLength);
      ChSignBigNbr(CalcAuxGcdT);
    }
    else {
      System.arraycopy(CalcAuxGcdU,0,CalcAuxGcdT,0,NumberLength);
    }
    do {
      while ((CalcAuxGcdT[0] & 1)==0) {                            // Step 4
        DivBigNbrByLong(CalcAuxGcdT,2,CalcAuxGcdT);                // Step 3
      }
      if (CalcAuxGcdT[NumberLength-1] < DosALa31) {                             // Step 5
        System.arraycopy(CalcAuxGcdT,0,CalcAuxGcdU,0,NumberLength);
      }
      else {
        System.arraycopy(CalcAuxGcdT,0,CalcAuxGcdV,0,NumberLength);
        ChSignBigNbr(CalcAuxGcdV);
      }
      SubtractBigNbr(CalcAuxGcdU,CalcAuxGcdV,CalcAuxGcdT);         // Step 6
      for (i=0; i<NumberLength; i++) {
        if (CalcAuxGcdT[i] != 0) {break;}
      }
    } while (i != NumberLength);
    System.arraycopy(CalcAuxGcdU,0,Gcd,0,NumberLength);            // Step 7
    while (k>0) {
      AddBigNbr(Gcd,Gcd,Gcd);
      k--;
    }
  }




  boolean BigNbrIsZero(long Nbr[]) {
    for(int i=0; i<NumberLength; i++) {
      if (Nbr[i] != 0) {
        return false;
      }
    }
    return true;
  }


  void AddBigNbr(long Nbr1[], long Nbr2[], long Sum[]) {
    int NumberLength = this.NumberLength;
    long Cy=0;
    for(int i=0; i<NumberLength; i++) {
      Cy += Nbr1[i]+Nbr2[i];
      Sum[i] = Cy & MaxUInt;
      Cy >>= 32;
    }
  }


  void SubtractBigNbr(long Nbr1[], long Nbr2[], long Diff[]) {
    int NumberLength = this.NumberLength;
    long Cy=0;
    for(int i=0; i<NumberLength; i++) {
      Cy += Nbr1[i]-Nbr2[i];
      Diff[i] = Cy & MaxUInt;
      Cy >>= 32;
    }
  }


  void ChSignBigNbr(long Nbr[]) {
    int NumberLength = this.NumberLength;
    long Cy=0;
    for(int i=0; i<NumberLength; i++) {
      Cy -= Nbr[i];
      Nbr[i]=(Cy>=0?Cy:Cy+DosALa32);
      Cy=(Cy>=0?0:-1);
    }
  }



  void DivBigNbrByLong(long Dividend[], long Divisor, long Quotient[]) {
    int i;
    boolean ChSignDivisor=false;
    long Divid, Rem = 0;

    if (Divisor < 0) {
      ChSignDivisor = true;
      Divisor = -Divisor;
    }
    if (Dividend[NumberLength-1] >= DosALa31) {
      Rem = Divisor - 1;
    }
    for (i=NumberLength-1; i>=0; i--) {
      Divid = Dividend[i] + (Rem << 32);
      Rem = Divid % Divisor;
      Quotient[i] = Divid / Divisor;
    }
    if (ChSignDivisor) {
      ChSignBigNbr(Quotient);
    }
  }



  /***********************************************************************/
  /* NAME: ModInvBigNbr                                                  */
  /*                                                                     */
  /* PURPOSE: Find the inverse multiplicative modulo v.                  */
  /*                                                                     */
  /* The algorithm terminates with u1 = u^(-1) MOD v.                    */
  /***********************************************************************/
  void ModInvBigNbr(long[] a, long [] inv, long[] b) {
    int i, j;
    int NumberLength = this.NumberLength;
    int Dif, E;
    int st1, st2;
    long Yaa, Yab;     // 2^E * A'     = Yaa A + Yab B
    long Yba, Ybb;     // 2^E * B'     = Yba A + Ybb B
    long Ygb0;         // 2^E * Mu'    = Yaa Mu + Yab Gamma + Ymb0 B0
    long Ymb0;         // 2^E * Gamma' = Yba Mu + Ybb Gamma + Ygb0 B0
    int Iaa, Iab, Iba, Ibb;
    long Tmp1, Tmp2, Tmp3, Tmp4, Tmp5;
    int B0l = (int)b[0];
    int invB0l;
    int Al, Bl, T1, Gl, Ml, P;
    long T;
    long Cy1, Cy2, Cy3, Cy4;
    int Yaah, Yabh, Ybah, Ybbh;
    int Ymb0h, Ygb0h;
    long Pr1, Pr2, Pr3, Pr4, Pr5, Pr6, Pr7;
    long [] CalcAuxModInvA = this.CalcAuxGcdU;
    long [] CalcAuxModInvB = this.CalcAuxGcdV;
    long [] CalcAuxModInvMu = inv;
    long [] CalcAuxModInvGamma = this.CalcAuxGcdT;

    invB0l = B0l;            // 2 least significant bits of inverse correct.
    invB0l = invB0l*(2-B0l*invB0l); // 4 LSB of inverse correct.
    invB0l = invB0l*(2-B0l*invB0l); // 8 LSB of inverse correct.
    invB0l = invB0l*(2-B0l*invB0l); // 16 LSB of inverse correct.
    invB0l = invB0l*(2-B0l*invB0l); // 32 LSB of inverse correct.

    for (i=NumberLength-1; i>=0; i--) {
      CalcAuxModInvA[i] = a[i];
      CalcAuxModInvB[i] = b[i];
      CalcAuxModInvGamma[i] = 0;
      CalcAuxModInvMu[i] = 0;
    }
    CalcAuxModInvMu[0] = 1;
    Dif = 0;
outer_loop:
    do {
      Iaa = Ibb = 1;
      Iab = Iba = 0;
      Al = (int)CalcAuxModInvA[0];
      Bl = (int)CalcAuxModInvB[0];
      E = 0; P = 1;
      if (Bl == 0) {
        for (i=NumberLength-1; i>=0; i--) {
          if (CalcAuxModInvB[i] != 0) break;
          }
        if (i<0) break;        // Go out of loop if CalcAuxModInvB = 0
        }
      do {
        T1 = 0;
        while ((Bl & 1) == 0) {
          if (E == 31) {
            Yaa = Iaa; Yab = Iab; Yba = Iba; Ybb = Ibb;
            Gl = (int)CalcAuxModInvGamma[0];
            Ml = (int)CalcAuxModInvMu[0];
            Dif++;
            T1++;
            Yaa <<= T1; Yab <<= T1;
            Ymb0 = (-(int)Yaa*Ml-(int)Yab*Gl)*invB0l;
            Ygb0 = (-Iba*Ml-Ibb*Gl)*invB0l;
            Cy1 = Cy2 = Cy3 = Cy4 = 0;
            Yaah = (int)(Yaa >> 32); Yabh = (int)(Yab >> 32);
            Ybah = (int)(Yba >> 32); Ybbh = (int)(Ybb >> 32);
            Ymb0h = (int)(Ymb0 >> 32);
            Ygb0h = (int)(Ygb0 >> 32);
            Yaa &= 0xFFFFFFFFL; Yab &= 0xFFFFFFFFL;
            Yba &= 0xFFFFFFFFL; Ybb &= 0xFFFFFFFFL;
            Ymb0 &= 0xFFFFFFFFL;
            Ygb0 &= 0xFFFFFFFFL;
      
            st1 = Yaah*6+Yabh*2+Ymb0h;
            st2 = Ybah*6+Ybbh*2+Ygb0h;
            for (i=0; i<NumberLength; i++) {
              Pr1 = Yaa * (Tmp1 = CalcAuxModInvMu[i]);
              Pr2 = Yab * (Tmp2 = CalcAuxModInvGamma[i]);
              Pr3 = Ymb0 * (Tmp3 = b[i]);
              Pr4 = (Pr1 & 0xFFFFFFFFL) + (Pr2 & 0xFFFFFFFFL) + (Pr3 & 0xFFFFFFFFL) + Cy3;
              Pr5 = Yaa * (Tmp4 = CalcAuxModInvA[i]);
              Pr6 = Yab * (Tmp5 = CalcAuxModInvB[i]);
              Pr7 = (Pr5 & 0xFFFFFFFFL) + (Pr6 & 0xFFFFFFFFL) + Cy1;
              switch (st1) {
                case  -9: Cy3 = -Tmp1-Tmp2-Tmp3; Cy1 = -Tmp4-Tmp5; break;
                case  -8: Cy3 = -Tmp1-Tmp2     ; Cy1 = -Tmp4-Tmp5; break;
                case  -7: Cy3 = -Tmp1     -Tmp3; Cy1 = -Tmp4     ; break;
                case  -6: Cy3 = -Tmp1          ; Cy1 = -Tmp4     ; break;
                case  -5: Cy3 = -Tmp1+Tmp2-Tmp3; Cy1 = -Tmp4+Tmp5; break;
                case  -4: Cy3 = -Tmp1+Tmp2     ; Cy1 = -Tmp4+Tmp5; break;
                case  -3: Cy3 =      -Tmp2-Tmp3; Cy1 =      -Tmp5; break;
                case  -2: Cy3 =      -Tmp2     ; Cy1 =      -Tmp5; break;
                case  -1: Cy3 =           -Tmp3; Cy1 =          0; break;
                case   0: Cy3 =               0; Cy1 =          0; break;
                case   1: Cy3 =       Tmp2-Tmp3; Cy1 =       Tmp5; break;
                case   2: Cy3 =       Tmp2     ; Cy1 =       Tmp5; break;
                case   3: Cy3 =  Tmp1-Tmp2-Tmp3; Cy1 =  Tmp4-Tmp5; break;
                case   4: Cy3 =  Tmp1-Tmp2     ; Cy1 =  Tmp4-Tmp5; break;
                case   5: Cy3 =  Tmp1     -Tmp3; Cy1 =  Tmp4     ; break;
                case   6: Cy3 =  Tmp1          ; Cy1 =  Tmp4     ; break;
                case   7: Cy3 =  Tmp1+Tmp2-Tmp3; Cy1 =  Tmp4+Tmp5; break;
                case   8: Cy3 =  Tmp1+Tmp2     ; Cy1 =  Tmp4+Tmp5; break;
                }
              Cy3 += (Pr1 >>> 32) + (Pr2 >>> 32) + (Pr3 >>> 32) + (Pr4 >> 32);
              Cy1 += (Pr5 >>> 32) + (Pr6 >>> 32) + (Pr7 >> 32);
              if (i>0) {
                CalcAuxModInvMu[i-1] = Pr4 & 0xFFFFFFFFL; 
                CalcAuxModInvA[i-1] = Pr7 & 0xFFFFFFFFL; 
                }
              Pr1 = Yba * Tmp1;
              Pr2 = Ybb * Tmp2;
              Pr3 = Ygb0 * Tmp3;
              Pr4 = (Pr1 & 0xFFFFFFFFL) + (Pr2 & 0xFFFFFFFFL) + (Pr3 & 0xFFFFFFFFL) + Cy4;
              Pr5 = Yba * Tmp4;
              Pr6 = Ybb * Tmp5;
              Pr7 = (Pr5 & 0xFFFFFFFFL) + (Pr6 & 0xFFFFFFFFL) + Cy2;
              switch (st2) {
                case  -9: Cy4 = -Tmp1-Tmp2-Tmp3; Cy2 = -Tmp4-Tmp5; break;
                case  -8: Cy4 = -Tmp1-Tmp2     ; Cy2 = -Tmp4-Tmp5; break;
                case  -7: Cy4 = -Tmp1     -Tmp3; Cy2 = -Tmp4     ; break;
                case  -6: Cy4 = -Tmp1          ; Cy2 = -Tmp4     ; break;
                case  -5: Cy4 = -Tmp1+Tmp2-Tmp3; Cy2 = -Tmp4+Tmp5; break;
                case  -4: Cy4 = -Tmp1+Tmp2     ; Cy2 = -Tmp4+Tmp5; break;
                case  -3: Cy4 =      -Tmp2-Tmp3; Cy2 =      -Tmp5; break;
                case  -2: Cy4 =      -Tmp2     ; Cy2 =      -Tmp5; break;
                case  -1: Cy4 =           -Tmp3; Cy2 =          0; break;
                case   0: Cy4 =               0; Cy2 =          0; break;
                case   1: Cy4 =       Tmp2-Tmp3; Cy2 =       Tmp5; break;
                case   2: Cy4 =       Tmp2     ; Cy2 =       Tmp5; break;
                case   3: Cy4 =  Tmp1-Tmp2-Tmp3; Cy2 =  Tmp4-Tmp5; break;
                case   4: Cy4 =  Tmp1-Tmp2     ; Cy2 =  Tmp4-Tmp5; break;
                case   5: Cy4 =  Tmp1     -Tmp3; Cy2 =  Tmp4     ; break;
                case   6: Cy4 =  Tmp1          ; Cy2 =  Tmp4     ; break;
                case   7: Cy4 =  Tmp1+Tmp2-Tmp3; Cy2 =  Tmp4+Tmp5; break;
                case   8: Cy4 =  Tmp1+Tmp2     ; Cy2 =  Tmp4+Tmp5; break;
                }
              Cy4 += (Pr1 >>> 32) + (Pr2 >>> 32) + (Pr3 >>> 32) + (Pr4 >> 32);
              Cy2 += (Pr5 >>> 32) + (Pr6 >>> 32) + (Pr7 >> 32);
              if (i>0) {
                CalcAuxModInvGamma[i-1] = Pr4 & 0xFFFFFFFFL; 
                CalcAuxModInvB[i-1] = Pr7 & 0xFFFFFFFFL; 
                }
              }

            if ((int)CalcAuxModInvA[i-1] < 0) {
              Cy1 -= Yaa;
              Cy2 -= Yba;
              }
            if ((int)CalcAuxModInvB[i-1] < 0) {
              Cy1 -= Yab;
              Cy2 -= Ybb;
              }
            if ((int)CalcAuxModInvMu[i-1] < 0) {
              Cy3 -= Yaa;
              Cy4 -= Yba;
              }
            if ((int)CalcAuxModInvGamma[i-1] < 0) {
              Cy3 -= Yab;
              Cy4 -= Ybb;
              }
            CalcAuxModInvA[i-1] = Cy1 & 0xFFFFFFFFL;
            CalcAuxModInvB[i-1] = Cy2 & 0xFFFFFFFFL;
            CalcAuxModInvMu[i-1] = Cy3 & 0xFFFFFFFFL;
            CalcAuxModInvGamma[i-1] = Cy4 & 0xFFFFFFFFL;
            continue outer_loop;
            }
          Bl >>= 1; Dif++;
          E++; P*=2; T1++;
          };            /* end while */
        Iaa <<= T1; Iab <<= T1;
        if (Dif >= 0) {
          Dif = -Dif;
          if (((Al+Bl) & 3) == 0) {
            T1 = Iba; Iba += Iaa; Iaa = T1;
            T1 = Ibb; Ibb += Iab; Iab = T1;
            T1 = Bl; Bl += Al; Al = T1;
            }
          else {
            T1 = Iba; Iba -= Iaa; Iaa = T1;
            T1 = Ibb; Ibb -= Iab; Iab = T1;
            T1 = Bl; Bl -= Al; Al = T1;
            }
          }
        else {
          if (((Al+Bl) & 3) == 0) {
            Iba += Iaa; Ibb += Iab; Bl += Al;
            }
          else {
            Iba -= Iaa; Ibb -= Iab; Bl -= Al;
            }
          }
        Dif--;
        } while (true);
      } while (true);
    if (CalcAuxModInvA[0] != 1) {
      SubtractBigNbr(b, CalcAuxModInvMu, CalcAuxModInvMu);
      }
    if ((int)CalcAuxModInvMu[i = NumberLength-1]<0) {
      AddBigNbr(b, CalcAuxModInvMu, CalcAuxModInvMu);
      }
    for (; i>=0; i--) {
      if (b[i] != CalcAuxModInvMu[i]) break;
      }
    if (i<0 || b[i] < CalcAuxModInvMu[i]) {                 // If B < Mu
      SubtractBigNbr(CalcAuxModInvMu, b, CalcAuxModInvMu);  // Mu <- Mu - B
      }
    System.arraycopy(CalcAuxModInvMu, 0, inv, 0, NumberLength);
    }



/* computes nP from P=(x:z) and puts the result in (x:z). Assumes n>2. */
void prac(int n, long [] x, long [] z,
          long [] xT, long [] zT, long [] xT2, long [] zT2)
{
   int d,e,r,i;
   long [] t;
   long [] xA = x, zA = z;
   long [] xB = fieldAux1, zB = fieldAux2;
   long [] xC = fieldAux3, zC = fieldAux4;
   double v[] = 
     {1.61803398875,1.72360679775,1.618347119656,1.617914406529,1.612429949509,
    1.632839806089,1.620181980807,1.580178728295,1.617214616534,1.38196601125};
   
   /* chooses the best value of v */
   r=lucas_cost(n,v[0]);i=0;
   for (d=1;d<10;d++) {
     e=lucas_cost(n,v[d]);
     if (e<r) { r=e; i=d; }
   }
   d=n;
   r=(int)((double)d/v[i]+0.5);
   /* first iteration always begins by Condition 3, then a swap */
   d=n-r; e=2*r-n;
   System.arraycopy(xA,0,xB,0,NumberLength);     // B = A
   System.arraycopy(zA,0,zB,0,NumberLength);
   System.arraycopy(xA,0,xC,0,NumberLength);     // C = A
   System.arraycopy(zA,0,zC,0,NumberLength);
   duplicate(xA,zA,xA,zA); /* A=2*A */
   while (d!=e) {
         if (d<e) { r=d; d=e; e=r; t=xA; xA=xB; xB=t; t=zA; zA=zB; zB=t; }
	 /* do the first line of Table 4 whose condition qualifies */
	 if (4*d<=5*e && ((d+e)%3)==0) { /* condition 1 */
	    r=(2*d-e)/3; e=(2*e-d)/3; d=r;
	    add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T = f(A,B,C) */
	    add3(xT2,zT2,xT,zT,xA,zA,xB,zB); /* T2 = f(T,A,B) */
	    add3(xB,zB,xB,zB,xT,zT,xA,zA); /* B = f(B,T,A) */
	    t=xA; xA=xT2; xT2=t; t=zA; zA=zT2; zT2=t; /* swap A and T2 */
	  } else
	 if (4*d<=5*e && (d-e)%6==0) { /* condition 2 */
	   d=(d-e)/2; 
	   add3(xB,zB,xA,zA,xB,zB,xC,zC); /* B = f(A,B,C) */
	   duplicate(xA,zA,xA,zA); /* A = 2*A */
	 } else
	 if (d<=(4*e)) { /* condition 3 */
	   d-=e; 
	   add3(xT,zT,xB,zB,xA,zA,xC,zC); /* T = f(B,A,C) */
	   t=xB; xB=xT; xT=xC; xC=t;
	   t=zB; zB=zT; zT=zC; zC=t; /* circular permutation (B,T,C) */
	 } else
	 if ((d+e)%2==0) { /* condition 4 */
	   d=(d-e)/2; 
	   add3(xB,zB,xB,zB,xA,zA,xC,zC); /* B = f(B,A,C) */
	   duplicate(xA,zA,xA,zA); /* A = 2*A */
	 } else
	 if (d%2==0) { /* condition 5 */
	   d/=2; 
	   add3(xC,zC,xC,zC,xA,zA,xB,zB); /* C = f(C,A,B) */
	   duplicate(xA,zA,xA,zA); /* A = 2*A */
	 } else
	 if (d%3==0) { /* condition 6 */
	   d=d/3-e; 
	   duplicate(xT,zT,xA,zA); /* T1 = 2*A */
	   add3(xT2,zT2,xA,zA,xB,zB,xC,zC); /* T2 = f(A,B,C) */
	   add3(xA,zA,xT,zT,xA,zA,xA,zA); /* A = f(T1,A,A) */
	   add3(xT,zT,xT,zT,xT2,zT2,xC,zC); /* T1 = f(T1,T2,C) */
	   t=xC; xC=xB; xB=xT; xT=t;
	   t=zC; zC=zB; zB=zT; zT=t; /* circular permutation (C,B,T) */
	 } else
	 if ((d+e)%3==0) { /* condition 7 */
	   d=(d-2*e)/3; 
	   add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T1 = f(A,B,C) */
	   add3(xB,zB,xT,zT,xA,zA,xB,zB); /* B = f(T1,A,B) */
	   duplicate(xT,zT,xA,zA); add3(xA,zA,xA,zA,xT,zT,xA,zA); /* A = 3*A */
	 } else
	 if ((d-e)%3==0) { /* condition 8 */
	   d=(d-e)/3; 
	   add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T1 = f(A,B,C) */
	   add3(xC,zC,xC,zC,xA,zA,xB,zB); /* C = f(A,C,B) */
	   t=xB; xB=xT; xT=t; t=zB; zB=zT; zT=t; /* swap B and T */
	   duplicate(xT,zT,xA,zA);
	   add3(xA,zA,xA,zA,xT,zT,xA,zA); /* A = 3*A */
	 } else
	 if (e%2==0) { /* condition 9 */
	   e/=2; 
	   add3(xC,zC,xC,zC,xB,zB,xA,zA); /* C = f(C,B,A) */
	   duplicate(xB,zB,xB,zB); /* B = 2*B */
         }
       }
   add3(x,z,xA,zA,xB,zB,xC,zC);
}


/* adds Q=(x2:z2) and R=(x1:z1) and puts the result in (x3:z3),
     using 5/6 mul, 6 add/sub and 6 mod. One assumes that Q-R=P or R-Q=P where P=(x:z).
     Uses the following global variables:
     - n : number to factor
     - x, z : coordinates of P
     - u, v, w : auxiliary variables
Modifies: x3, z3, u, v, w.
(x3,z3) may be identical to (x2,z2) and to (x,z)
*/
void add3(long [] x3, long [] z3, long [] x2, long [] z2,
          long [] x1, long [] z1, long [] x, long [] z)
{
   long [] t = fieldTX;
   long [] u = fieldTZ;
   long [] v = fieldUX;
   long [] w = fieldUZ;
   SubtractBigNbrModN(x2,z2,v);        // v = x2-z2
   AddBigNbrModN(x1,z1,w);             // w = x1+z1
   MontgomeryMult(v,w,u);              // u = (x2-z2)*(x1+z1)
   AddBigNbrModN(x2,z2,w);             // w = x2+z2
   SubtractBigNbrModN(x1,z1,t);        // t = x1-z1
   MontgomeryMult(t,w,v);              // v = (x2+z2)*(x1-z1)
   AddBigNbrModN(u,v,t);               // t = 2*(x1*x2-z1*z2)
   MontgomeryMult(t,t,w);              // w = 4*(x1*x2-z1*z2)^2
   SubtractBigNbrModN(u,v,t);          // t = 2*(x2*z1-x1*z2)
   MontgomeryMult(t,t,v);              // v = 4*(x2*z1-x1*z2)^2
   if (BigNbrAreEqual(x,x3)) {
     System.arraycopy(x,0,u,0,NumberLength);
     System.arraycopy(w,0,t,0,NumberLength);
     MontgomeryMult(z,t,w);
     MontgomeryMult(v,u,z3);
     System.arraycopy(w,0,x3,0,NumberLength);
   }
   else {
     MontgomeryMult(w,z,x3);           // x3 = 4*z*(x1*x2-z1*z2)^2
     MontgomeryMult(x,v,z3);           // z3 = 4*x*(x2*z1-x1*z2)^2
   }
}



/* computes 2P=(x2:z2) from P=(x1:z1), with 5 mul, 4 add/sub, 5 mod.
     Uses the following global variables:
     - n : number to factor
     - b : (a+2)/4 mod n
     - u, v, w : auxiliary variables
Modifies: x2, z2, u, v, w
*/
void duplicate(long [] x2,long [] z2,long [] x1,long [] z1)
{
   long [] u = fieldUZ;
   long [] v = fieldTX;
   long [] w = fieldTZ;
   AddBigNbrModN(x1,z1,w);             // w = x1+z1
   MontgomeryMult(w,w,u);              // u = (x1+z1)^2
   SubtractBigNbrModN(x1,z1,w);        // w = x1-z1
   MontgomeryMult(w,w,v);              // v = (x1-z1)^2
   MontgomeryMult(u,v,x2);             // x2 = u*v
   SubtractBigNbrModN(u,v,w);          // w = u-v = 4*x1*z1
   MontgomeryMult(fieldAA,w,u);
   AddBigNbrModN(u,v,u);               // u = (v+b*w)
   MontgomeryMult(w,u,z2);             // z2 = (w*u)
}


final static int ADD=6; /* number of multiplications in an addition */
final static int DUP=5; /* number of multiplications in a duplicate */

/* returns the number of modular multiplications */
static int lucas_cost(int n, double v)
{
  int c,d,e,r;

  d=n; r=(int)((double)d/v+0.5);
  if (r>=n) return(ADD*n);
  d=n-r; e=2*r-n; c=DUP+ADD; /* initial duplicate and final addition */
  while (d!=e) {
    if (d<e) { r=d; d=e; e=r; }
    if (4*d<=5*e && ((d+e)%3)==0) { /* condition 1 */
      r=(2*d-e)/3; e=(2*e-d)/3; d=r; c+=3*ADD; /* 3 additions */
    } else
    if (4*d<=5*e && (d-e)%6==0) { /* condition 2 */
      d=(d-e)/2; c+=ADD+DUP; /* one addition, one duplicate */
    } else
    if (d<=(4*e)) { /* condition 3 */
      d-=e; c+=ADD; /* one addition */
    } else
    if ((d+e)%2==0) { /* condition 4 */
      d=(d-e)/2; c+=ADD+DUP; /* one addition, one duplicate */
    } else
    if (d%2==0) { /* condition 5 */
      d/=2; c+=ADD+DUP; /* one addition, one duplicate */
    } else
    if (d%3==0) { /* condition 6 */
      d=d/3-e; c+=3*ADD+DUP; /* three additions, one duplicate */
    } else
    if ((d+e)%3==0) { /* condition 7 */
      d=(d-2*e)/3; c+=3*ADD+DUP; /* three additions, one duplicate */
    } else
    if ((d-e)%3==0) { /* condition 8 */
      d=(d-e)/3; c+=3*ADD+DUP; /* three additions, one duplicate */
    } else 
    if (e%2==0) { /* condition 9 */
      e/=2; c+=ADD+DUP; /* one addition, one duplicate */
    }
  }
  return(c);
}


  void MultBigNbrByLongModN(long Nbr1[], long Nbr2, long Prod[]) {
    int NumberLength = this.NumberLength;
    long Pr;
    int j;

    Pr = 0;
    for (j=0; j<NumberLength; j++) {
      Pr = (Pr >>> 32) + Nbr2*Nbr1[j];
      Prod[j] = Pr & MaxUInt;
    }
    Prod[j] = (Pr >>> 32);
    AdjustModN(Prod);
  }



  void MultBigNbrModN(long Nbr1[], long Nbr2[], long Prod[]) {
    int NumberLength = this.NumberLength;
    int i,j;
    long Pr, Nbr;

    i=NumberLength;
    do {
      Prod[--i] = 0;
    } while (i>0);
    i=NumberLength;
    do {
      Nbr=Nbr1[--i];
      j=NumberLength;
      do {
        Prod[j] = Prod[j-1];
        j--;
      } while (j>0);
      Prod[0] = 0;
      Pr = 0;
      for (j=0; j<NumberLength; j++) {
        Pr = (Pr >>> 32) + Nbr*Nbr2[j] + Prod[j];
        Prod[j] = Pr & MaxUInt;
      }
      Prod[j] += (Pr >>> 32);
      AdjustModN(Prod);
    } while (i>0);
  }



  void AdjustModN(long Nbr[]) {
    int NumberLength = this.NumberLength;
    long TrialQuotient;
    long Cy;
    int i;
    double dAux;

    dAux = (double)Nbr[NumberLength]*dDosALa32 + (double)Nbr[NumberLength-1];
    if (NumberLength > 1) {
      dAux += (double)Nbr[NumberLength-2]/dDosALa32;
    }
    TrialQuotient = (long)Math.ceil(dAux / dN) + 2;
    if (TrialQuotient >= DosALa32) {
      Cy = 0;
      for(i=0; i<NumberLength; i++) {
        Cy = Nbr[i+1]-(TrialQuotient>>>32)*TestNbr[i] - Cy;
        Nbr[i+1] = Cy & MaxUInt;
        Cy = (MaxUInt - Cy) >>> 32;
      }
      TrialQuotient &= MaxUInt;
    }
    Cy = 0;
    for (i=0; i<NumberLength; i++) {
      Cy = Nbr[i] - TrialQuotient*TestNbr[i] - Cy;
      Nbr[i] = Cy & MaxUInt;
      Cy = (MaxUInt - Cy) >>> 32;
    }
    Nbr[NumberLength] -= Cy;
    while ((Nbr[NumberLength] & MaxUInt) != 0) {
      Cy=0;
      for (i=0; i<NumberLength; i++) {
        Cy += Nbr[i]+TestNbr[i];
        Nbr[i] = Cy & MaxUInt;
        Cy >>= 32;
      }
      Nbr[NumberLength] += Cy;
    }
  }



  void LongToBigNbr(long Nbr, long Out[]) {
    int i;

    Out[0] = Nbr & MaxUInt;
    Out[1] = Nbr >>> 32;
    for (i=2; i<NumberLength; i++) {
      Out[i] = (Nbr<0? MaxUInt: 0);
    }
  }



  void BigNbrToBigInt(BigInteger N) {
    byte [] Result;
    int I,J;
    long K,P;

    Result = N.toByteArray();
    NumberLength = (Result.length+3)/4;
    J = 0;
    K = 1;
    P = 0;
    for (I=Result.length-1; I>=0; I--) {
      P += K*(long)(Result[I]>=0?Result[I]:Result[I]+256);
      K *= 0x100;
      if (K == DosALa32) {
        TestNbr[J] = P;
        J++;
        K = 1;
        P = 0;
      }
    }
    TestNbr[J] = P;
    if (TestNbr[NumberLength-1] > Mi) {
      TestNbr[NumberLength] = 0;
      NumberLength++;
    }
    TestNbr[NumberLength] = 0;
  }



  void GetMontgomeryParms() {
    int NumberLength = this.NumberLength;
    int N, x, j;

    dN = (double)TestNbr[NumberLength-1];
    if (NumberLength > 1) {
      dN += (double)TestNbr[NumberLength-2]/dDosALa32;
    }
    if (NumberLength > 2) {
      dN += (double)TestNbr[NumberLength-3]/dDosALa64;
    }

    x = N = (int)TestNbr[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.
    MontgomeryMultN = (-x) & MaxUInt;
    j=NumberLength;
    MontgomeryMultR1[j] = 1;
    do {
      MontgomeryMultR1[--j] = 0;
    } while (j>0);
    AdjustModN(MontgomeryMultR1);
    MultBigNbrModN(MontgomeryMultR1,MontgomeryMultR1,MontgomeryMultR2);
    MontgomeryMult(MontgomeryMultR2,MontgomeryMultR2,MontgomeryMultAfterInv);
    AddBigNbrModN(MontgomeryMultR1,MontgomeryMultR1,MontgomeryMultR2);
  }
}