superfac9.java

[bearnol]

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

// Added JGW 2005-05-05
// Wanless Factorize, Wanless Extended Factorize
// copyright 2000 James Wanless. All rights reserved.

// Added DT 2006-02-15
// Fermat factorization method

// Added JGW 2006-04-25
// Fermat probable primality test

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

public class superfac9 {

// Declare ZERO so we don't have to write BigInteger.ZERO every time,
etc
static final BigInteger ZERO = new BigInteger("0");
static final BigInteger ONE = new BigInteger("1");
static final BigInteger TWO = new BigInteger("2");
static final BigInteger THREE = new BigInteger("3");
static final BigInteger FIVE = new BigInteger("5");
static final BigInteger TEN = new BigInteger("10");
static final BigInteger HUNDRED = new BigInteger("100");
static final BigInteger TENTHOUSAND = new BigInteger("10000");

static BigInteger n = new BigInteger("0");

static String s;

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

static final int PWmax=32, Qmax=30241, LEVELmax=11;
static final int NLen = 1200;

final int aiIndx[] = new int[Qmax];
final int aiF[] = new int[Qmax];

static final int aiP[] = {2, 3, 5, 7, 11, 13};
static final int aiQ[] = {
2, 3, 5, 7, 13, 11, 31, 61, 19, 37, 181,
29, 43, 71, 127, 211, 421, 631, 41, 73, 281,
2521, 17, 113, 241, 337, 1009, 109, 271, 379,
433, 541, 757, 2161, 7561, 15121, 23, 67, 89,
199, 331, 397, 463, 617, 661, 881, 991, 1321,
2311, 2377, 2971, 3697, 4159, 4621, 8317, 9241,
16633, 18481, 23761, 101, 151, 401, 601, 701,
1051, 1201, 1801, 2801, 3301, 3851, 4201, 4951,
6301, 9901, 11551, 12601, 14851, 15401, 19801,
97, 353, 673, 2017, 3169, 3361, 5281, 7393,
21601, 30241, 53, 79, 131, 157, 313, 521, 547,
859, 911, 937, 1093, 1171, 1249, 1301, 1873,
1951, 2003, 2081, 41, 2731, 2861, 3121, 3433,
3511, 5851, 6007, 6553, 7151, 7723, 8009, 8191,
8581, 8737, 9829, 11701, 13729, 14561, 15601,
16381, 17551, 20021, 20593, 21841, 24571, 25741,
26209, 28081};
static final int aiG[] = {
1, 2, 2, 3, 2, 2, 3, 2, 2, 2, 2, 2, 3, 7,
3, 2, 2, 3, 6, 5, 3, 17, 3, 3, 7, 10, 11,
6, 6, 2, 5, 2, 2, 23, 13, 11, 5, 2, 3, 3,
3, 5, 3, 3, 2, 3, 6, 13, 3, 5, 10, 5, 3,
2, 6, 13, 15, 13, 7, 2, 6, 3, 7, 2, 7, 11,
11, 3, 6, 2, 11, 6, 10, 2, 7, 11, 2, 6, 13,
5, 3, 5, 5, 7, 22, 7, 5, 7, 11, 2, 3, 2,
5, 10, 3, 2, 2, 17, 5, 5, 2, 7, 2, 10, 3,
5, 3, 7, 3, 2, 7, 5, 7, 2, 3, 10, 7, 3, 3,
17, 6, 5, 10, 6, 23, 6, 23, 2, 3, 3, 5, 11,
7, 6, 11, 19};

final int aiInv[] = new int[PWmax];

static final int aiNP[] = {2,3,3,4,4,4,4,5,5,5,6};
static final int aiNQ[] = {5,8,11,18,22,27,36,59,79,89,136};
static final int aiT[] = {2*2*3, 2*2*3*5, 2*2*3*3*5, 2*2*3*3*5*7,
2*2*2*3*3*5*7, 2*2*2*2*3*3*5*7, 2*2*2*2*3*3*3*5*7,
2*2*2*2*3*3*3*5*7*11, 2*2*2*2*3*3*3*5*5*7*11,
2*2*2*2*2*3*3*3*5*5*7*11, 2*2*2*2*2*3*3*3*5*5*7*11*13};

final long biTmp[] = new long[NLen];
final long biExp[] = new long[NLen];
final long biN[] = new long[NLen];
final long biR[] = new long[NLen];
final long biS[] = new long[NLen];
final long biT[] = new long[NLen];

final long biU[] = new long[NLen]; /* Temp */
final long biV[] = new long[NLen]; /* Temp */
final long biW[] = new long[NLen]; /* Temp */

final long aiJS[][] = new long[PWmax][NLen];
final long aiJW[][] = new long[PWmax][NLen];
final long aiJX[][] = new long[PWmax][NLen];
final long aiJ0[][] = new long[PWmax][NLen];
final long aiJ1[][] = new long[PWmax][NLen];
final long aiJ2[][] = new long[PWmax][NLen];
final long aiJ00[][] = new long[PWmax][NLen];
final long aiJ01[][] = new long[PWmax][NLen];

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 TestNbr2[] = 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 long DosALa31_1 = DosALa31 - 1;

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

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;

// Quick lookup tables to eliminate most non-squares without needing
a
// full sqrt calculation, from Cohen 'A Course in Computational
Algebraic
// Number Theory' pp. 40-41. This fails to detect non-squares only
6/715
// of the time.
static final BigInteger SIXTYTHREE = new BigInteger("63"); // bit and
mask
static boolean[] q64 = new boolean[64];
static final BigInteger TABLESIZE = new BigInteger("45045"); // =
11*63*65
static boolean[] q11 = new boolean[11];
static boolean[] q63 = new boolean[63];
static boolean[] q65 = new boolean[65];

public static void main (String args[])
throws java.io.IOException {

char c;
String sInput;

Date d;
long starttime;
long finishtime;
long duration;

// initialize tables with 'true' for all quadratic residues
for(int i=0;i<64;i++)q64[i]=false;for(int i=0;i<=32;i++)q64[i*i
%64]=true;
for(int i=0;i<11;i++)q11[i]=false;for(int i=0;i<=5;i++)q11[i*i
%11]=true;
for(int i=0;i<63;i++)q63[i]=false;for(int i=0;i<=31;i++)q63[i*i
%63]=true;
for(int i=0;i<65;i++)q65[i]=false;for(int i=0;i<=32;i++)q65[i*i
%65]=true;

while (true) {

superfac9 superfac9inst = new superfac9();

StringBuffer sbInput = new StringBuffer("");

while ((c = (char)System.in.read()) != '\n' && c != '\r')
sbInput.append(c);
System.in.read();
sInput = sbInput.toString().trim();

if (sInput.charAt(0) == 'f' || sInput.charAt(0) == 'F') {
s = sInput.substring(1).trim();

s1 = 0;
olds1 = 0;
n = superfac9inst.eval(s);

System.out.println('[' + n.toString() + ']');
}
else {
n = new BigInteger(sInput);
}

d = new Date();
starttime = d.getTime();

superfac9inst.factorize(n);

d = new Date();
finishtime = d.getTime();

duration = (finishtime-starttime)/1000;

System.out.println("duration: " + duration + " seconds");

System.out.println();
}
}

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;
}

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 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;
}

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;
}

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;
}

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;
}

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 boolean isprime(BigInteger n) {
superfac9 superfac9inst = new superfac9();

if (n.bitLength() < 10)
return superfac9inst.brutep(n);
else {
if (n.bitLength() < 2500) {
if (superfac9inst.aprtcle(n) == 0)
return true;
else
return false;
}
else
return superfac9inst.fermat(n);
}
}

boolean brutep(BigInteger n) {
BigInteger a = new BigInteger("2");

System.out.print("brutep: ");

while (a.compareTo(HUNDRED) < 0 && a.multiply(a).compareTo(n) <= 0)
{
if (n.remainder(a).compareTo(ZERO) == 0)
return false;
else
a = a.add(ONE);
}

return true;
}

// Prime checking routine
// Return codes: 0 = Number is prime.
// 1 = Number is composite.
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;

System.out.print("aprtcle: ");

BigIntToBigNbr(N);
GetMontgomeryParms();
j = PK = PL = PM = 0;
for (I=0; I<NumberLength; I++) {
biS[I] = 0;
for (J=0; J<PWmax; J++) {
aiJX[J][I] = 0;
}
}
GetPrimes2Test:
for (i=0; i<LEVELmax; i++) {
biS[0] = 2;
for (I=1; I<NumberLength; I++) {
biS[I] = 0;
}
for (j=0; j<aiNQ[i]; j++) {
Q = aiQ[j];
U = aiT[i]*Q;
do {
U /= Q;
MultBigNbrByLong(biS, Q, biS);
} while (U % Q == 0);

// Exit loop if S^2 > N.

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

if (i == LEVELmax) { /* too big */
System.out.println("Too Big!");
return 1;
}

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)BigNbrModLong(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;}
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<NumberLength; J++) {
aiJ0[I][J] = aiJ1[I][J] = 0;
}
}
if (P>2) {
JacobiSum(1,1,P,PK,PL,PM,Q);
}
else {
if (K != 1) {
JacobiSum(1,1,P,PK,PL,PM,Q);
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = 0;
}
}
if (K != 2) {
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = aiJ0[I][J];
}
}
JacobiSum(2,1,P,PK,PL,PM,Q);
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ0[I][J];
}
}
JS_JW(PK, PL, PM, P);
NormalizeJS(PK, PL, PM, P);
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJ1[I][J] = aiJS[I][J];
}
}
JacobiSum(3 << (K-3),1 << (K-3),P,PK,PL,PM,Q);
for (J=0; J<NumberLength; J++) {
for (I=0; I<PK; I++) {
aiJW[I][J] = 0;
}
for (I=0; I<PM; I++) {
aiJS[I][J] = aiJ0[I][J];
}
}
JS_2(PK, PL, PM, P);
NormalizeJS(PK, PL, PM, P);
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJ2[I][J] = aiJS[I][J];
}
}
}
}
}
for (J=0; J<NumberLength; J++) {
aiJ00[0][J] = aiJ01[0][J] = MontgomeryMultR1[J];
for (I=1; I<PK; I++) {
aiJ00[I][J] = aiJ01[I][J] = 0;
}
}
VK = (int)BigNbrModLong(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<NumberLength; J++) {
aiJS[I][J] = aiJ0[I][J];
}
}
if (X%P == 0) {continue;}
if (IV == 0) {
LongToBigNbr(X, biExp);
}
else {
LongToBigNbr(VK * X / PK, biExp);
if (VK * X / PK == 0) {continue;}
}
JS_E(PK, PL, PM, P);
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = 0;
}
}
InvX = aiInv[X];
for (I=0; I<PK; I++) {
J = I*InvX % PK;
AddBigNbrModN(aiJW[J],aiJS[I],aiJW[J]);
}
NormalizeJW(PK, PL, PM, P);
if (IV == 0) {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ00[I][J];
}
}
}
else {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ01[I][J];
}
}
}
JS_JW(PK, PL, PM, P);
if (IV == 0) {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJ00[I][J] = aiJS[I][J];
}
}
}
else {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJ01[I][J] = aiJS[I][J];
}
}
}
} /* end for X */
} /* end for IV */
}
else {
if (K == 1) {
MultBigNbrByLongModN(MontgomeryMultR1,Q,aiJ00[0]);
for (J=0; J<NumberLength; J++) {
aiJ01[0][J] = MontgomeryMultR1[J];
}
}
else {
if (K == 2) {
if (VK == 1) {
for (J=0; J<NumberLength; J++) {
aiJ01[0][J] = MontgomeryMultR1[J];
}
}
for (J=0; J<NumberLength; J++) {
aiJS[0][J] = aiJ0[0][J];
aiJS[1][J] = aiJ0[1][J];
}
JS_2(PK, PL, PM, P);
if (VK == 3) {
for (J=0; J<NumberLength; J++) {
aiJ01[0][J] = aiJS[0][J];
aiJ01[1][J] = aiJS[1][J];
}
}
MultBigNbrByLongModN(aiJS[0], Q, aiJ00[0]);
MultBigNbrByLongModN(aiJS[1], Q, aiJ00[1]);
}
else {
for (IV=0; IV<=1; IV++) {
for (X=1; X<PK; X+=2) {
for (I=0; I<=PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ1[I][J];
}
}
if (X%8 == 5 || X%8 == 7) {continue;}
if (IV == 0) {
LongToBigNbr(X, biExp);
}
else {
LongToBigNbr(VK * X / PK, biExp);
if (VK * X / PK == 0) {continue;}
}
JS_E(PK, PL, PM, P);
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = 0;
}
}
InvX = aiInv[X];
for (I=0; I<PK; I++) {
J = I*InvX%PK;
AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J]);
}
NormalizeJW(PK, PL, PM, P);
if (IV == 0) {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ00[I][J];
}
}
}
else {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ01[I][J];
}
}
}
NormalizeJS(PK, PL, PM, P);
JS_JW(PK, PL, PM, P);
if (IV == 0) {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJ00[I][J] = aiJS[I][J];
}
}
}
else {
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJ01[I][J] = aiJS[I][J];
}
}
}
} /* end for X */
if (IV==0 || VK%8==1 || VK%8 == 3) {continue;}
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = aiJ2[I][J];
aiJS[I][J] = aiJ01[I][J];
}
}
for (; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = aiJS[I][J] = 0;
}
}
JS_JW(PK, PL, PM, P);
for (I=0; I<PM; I++) {
for (J=0; J<NumberLength; J++) {
aiJ01[I][J] = aiJS[I][J];
}
}
} /* end for IV */
}
}
}
for (I=0; I<PL; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJ00[I][J];
}
}
for (; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = 0;
}
}
DivBigNbrByLong(TestNbr, PK, biExp);
JS_E(PK, PL, PM, P);
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = 0;
}
}
for (I=0; I<PL; I++) {
for (J=0; J<PL; J++) {
MontgomeryMult(aiJS[I], aiJ01[J], biTmp);
AddBigNbrModN(biTmp, aiJW[(I+J)%PK], aiJW[(I+J)%PK]);
}
}
NormalizeJW(PK, PL, PM, P);
MatchingRoot:
do {
H = -1;
W = 0;
for (I=0; I<PL; I++) {
if (BigNbrIsZERO(aiJW[I]) == false) {
if (H == -1 &&
BigNbrAreEqual(aiJW[I], MontgomeryMultR1) ==
true) {
H = I;
}
else {
H = -2;
AddBigNbrModN(aiJW[I], MontgomeryMultR1, biTmp);
if (BigNbrIsZERO(biTmp)) {W++;}
}
}
}
if (H >= 0) {
break MatchingRoot;
}
if (W != P-1) {
return 1; /* Not prime */
}
for (I=0; I<PM; I++) {
AddBigNbrModN(aiJW[I], MontgomeryMultR1, biTmp);
if (BigNbrIsZERO(biTmp) == true) {
break;
}
}
if (I == PM) {
return 1; /* Not prime */
}
for (J=1; J<=P-2; J++) {
AddBigNbrModN(aiJW[I+J*PM], MontgomeryMultR1, biTmp);
if (BigNbrIsZERO(biTmp) == false) {
return 1; /* Not prime */
}
}
H = I+PL;
} while (false);
if (SW == 1 || H%P == 0) {continue;}
if (P != 2) {SW = 1; continue;}
if (K == 1) {
if ((TestNbr[0] & 3) == 1) {
SW = 1;
}
continue;
}

// if (Q^((N-1)/2) mod N != N-1), N is not prime.

MultBigNbrByLongModN(MontgomeryMultR1, Q, biTmp);
for (I = 0; I < NumberLength; I++) {
biR[I] = biTmp[I];
}
I = NumberLength - 1;
Mask = DosALa31;
while ((TestNbr[I] & Mask) == 0) {
Mask /= 2;
if (Mask == 0) {
I--;
Mask = DosALa31;
}
}
do {
Mask /= 2;
if (Mask == 0) {
I--;
Mask = DosALa31;
}
MontgomeryMult(biR,biR,biT);
for (J=0; J<NumberLength; J++) {
biR[J] = biT[J];
}
if ((TestNbr[I] & Mask) != 0) {
MontgomeryMult(biR,biTmp,biT);
for (J=0; J<NumberLength; J++) {
biR[J] = biT[J];
}
}
} while (I>0 || Mask>2);
AddBigNbrModN(biR, MontgomeryMultR1, biTmp);
if (BigNbrIsZERO(biTmp) == false) {
return 1; /* Not prime */
}
SW = 1;
} /* end for j */
if (SW == 0) {
TestedQs = TestingQs+1;
if (TestingQs < aiNQ[LEVELnow]-1) {
TestingQs++;
Q = aiQ[TestingQs];
U = T*Q;
do {
MultBigNbrByLong(biS, Q, biS);
U /= Q;
} while (U%Q == 0);
continue; /* Retry */
}
LEVELnow++;

if (LEVELnow==LEVELmax) {
System.out.println("Cannot Tell!");
return 1; /* Cannot tell */
}

T = aiT[LEVELnow];
NP = aiNP[LEVELnow];
biS[0] = 2;
for (J=1; J<NumberLength; J++) {
biS[J] = 0;
}
for (J=0; J<=aiNQ[LEVELnow]; J++) {
Q = aiQ[J];
U = T*Q;
do {
MultBigNbrByLong(biS, Q, biS);
U /= Q;
} while (U%Q == 0);
if (CompareSquare(biS, TestNbr) > 0) {
TestingQs = J;
continue MainStart; /* Retry from the beginning */
}
} /* end for J */

System.out.println("Program error!");
return 1;

} /* end if */
break;
} while (true); /* end do */
} /* end for i */

// Final Test

LengthN = NumberLength;
for (I=0; I<NumberLength; I++) {
biN[I] = TestNbr[I];
TestNbr[I] = biS[I];
biR[I] = 0;
}
while (true) {
if (TestNbr[NumberLength-1] != 0) {break;}
NumberLength--;
}
dN = (double)TestNbr[NumberLength-1];
if (NumberLength > 1) {
dN += (double)TestNbr[NumberLength-2]/dDosALa32;
}
if (NumberLength > 2) {
dN += (double)TestNbr[NumberLength-3]/dDosALa64;
}
LengthS = NumberLength;
dS = dN;
MontgomeryMultR1[0] = 1;
for (I=1; I<NumberLength; I++) {
MontgomeryMultR1[I] = 0;
}

biR[0] = 1;
BigNbrModN(biN, LengthN, biT); /* Compute N mod S */
for (J=1; J<=T; J++) {
MultBigNbrModN(biR, biT, biTmp);
for (i=NumberLength-1; i>0; i--) {
if (biTmp[i] != 0)
{
break;
}
}
if (i==0 && biTmp[0] != 1) {
return 0; /* Number is prime */
}
while (true) {
if (biTmp[NumberLength-1] != 0) {break;}
NumberLength--;
}
for (I=0; I<NumberLength; I++) {
TestNbr[I] = biTmp[I];
}
dN = (double)TestNbr[NumberLength-1];
if (NumberLength > 1) {
dN += (double)TestNbr[NumberLength-2]/dDosALa32;
}
if (NumberLength > 2) {
dN += (double)TestNbr[NumberLength-3]/dDosALa64;
}
for (i=NumberLength-1; i>0; i--) {
if (TestNbr[i] != biTmp[i]) {
break;
}
}
if (TestNbr[i] > biTmp[i]) {
BigNbrModN(biN, LengthN, biTmp); /* Compute N mod R
*/
if (BigNbrIsZERO(biTmp) == true) { /* If N is
multiple of R.. */
return 1; /* Number is
composite */
}
}
dN = dS;
NumberLength = LengthS;
for (I=0; I<NumberLength; I++) {
biR[I] = TestNbr[I];
TestNbr[I] = biS[I];
}
} /* End for J */
return 0; /* Number is prime */
} while (true);
}

boolean fermat(BigInteger n) {
boolean prime = false;
BigInteger T = new BigInteger("1");
BigInteger a = new BigInteger("2");

System.out.print("fermat: ");

// workaround - apparent java bug in modPow - JGW
if (n.compareTo(TWO) < 0)
return false;

if (n.remainder(TWO).compareTo(ZERO) == 0) {
return false;
}
// end workaround

T = n;

while (a.compareTo(HUNDRED) < 0 && T.compareTo(n) == 0) {
T = n.gcd(a.modPow(n, n).subtract(a));
a = a.add(ONE);
}

if (T.compareTo(n) < 0)
prime = false;
else
prime = true;

return prime;
}

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

if (isprime(n)) {
prime = true;
System.out.println(n);
return prime;
}

System.out.print('\r' + "wanless... ");
factor = factorizewanless(n, THREE);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "brute... ");
factor = factorizebrute(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

/* else {

System.out.print('\r' + "rho... ");
factor = factorizerho(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "brent... ");
factor = factorizebrent(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "pollard... ");
factor = factorizepollard(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "williams... ");
factor = factorizewilliams(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "fermat... ");
factor = factorizefermat(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "shanks... ");
factor = factorizeshanks(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}
*/
else {

System.out.print('\r' + "siqs... ");
factor = factorizesiqs(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

System.out.print('\r' + "ecm... ");
factor = factorizeecm(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

/* else {

System.out.print('\r' + "wanless... ");
factor = factorizewanless(n, TENTHOUSAND);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

*/ else {

System.out.print('\r' + "wanless ext... ");
factor = factorizewe(n);

if (factor.compareTo(ONE) > 0 && factor.compareTo(n) < 0) {
System.out.println();
factorize(factor);
factorize(n.divide(factor));
}

else {

prime = false;
System.out.println("composite" + '\t' + n);

/* }
}
}
}
}
}
}
*/ }
}
}
}
}

return prime;

}

BigInteger factorizewanless(BigInteger n, BigInteger exTENt) {
BigInteger T = new BigInteger("1");
BigInteger a = new BigInteger("2");

// workaround - apparent java bug in modPow - JGW
if (n.compareTo(TWO) < 0)
return ONE;

if (n.remainder(TWO).compareTo(ZERO) == 0)
return TWO;
// end workaround

while (a.multiply(a).compareTo(n) <= 0 && (T.compareTo(ONE) == 0 ||
T.compareTo(n) == 0) && a.compareTo(exTENt)<0) {
T = n.gcd(a.modPow(n, n).subtract(a));
a = a.add(TWO);
}

return T;
}

BigInteger factorizebrute(BigInteger n) {
BigInteger a = new BigInteger("2");

while (a.compareTo(TENTHOUSAND) < 0 && a.multiply(a).compareTo(n) <=
0) {
if (n.remainder(a).compareTo(ZERO) == 0)
return a;
else
a = a.add(ONE);
}

return n;
}

BigInteger factorizerho(BigInteger n) {
BigInteger loop = new BigInteger("1");
BigInteger x = new BigInteger("5");
BigInteger y = new BigInteger("2");

while (n.gcd(x.subtract(y)).compareTo(ONE) == 0 &&
loop.compareTo(TENTHOUSAND) < 0) {
x = x.multiply(x).add(ONE).mod(n);
x = x.multiply(x).add(ONE).mod(n);
y = y.multiply(y).add(ONE).mod(n);
loop = loop.add(ONE);
}

return n.gcd(x.subtract(y));

}

BigInteger factorizebrent(BigInteger n) {
BigInteger k = new BigInteger("1");
BigInteger r = new BigInteger("1");
BigInteger i = new BigInteger("1");
BigInteger m = new BigInteger("1");
BigInteger iter = new BigInteger("1");
BigInteger z = new BigInteger("1");
BigInteger x = new BigInteger("1");
BigInteger y = new BigInteger("1");
BigInteger q = new BigInteger("1");
BigInteger ys = new BigInteger("1");

m=TEN;
r=ONE;
iter=ZERO;
z=ZERO;
y=z;
q=ONE;

do {
x=y;
for (i=ONE;i.compareTo(r)<=0;i=i.add(ONE))
y=((y.multiply(y)).add(THREE)).mod(n);
k=ZERO;
do {
iter=iter.add(ONE);
// System.out.print("iter=" + iter.toString() + '\r');
ys=y;
for (i=ONE;i.compareTo(mr_min(m,r.subtract(k)))<=0;i=i.add(ONE)) {
y=((y.multiply(y)).add(THREE)).mod(n);
q=((y.subtract(x)).multiply(q)).mod(n);
}
z=n.gcd(q);
k=k.add(m);
} while (k.compareTo(r)<0 && z.compareTo(ONE)==0);
r=r.multiply(TWO);
} while (z.compareTo(ONE)==0 && iter.compareTo(TENTHOUSAND)<0);

if (z.compareTo(n)==0) do {
ys=((ys.multiply(ys)).add(THREE)).mod(n);
z=n.gcd(ys.subtract(x));
} while (z.compareTo(ONE)==0);

return z;
}

BigInteger mr_min(BigInteger a, BigInteger b) {
if (a.compareTo(b) < 0)
return a;
else
return b;
}

BigInteger factorizepollard(BigInteger numb) {
BigInteger bound = new BigInteger("0");
BigInteger i = new BigInteger("2");
BigInteger thegcd = new BigInteger("0");
BigInteger base = new BigInteger("2");
BigInteger count = new BigInteger("1");

count = TENTHOUSAND;

bound = root3(numb);

for (i = TWO; i.compareTo(bound) < 0 && i.compareTo(count) < 0; i =
i.add(ONE)) {
base = i.modPow(base, numb);
if ((thegcd = numb.gcd(base.subtract(ONE))).compareTo(ONE) != 0)
break;
}

return thegcd;
}

BigInteger root3(BigInteger in) {
BigInteger out = new BigInteger("0");

BigInteger factor1 = new BigInteger("0");
BigInteger factor2 = new BigInteger("0");
BigInteger swapfactor = new BigInteger("0");

factor1 = ONE;
factor2 = in;

do {
factor1 = (factor1.add(factor2)).divide(TWO);
factor2 = (in.divide(factor1)).divide(factor1);

if (factor1.compareTo(factor2) > 0) {
swapfactor = factor1;
factor1 = factor2;
factor2 = swapfactor;
}
} while ((factor2.subtract(factor1)).compareTo(ONE) > 0);

out = factor1;

return out;
}

BigInteger factorizewilliams(BigInteger m) {
BigInteger f = new BigInteger("1");
BigInteger n = new BigInteger("1");
BigInteger P = new BigInteger("0");
BigInteger v = new BigInteger("0");
BigInteger count = new BigInteger("1");

for (n = ONE; n.compareTo(FIVE) <= 0; n = n.add(ONE)) {
P = n.add(THREE);

count = ONE;
v = P;

while (count.compareTo(TENTHOUSAND) <= 0) {
if ((f = m.gcd(v.subtract(TWO))).compareTo(ONE) > 0)
break;

for (int i = 1; i <= 10; i++) {
v = nextV(ONE, P, count, m);
P = v;
count = count.add(ONE);
}
}

if (f.compareTo(ONE) > 0)
break;
}

return f;
}

BigInteger nextV(BigInteger n, BigInteger h, BigInteger j, BigInteger
p) {
BigInteger x = new BigInteger("0");
BigInteger m = new BigInteger("0");
BigInteger v = new BigInteger("0");
BigInteger w = new BigInteger("0");
BigInteger tt = new BigInteger("0");
int t = 0;

m = n;
v = h;
w = ((h.multiply(h)).subtract(TWO.multiply(m))).mod(p);

t = bitLength(j);

for (int k = 0; k < t; k++) {
x = ((v.multiply(w)).subtract(h.multiply(m))).mod(p);

tt = TWO.multiply(m);

v = ((v.multiply(v)).subtract(tt)).mod(p);
w = ((w.multiply(w)).subtract(tt.multiply(n))).mod(p);
m = (m.multiply(m)).mod(p);

if (getBit(j, k) == 0)
w = x;
else {
v = x;
m = (n.multiply(m)).mod(p);
}
}

return v;
}

int bitLength(BigInteger j) {
int result = 0;

result = j.bitLength();

return result;
}

int getBit(BigInteger j, int k) {
int result = 0;

if (j.testBit(k) == false)
result = 0;
else
result = 1;

return result;
}

// Original Fermat method, with fast square check
BigInteger factorizefermat(BigInteger n) {
BigInteger a, bSq;
long iterations = 1L;

a = sqrt(n);
if(n.mod(a).compareTo(ZERO) == 0) return a;

while(iterations<10000000) {
bSq = a.multiply(a).subtract(n);
if(fastsquareQ(bSq))
return a.subtract(sqrt(bSq)); // bSq is a square,
factorization found
a = a.add(ONE);
iterations++;
}

return n;
}

// fast table driven test whether n is a square of an integer or
not
boolean fastsquareQ(BigInteger n) {
BigInteger root; int j;

// Try to quickly eliminate most non-squares to avoid costly
sqrt
if(q64[(n.and(SIXTYTHREE)).intValue()]) { // Might be a square
j= n.mod(TABLESIZE).intValue(); // ONE mod 63*65*11 then work
with ints
if(q63[j%63] && q65[j%65] && q11[j%11]) {
// Still might be a non-square, now have to do a full
sized sqrt
root = sqrt(n);
return root.multiply(root).compareTo(n) == 0;
};
};
return false;
}

BigInteger sqrt(BigInteger in) {
BigInteger out = new BigInteger("0");

BigInteger factor1 = new BigInteger("0");
BigInteger factor2 = new BigInteger("0");
BigInteger swapfactor = new BigInteger("0");

factor1 = ONE;
factor2 = in;

do {
factor1 = (factor1.add(factor2)).divide(TWO);
factor2 = in.divide(factor1);

if (factor1.compareTo(factor2) > 0) {
swapfactor = factor1;
factor1 = factor2;
factor2 = swapfactor;
}
} while ((factor2.subtract(factor1)).compareTo(ONE) > 0);

out = factor1;

return out;
}

BigInteger factorizeshanks(BigInteger n) {
BigInteger a = new BigInteger("0");
BigInteger f = new BigInteger("0");
BigInteger h1 = new BigInteger("0");
BigInteger h2 = new BigInteger("0");
BigInteger k = new BigInteger("0");
BigInteger p = new BigInteger("0");
BigInteger pp = new BigInteger("0");
BigInteger q = new BigInteger("0");
BigInteger qq = new BigInteger("0");
BigInteger qqq = new BigInteger("0");
BigInteger r = new BigInteger("0");
BigInteger te = new BigInteger("0");
BigInteger i = new BigInteger("0");
BigInteger count = new BigInteger("0");

k = sqrt(n);

if (fastsquareQ(n)) return k;

a=k; h1=k; h2=ONE; pp=ZERO; qq=ONE; qqq=n; r=ZERO;

for (count=ONE;count.compareTo(TENTHOUSAND)<0;count=count.add(ONE))
{
p=k.subtract(r);
q=qqq.add(a.multiply(pp.subtract(p)));
a=(p.add(k)).divide(q);
r=(p.add(k)).remainder(q);
te=(a.multiply(h1)).add(h2);
h2=h1;
h1=te;
pp=p;
qqq=qq;
qq=q;
te = sqrt(q);
i=i.add(ONE);
if ((i.remainder(TWO).compareTo(ZERO))!=0 || !fastsquareQ(q))
continue;
te=h2.subtract(te);
f=n.gcd(te);
if (f.compareTo(ONE) > 0 && f.compareTo(n) < 0)
return f;
}

return f;
}

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

// added JGW
if (NbrToFactor.bitLength()<100 || NbrToFactor.bitLength()>266)
return NbrToFactor;

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

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

top=-10.0e0;
nk=0;
bk=0;
prime[0]=1;
prime[1]=2;
currentPrime = 3;

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

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

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

// Compute the leading coefficient in biS.

LongToBigNbr(afact[0], biS);
for (index = 1; index < s; index++) {
MultBigNbrByLong(biS, afact[index], biS);
}
for (index = 0; index < s; index++) {
D = 1;
E = afact[index];
for (index2 = 0; index2 < s; index2++) {
if (index != index2) {
D = D * afact[index2] % E;
}
}
amodq[index] = (int)D;
Q = modsqrt[aindex[index]] * modInv(D, E) % E;
if (Q > E/2) {
Q = E - Q;
}
DivBigNbrByLong(biS, E, biR);
MultBigNbrByLong(biR, Q, aiJS[index]);
}
for (index = 0; index < NumberLength; index++) {
biN[index] = aiJS[0][index];
}
for (index2 = 1; index2 < s; index2++) {
AddBigNbr(biN, aiJS[index2], biN);
}
for (index = 1; index < NbrPrimes; index++) {
D = 1;
E = prime[index];
for (index2 = 0; index2 < s; index2++) {
D = D * afact[index2] % E;
}
ainv[index] = modInv(D, E);
for (index2 = 0; index2 < s; index2++) {
Bainv2[index2][index] = (int)
(RemDivBigNbrByLong(aiJS[index2], E) * 2 *
ainv[index] % E);
}
D = RemDivBigNbrByLong(biN, E);
soln1[index] = (int)((SieveLimit + ainv[index]*(E +
modsqrt[index] - D)) % E);
soln2[index] = (int)((SieveLimit + ainv[index]*(E + E -
modsqrt[index] - D)) % E);
}
do {

/***************/
/* Sieve stage */
/***************/
D = PolynomialIndex;
index2 = 0;
while ((D & 1) == 0) {
D /= 2;
index2++;
}
if (polyadd = ((D & 2) != 0)) {
AddBigNbr(biN, aiJS[index2], biN);
AddBigNbr(biN, aiJS[index2], biN);
}
else {
SubtractBigNbr(biN, aiJS[index2], biN);
SubtractBigNbr(biN, aiJS[index2], biN);
}
int [] rowBainv2 = Bainv2[index2];
// Compute solutions for divisors of the leading coefficients
for (index2 = 0; index2 < s; index2++) {
index = aindex[index2];
E = prime[index];
D = RemDivBigNbrByLong(TestNbr, E*E);
Q = RemDivBigNbrByLong(biN, E*E);
soln1[index] = (int)(((D - Q*Q)/E * modInv(amodq[index2]*Q%E,
E)%E+E+SieveLimit)%E);
soln2[index] = 2*SieveLimit+1;
}
X1 = 2 * SieveLimit;
F1 = polyadd? 2 - rowBainv2[1]: rowBainv2[1];
if ((soln1[1] += F1)%2 == 0) {
SieveArray[0] = (byte)(logar[1] - G);
SieveArray[1] = (byte)(-G);
}
else {
SieveArray[0] = (byte)(-G);
SieveArray[1] = (byte)(logar[1] - G);
}
F2 = 2;
index = 2;
while (true) {
F = (int)prime[index];
F3 = F2 * F;
if (X1+1 < F3) {
F3 = X1+1;
}
F4 = F2;
while (F4*2 <= F3) {
System.arraycopy(SieveArray, 0, SieveArray, F4, F4);
F4 *= 2;
}
System.arraycopy(SieveArray, 0, SieveArray, F4, F3-F4);
if (F3 == X1+1) {
break;
}
Q1 = logar[index];
F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
for (index2=(soln1[index] += F1) % F; index2 < F3; index2 +=
F) {
SieveArray[index2] += Q1;
}
for (index2=(soln2[index] += F1) % F; index2 < F3; index2 +=
F) {
SieveArray[index2] += Q1;
}
index++;
F2 *= F;
}

for (; index<smallPrimeUpperLimit; index++) {
F = (int)prime[index];
F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
soln1[index] += F1;
soln2[index] += F1;
}

int primesmall = (int)prime[index-1];
int soln1small = soln1[index-1];
int soln2small = soln2[index-1];
for (; index<firstLimit; index++) {
F = (int)prime[index];
F1 = polyadd? F - rowBainv2[index]: rowBainv2[index];
F2 = F+F; F3=F2+F; F4=F3+F;
S1 = (soln1[index] += F1) % F;
S2 = (soln2[index] += F1) % F;
if ((G0 = S2 - S1) < 0) {
S1 = S2;
G0 = -G0;
}
G1 = G0+F; G2 = G1+F; G3 = G2+F;
Q1 = logar[index];
index2=X1/F4*F4+S1;
do {
SieveArray[index2] += Q1;
SieveArray[index2+F] += Q1;
SieveArray[index2+F2] += Q1;
SieveArray[index2+F3] += Q1;
SieveArray[index2+G0] += Q1;
SieveArray[index2+G1] += Q1;
SieveArray[index2+G2] += Q1;
SieveArray[index2+G3] += Q1;
} while ((index2 -= F4) >= 0);
}
for (; index<secondLimit; index++) {
F = (int)prime[index];
F1 = (polyadd? F - rowBainv2[index]: rowBainv2[index]);
X2 = X1 - 4*F;
F2 = F+F; F3=F2+F; F4=F2+F2;
Q1 = logar[index];
for (index2=(soln1[index] += F1)%F; index2 <= X2; index2 +=
F4) {
SieveArray[index2] += Q1;
SieveArray[index2+F] += Q1;
SieveArray[index2+F2] += Q1;
SieveArray[index2+F3] += Q1;
}
for (; index2 <= X1; index2 += F) {
SieveArray[index2] += Q1;
}
for (index2=(soln2[index] += F1)%F; index2 <= X2; index2 +=
F4) {
SieveArray[index2] += Q1;
SieveArray[index2+F] += Q1;
SieveArray[index2+F2] += Q1;
SieveArray[index2+F3] += Q1;
}
for (; index2 <= X1; index2 += F) {
SieveArray[index2] += Q1;
}
}
for (; index<thirdLimit; index++) {
F = (int)prime[index];
Q1 = logar[index];
F1 = (polyadd? F - rowBainv2[index]: rowBainv2[index]);
for (index2=(soln1[index] += F1)%F; index2 <= X1; index2 += F)
{
SieveArray[index2] += Q1;
}
for (index2=(soln2[index] += F1)%F; index2 <= X1; index2 += F)
{
SieveArray[index2] += Q1;
}
}
if (polyadd) {
for (; index<NbrPrimes2; index++) {
Q1 = logar[index];
F1 = (F = (int)prime[index]) - rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F1 = (F = (int)prime[++index]) - rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F1 = (F = (int)prime[++index]) - rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F1 = (F = (int)prime[++index]) - rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
}
for (; index<NbrPrimes; index++) {
Q1 = logar[index];
F1 = (F = (int)prime[index]) - rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
}
}
else {
for (; index<NbrPrimes2; index++) {
Q1 = logar[index];
F = (int)prime[index];
F1 = rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F = (int)prime[++index];
F1 = rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F = (int)prime[++index];
F1 = rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
F = (int)prime[++index];
F1 = rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
}
for (; index<NbrPrimes; index++) {
Q1 = logar[index];
F = (int)prime[index];
F1 = rowBainv2[index];
if ((F2 = (soln1[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
if ((F2 = (soln2[index] += F1) % F) < X1) {
SieveArray[F2] += Q1;
}
}
}

/************************/
/* Trial division stage */
/************************/
index2 = 2*SieveLimit+1;
do {
if (SieveArray[--index2] >= 0 && SieveArray[index2]<64) {
if (SieveArray[index2] < logarsmall) {
if ((index2-soln1small)%primesmall != 0) {
if ((index2-soln2small)%primesmall != 0) {continue;}
}
}

MultBigNbrByLong(biS, index2 - SieveLimit, biT);
AddBigNbr(biT, biN, biT);
MultBigNbr(biT, biT, biR);
SubtractBigNbr(biR, TestNbr, biR); // Number to factor: (Ax
+B)^2-N
for (i=0; i<NumberLength; i++) {
biT[i] = biR[i];
}
/* factor biR */

if (kk>1) {
while (RemDivBigNbrByLong(biR, kk*kk)==0) {
DivBigNbrByLong(biR, kk*kk, biR);
}
}
F = NumberLength; /* Back up NumberLength */
for (index=0; index<s; index++) {
DivBigNbrByLong(biR, afact[index], biR);
if (biR[NumberLength-1] == 0 && biR[NumberLength-2] <
DosALa31 ||
biR[NumberLength-1] == MaxUInt && biR[NumberLength-2]
>= DosALa31) {
NumberLength--;
}
}
switch (NumberLength) {
case 6:
biR5 = biR[5];
case 5:
biR4 = biR[4];
case 4:
biR3 = biR[3];
case 3:
biR2 = biR[2];
case 1:
case 2:
biR1 = biR[1];
biR0 = biR[0];
}
nbrColumns = 0;
if (NumberLength <= 2) {
Divid = (biR1 << 32) | biR0;
for (index=1; index<NbrPrimes; index++) {
Divisor = prime[index];
while (Divid % Divisor == 0) {
Divid /= Divisor;
}
}
}
else {
Divid = 0;
if (biR[NumberLength-1] < DosALa31) { /* Positive */
for (index=1; index<NbrPrimes; index++) {
while (true) {
Divisor = prime[index];
Rem = 0;
switch (NumberLength) {
case 6:
Rem = (biR5 + (Rem << 32)) % Divisor;
case 5:
Rem = (biR4 + (Rem << 32)) % Divisor;
case 4:
Rem = (biR3 + (Rem << 32)) % Divisor;
case 3:
Rem = (biR2 + (Rem << 32)) % Divisor;
Rem = (biR1 + (Rem << 32)) % Divisor;
Rem = (biR0 + (Rem << 32)) % Divisor;
}
if (Rem != 0) {
break;
}
switch (NumberLength) {
case 6:
Divid = biR5 + (Rem << 32);
Rem = Divid % Divisor;
biR5 = Divid / Divisor;
case 5:
Divid = biR4 + (Rem << 32);
Rem = Divid % Divisor;
biR4 = Divid / Divisor;
case 4:
Divid = biR3 + (Rem << 32);
Rem = Divid % Divisor;
biR3 = Divid / Divisor;
case 3:
Divid = biR2 + (Rem << 32);
Rem = Divid % Divisor;
biR2 = Divid / Divisor;
Divid = biR1 + (Rem << 32);
biR1 = Divid / Divisor;
biR0 = (biR0 + ((Divid % Divisor) << 32)) /
Divisor;
}
switch (NumberLength) {
case 6:
cond = (biR5 == 0 && biR4 < DosALa31);
break;
case 5:
cond = (biR4 == 0 && biR3 < DosALa31);
break;
case 4:
cond = (biR3 == 0 && biR2 < DosALa31);
break;
case 3:
cond = (biR2 == 0 && biR1 < DosALa31);
break;
}
if (cond) {
NumberLength--;
if (NumberLength == 2) {
Divid = (biR1 << 32) | biR0;
for (; index<NbrPrimes; index++) {
Divisor = prime[index];
while (Divid % Divisor == 0) {
Divid /= Divisor;
}
}
break;
}
}
}
}
}
else { /* Negative */
for (index=1; index<NbrPrimes; index++) {
while (true) {
Divisor = prime[index];
Rem = Divisor - 1;
switch (NumberLength) {
case 6:
Rem = (biR5 + (Rem << 32)) % Divisor;
case 5:
Rem = (biR4 + (Rem << 32)) % Divisor;
case 4:
Rem = (biR3 + (Rem << 32)) % Divisor;
case 3:
Rem = (biR2 + (Rem << 32)) % Divisor;
Rem = (biR1 + (Rem << 32)) % Divisor;
Rem = (biR0 + (Rem << 32)) % Divisor;
}
if (Rem != 0) {
break;
}
Rem = Divisor - 1;
switch (NumberLength) {
case 6:
Divid = biR5 + (Rem << 32);
Rem = Divid % Divisor;
biR5 = Divid / Divisor;
case 5:
Divid = biR4 + (Rem << 32);
Rem = Divid % Divisor;
biR4 = Divid / Divisor;
case 4:
Divid = biR3 + (Rem << 32);
Rem = Divid % Divisor;
biR3 = Divid / Divisor;
case 3:
Divid = biR2 + (Rem << 32);
Rem = Divid % Divisor;
biR2 = Divid / Divisor;
Divid = biR1 + (Rem << 32);
biR1 = Divid / Divisor;
biR0 = (biR0 + ((Divid % Divisor) << 32)) /
Divisor;
}
switch (NumberLength) {
case 6:
cond = (biR5 == MaxUInt && biR4 >= DosALa31);
break;
case 5:
cond = (biR4 == MaxUInt && biR3 >= DosALa31);
break;
case 4:
cond = (biR3 == MaxUInt && biR2 >= DosALa31);
break;
case 3:
cond = (biR2 == MaxUInt && biR1 >= DosALa31);
break;
}
if (cond) {
NumberLength--;
if (NumberLength == 2) {
Divid = (biR1 << 32) | biR0;
for (; index<NbrPrimes; index++) {
Divisor = prime[index];
while (Divid % Divisor == 0) {
Divid /= Divisor;
}
}
break;
}
}
}
}
}
} /* end if */
F2 = NumberLength; NumberLength = F;
if (F2==2 && (Divid==1 || Divid == -1)) { // Smooth relation
found.
if (Divid == -1) {
rowMatrixB[nbrColumns++] = 0; // Insert -1 as a
factor.
}
/* factor biT (backup) storing in BiR the square part */
LongToBigNbr(1, biR);
while (RemDivBigNbrByLong(biT, 3*3)==0) {
DivBigNbrByLong(biT, 3*3, biT);
if (RemDivBigNbrByLong(TestNbr, 3) == 0) {
DivBigNbrByLong(biU, 3, biU);
}
else {
MultBigNbrByLong(biR, 3, biR);
}
}
for (index=5; index<100; index+=4) {
while (RemDivBigNbrByLong(biT, index*index)==0) {
DivBigNbrByLong(biT, index*index, biT);
if (RemDivBigNbrByLong(TestNbr, index) == 0) {
DivBigNbrByLong(biU, index, biU);
}
else {
MultBigNbrByLong(biR, index, biR);
}
}
index += 2;
while (RemDivBigNbrByLong(biR, index*index)==0) {
DivBigNbrByLong(biR, index*index, biR);
if (RemDivBigNbrByLong(TestNbr, index) == 0) {
DivBigNbrByLong(biU, index, biU);
}
else {
MultBigNbrByLong(biR, index, biR);
}
}
}
MultBigNbrByLong(biS, index2 - SieveLimit, biU);
AddBigNbr(biU, biN, biU);
for (index=1; index<NbrPrimes; index++) {
expParity = 0;
D = prime[index];
while (RemDivBigNbrByLong(biT, D)==0) {
DivBigNbrByLong(biT, D, biT);
expParity = 1 - expParity;
if (expParity == 0) {
if (RemDivBigNbrByLong(TestNbr, D) == 0) {
DivBigNbrByLong(biU, D, biU);
}
else {
MultBigNbrByLong(biR, D, biR);
}
}
}
if (expParity != 0) {
rowMatrixB[nbrColumns++] = index;
expParity = 0;
}
}

if (InsertNewRelation(biR, biT, biU, nbrColumns, matrixB,
rowMatrixB, pp,
vectLeftHandSide)) {

pp++;

if (pp == matrixB.length) {
i = EraseSingletons(matrixB, vectLeftHandSide,
vectExpParity);
if (i!=0) {

pp -= i;
}
else {

if (LinearAlgebraPhase(NbrPrimes, matrixB, prime,
biT, biR, biU, vectExpParity,
vectLeftHandSide)) {
return BigNbrToBigInt(biT); /* Factor found
*/
}
else {

pp -= 50; /* Factor not found */
}
}
}
}
continue; /* Continue sieving */
}
else {
if (F2==2 && (Divid<64*FactorBase && Divid >
-64*FactorBase)) {
// Partial relation
found.

// Check if there is already another relation with the
same
// factor outside the prime base.
// Calculate hash index
i = matrixPartialHashIndex[(int)(Divid & 0x7FE)/2];
prev = -1;
while (i >= 0) {
if ((int)Divid == matrixPartial[i][0]) {
// Match of partials.
for (index=0; index<NumberLength; index++) {
biV[index] = matrixPartial[i][index+2] & MaxUInt;
}
MultBigNbr(biV, biV, biT);
SubtractBigNbr(biT, TestNbr, biT);

/* factor biT */

D = (long)matrixPartial[i][0];
E = RemDivBigNbrByLong(TestNbr, D);
t1 = D;
t2 = E;
while (t1 != 0) {
t3 = t2%t1;
t2 = t1;
t1 = t3;
} // t2 = GCD(D, E)
LongToBigNbr(D, biR);
if (D<0) {
AddBigNbr(biR, TestNbr, biR);
}
if (t2 < 0) {t2=-t2;}
if (t2 != 1) {
DivBigNbrByLong(TestNbr, t2<0?-t2:t2, biU);
AddBigNbrModN(biR, biU, biR);
}
if (kk>1) {
while (RemDivBigNbrByLong(biT, kk*kk)==0) {
DivBigNbrByLong(biT, kk*kk, biT);
NumberLength--;
MultBigNbrByLongModN(biR, kk, biR);
NumberLength++;
if (RemDivBigNbrByLong(TestNbr, kk) == 0) {
DivBigNbrByLong(TestNbr, kk, biU);
AddBigNbrModN(biR, biU, biR);
}
}
}
nbrColumns = 0;
for (index=1; index<NbrPrimes; index++) {
D = prime[index];
expParity = 0;
while (RemDivBigNbrByLong(biT, D)==0) {
expParity = 1 - expParity;
DivBigNbrByLong(biT, D, biT);
if (expParity == 0) {
NumberLength--;
MultBigNbrByLongModN(biR, D, biR);
NumberLength++;
if (RemDivBigNbrByLong(TestNbr, D) == 0) {
DivBigNbrByLong(TestNbr, D, biU);
AddBigNbrModN(biR, biU, biR);
}
}
}
if (expParity != 0) {
rowMatrixB[nbrColumns++] = index;
}
}

MultBigNbrByLong(biS, index2 - SieveLimit, biT);
AddBigNbr(biT, biN, biW);
MultBigNbr(biW, biW, biT);
SubtractBigNbr(biT, TestNbr, biT);

/* factor biT */

if (kk>1) {
while (RemDivBigNbrByLong(biT, kk*kk)==0) {
DivBigNbrByLong(biT, kk*kk, biT);
NumberLength--;
MultBigNbrByLongModN(biR, kk, biR);
NumberLength++;
if (RemDivBigNbrByLong(TestNbr, kk) == 0) {
DivBigNbrByLong(TestNbr, kk, biU);
AddBigNbrModN(biR, biU, biR);
}
}
}
for (index=1; index<NbrPrimes; index++) {
expParity = 0;
D = prime[index];
while (RemDivBigNbrByLong(biT, D)==0) {
expParity = 1 - expParity;
DivBigNbrByLong(biT, D, biT);
if (expParity == 0) {
NumberLength--;
MultBigNbrByLongModN(biR, D, biR);
NumberLength++;
if (RemDivBigNbrByLong(TestNbr, D) == 0) {
DivBigNbrByLong(TestNbr, D, biU);
AddBigNbrModN(biR, biU, biR);
}
}
}
if (expParity != 0) {
// Check if the index is already in the row.
for (j=0; j<nbrColumns; j++) {
if (index <= rowMatrixB[j]) {
break;
}
}
if (j<nbrColumns && index == rowMatrixB[j]) {
// Index already in row.
D = prime[rowMatrixB[j]];
NumberLength--;
MultBigNbrByLongModN(biR, D, biR);
NumberLength++;
if (RemDivBigNbrByLong(TestNbr, D) == 0) {
DivBigNbrByLong(TestNbr, D, biU);
AddBigNbrModN(biR, biU, biR);
}
// Delete entry from row.
for (k=j+1; k<nbrColumns; k++) {
rowMatrixB[k-1] = rowMatrixB[k];
}
nbrColumns--;
}
else { // Index not in row.
// Insert entry in row.
for (k=nbrColumns; k>j; k--) {
rowMatrixB[k] = rowMatrixB[k-1];
}
rowMatrixB[j] = index;
nbrColumns++;
}
}
}
if ((biV[NumberLength-1] & 0x80000000L) != 0) {
AddBigNbr(biV, TestNbr, biV);
}
if ((biW[NumberLength-1] & 0x80000000L) != 0) {
AddBigNbr(biW, TestNbr, biW);
}
NumberLength--;
MultBigNbrModN(biV, biW, biU);
NumberLength++;
if (InsertNewRelation(biR, biT, biU, nbrColumns,
matrixB, rowMatrixB, pp, vectLeftHandSide))
{

pp++;

if (pp == matrixB.length) {
i = EraseSingletons(matrixB, vectLeftHandSide,
vectExpParity);
if (i!=0) {

pp -= i;
}
else {

if (LinearAlgebraPhase(NbrPrimes, matrixB,
prime,
biT, biR, biU, vectExpParity,
vectLeftHandSide)) {
return BigNbrToBigInt(biT); /* Factor
found */
}
else {

pp -= 50; /* Factor not found
*/
}
}
}
break;
}
i = -2; // Do not execute next if.
break;
}
prev = i;
i = matrixPartial[i][1]; // Get next index for same
hash.
} /* end while */
if (i==-1 && nbrPartials<matrixPartial.length) { // No
match.
// Add partial to table of partials.
if (prev >= 0) {
matrixPartial[prev][1] = nbrPartials;
}
else {
matrixPartialHashIndex[(int)(Divid & 0x7FE)/2] =
nbrPartials;
}
matrixPartial[nbrPartials][0] = (int)Divid;
matrixPartial[nbrPartials][1] = -1; // Indicate last
index
// with this
hash.
MultBigNbrByLong(biS, index2 - SieveLimit, biT);
AddBigNbr(biT, biN, biT); // biT = Ax+B
for (index=0; index<NumberLength; index++) {
matrixPartial[nbrPartials][index+2] =
(int)biT[index];
}
nbrPartials++;
}
}
}
}
} while (index2 > 0);

/*******************/
/* Next polynomial */
/*******************/
PolynomialIndex++;
} while (PolynomialIndex < NbrPolynomials);
PolynomialSetNbr++;
} while (true);
}

void BigIntToBigNbr(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;
}

BigInteger BigNbrToBigInt(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));
}

/************************/
/* Linear algebra phase */
/************************/
boolean LinearAlgebraPhase(int NbrPrimes, int [][] matrixB, long []
prime,
long [] biT, long [] biR, long [] biU,
int [] vectExpParity, long[][] vectLeftHandSide) {

int mask, i, j, index;
int [] rowMatrixB;

int [] matrixV = BlockLanczos(matrixB);
for (mask=1; mask != 0; mask*=2) {
LongToBigNbr(1, biT);
LongToBigNbr(1, biR);
for (i=matrixB.length-1; i>=0; i--) {
vectExpParity[i] = 0;
}
for (i=matrixB.length-1; i>=0; i--) {
if ((matrixV[i] & mask) != 0) {
NumberLength--;
MultBigNbrModN(vectLeftHandSide[i], biR, biU);
NumberLength++;
for (j=0; j<NumberLength; j++) {
biR[j] = biU[j];
}
rowMatrixB = matrixB[i];
NumberLength--;
for (j=rowMatrixB.length-1; j>=0; j--) {
vectExpParity[rowMatrixB[j]] ^= 1;
if (vectExpParity[rowMatrixB[j]] == 0) {
if (rowMatrixB[j] == 0) {
SubtractBigNbr(TestNbr, biT, biT); // Multiply biT
by -1.
}
else {
MultBigNbrByLongModN(biT, prime[rowMatrixB[j]], biT);
}
}
}
NumberLength++;
}
}
for (i=NbrPrimes-1; i>=0; i--) {
if (vectExpParity[i] != 0) {
break;
}
}
SubtractBigNbrModN(biR, biT, biR);
GcdBigNbr(biR, TestNbr2, biT);
index = 0;
if (biT[0] == 1) {
for (index=1; index<NumberLength; index++) {
if (biT[index] != 0) {
break;
}
}
}
if (index<NumberLength) { /* GCD not 1 */
for (index=0; index<NumberLength; index++) {
if (biT[index] != TestNbr2[index]) {
break;
}
}
if (index<NumberLength) { /* GCD not 1 */
return true;
}
}
}
return false;
}

// 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--;
}
}

int [] BlockLanczos(int [][] matrixB) {
int i, j, k=matrixB.length;
int oldDiagonalSSt,newDiagonalSSt;
int index, indexC, mask;
int [] matrixD = new int[32];
int [] matrixE = new int[32];
int [] matrixF = new int[32];
int [] matrixWinv = new int[32];
int [] matrixWinv1 = new int[32];
int [] matrixWinv2 = new int[32];
int [] matrixVtV0 = new int[32];
int [] matrixVt1V0 = new int[32];
int [] matrixVt2V0 = new int[32];
int [] matrixVtAV = new int[32];
int [] matrixVt1AV1 = new int[32];
int [] matrixAV = new int[k];
int [] matrixCalcParenD = new int[32];
int [] vectorIndex = new int[64];
int [] matrixV = new int[k];
int [] matrixV1 = new int[k];
int [] matrixV2 = new int[k];
int [] matrixXmY = new int[k];
int [] matrixCalc3 = new int[k]; // Matrix that holds
temporary data
int [] matrixTemp;
int [] matrixCalc1 = new int[32]; // Matrix that holds
temporary data
int [] matrixCalc2 = new int[32]; // Matrix that holds
temporary data
int [] matr;
long seed;
int Temp, Temp1;
int dimension = 0;
int stepNbr = 0;
int currentOrder, currentMask;
int dim, maxdim;
int minind, min, minanz;
int [] rowMatrixB;

newDiagonalSSt = oldDiagonalSSt = -1;

/* Initialize matrix X-Y with random data */
seed = 123456789L;
for (i=matrixXmY.length-1; i>=0; i--) {
matrixXmY[i] = (int)seed;
seed = (seed * 62089911L + 54325442L) % DosALa31_1;
matrixXmY[i] += (int)(seed * 6543265L);
seed = (seed * 62089911L + 54325442L) % DosALa31_1;
}
// Compute matrix V(0) = A(X-Y)
MultiplyAByMatrix(matrixB, matrixXmY, matrixCalc3, matrixV);
// Compute matrix Vt(0) * V(0)
MatrTranspMult(matrixV, matrixV, matrixVtV0);
while (true) {
oldDiagonalSSt = newDiagonalSSt;
stepNbr++;
// Compute matrix A * V(i)
MultiplyAByMatrix(matrixB, matrixV, matrixCalc3, matrixAV);
// Compute matrix Vt(i) * A * V(i)
MatrTranspMult(matrixV, matrixAV, matrixVtAV);

/* If Vt(i) * A * V(i) = 0, end of loop */
for (i=matrixVtAV.length-1; i>=0; i--) {
if (matrixVtAV[i] != 0) {break;}
}
if (i<0) {break;} /* End X-Y calculation loop */

/* Selection of S(i) and W(i) */

matrixTemp = matrixWinv2;
matrixWinv2 = matrixWinv1;
matrixWinv1 = matrixWinv;
matrixWinv = matrixTemp;

mask = 1;
for (j=31; j>=0; j--) {
matrixD[j] = matrixVtAV[j]; /* D = VtAV */
matrixWinv[j] = mask; /* Winv = I */
mask *= 2;
}

index = 31;
indexC = 31;
for (mask = 1; mask != 0; mask *= 2) {
if ((oldDiagonalSSt & mask) != 0) {
matrixE[index] = indexC;
matrixF[index] = mask;
index--;
}
indexC--;
}
indexC = 31;
for (mask = 1; mask != 0; mask *= 2) {
if ((oldDiagonalSSt & mask) == 0) {
matrixE[index] = indexC;
matrixF[index] = mask;
index--;
}
indexC--;
}
newDiagonalSSt = 0;
for (j=0; j<32; j++) {
currentOrder = matrixE[j];
currentMask = matrixF[j];
for (k=j; k<32; k++) {
if ((matrixD[matrixE[k]] & currentMask) != 0) {
break;
}
}
if (k<32) {
i = matrixE[k];
Temp = matrixWinv[i];
matrixWinv[i] = matrixWinv[currentOrder];
matrixWinv[currentOrder] = Temp;
Temp1 = matrixD[i];
matrixD[i] = matrixD[currentOrder];
matrixD[currentOrder] = Temp1;
newDiagonalSSt |= currentMask;
dimension++;
for (k=31; k>=0; k--) {
if (k != currentOrder) {
if ((matrixD[k] & currentMask) != 0) {
matrixWinv[k] ^= Temp;
matrixD[k] ^= Temp1;
}
}
} /* end for k */
}
else {
for (k=j; k<32; k++) {
if ((matrixWinv[matrixE[k]] & currentMask) != 0) {
break;
}
}
i = matrixE[k];
Temp = matrixWinv[i];
matrixWinv[i] = matrixWinv[currentOrder];
matrixWinv[currentOrder] = Temp;
Temp1 = matrixD[i];
matrixD[i] = matrixD[currentOrder];
matrixD[currentOrder] = Temp1;
for (k=31; k>=0; k--) {
if ((matrixWinv[k] & currentMask) != 0) {
matrixWinv[k] ^= Temp;
matrixD[k] ^= Temp1;
}
} /* end for k */
} /* end if */
} /* end for j */
/* Compute D(i), E(i) and F(i) */
if (oldDiagonalSSt != -1) { /* S=I => F=0 */
// F = -Winv(i-2) * (I - Vt(i-1)*A*V(i-1)*Winv(i-1)) * ParenD
* S*St
MatrixMultiplication(matrixVt1AV1, matrixWinv1, matrixCalc2);
index = 31; /* Add identity matrix */
for (mask = 1; mask != 0; mask *= 2) {
matrixCalc2[index] ^= mask;
index--;
}
MatrixMultiplication(matrixWinv2, matrixCalc2, matrixCalc1);
MatrixMultiplication(matrixCalc1, matrixCalcParenD, matrixF);
MatrMultBySSt(matrixF, newDiagonalSSt, matrixF);
}
// E = -Winv(i-1) * Vt(i)*A*V(i) * S*St
MatrixMultiplication(matrixWinv1, matrixVtAV, matrixE);
MatrMultBySSt(matrixE, newDiagonalSSt, matrixE);

// ParenD = Vt(i)*A*A*V(i) * S*St + Vt(i)*A*V(i)
// D = I - Winv(i) * ParenD
MatrTranspMult(matrixAV, matrixAV, matrixCalc1); //
Vt(i)*A*A*V(i)
MatrMultBySSt(matrixCalc1, newDiagonalSSt, matrixCalc1);
MatrixAddition(matrixCalc1, matrixVtAV, matrixCalcParenD);
MatrixMultiplication(matrixWinv, matrixCalcParenD, matrixD);
index = 31; /* Add identity matrix */
for (mask = 1; mask != 0; mask *= 2) {
matrixD[index] ^= mask;
index--;
}

/* Update value of X - Y */
MatrixMultiplication(matrixWinv, matrixVtV0, matrixCalc1);
MatrixMultAdd(matrixV, matrixCalc1, matrixXmY);

/* Compute value of new matrix V(i) */
// V(i+1) = A * V(i) * S * St + V(i) * D + V(i-1) * E + V(i-2) *
F
MatrMultBySSt(matrixAV, newDiagonalSSt, matrixCalc3);
MatrixMultAdd(matrixV, matrixD, matrixCalc3);
MatrixMultAdd(matrixV1, matrixE, matrixCalc3);
if (oldDiagonalSSt != -1) { // F != 0
MatrixMultAdd(matrixV2, matrixF, matrixCalc3);
}

/* Compute value of new matrix Vt(i)V0 */
if (stepNbr > 3) {
// Vt(i+1)V(0) = Dt * Vt(i)V(0) + Et * Vt(i-1)V(0) + Ft *
Vt(i-2)V(0)
MatrTranspMult(matrixD, matrixVtV0, matrixCalc1);
MatrTranspMult(matrixE, matrixVt1V0, matrixCalc2);
MatrixAddition(matrixCalc1, matrixCalc2, matrixCalc2);
if (oldDiagonalSSt != -1) { // F != 0
MatrTranspMult(matrixF, matrixVt2V0, matrixCalc1);
MatrixAddition(matrixCalc1, matrixCalc2, matrixCalc2);
}
}
else {
if (stepNbr == 1) {
MatrTranspMult(matrixCalc3, matrixV, matrixCalc2); //
V(1)t * V(0)
}
else { // if stepNbr == 2 ...
MatrTranspMult(matrixCalc3, matrixV1, matrixCalc2); //
V(2)t * V(0)
}
}
matrixTemp = matrixV2;
matrixV2 = matrixV1;
matrixV1 = matrixV;
matrixV = matrixCalc3;
matrixCalc3 = matrixTemp;
matrixTemp = matrixVt2V0;
matrixVt2V0 = matrixVt1V0;
matrixVt1V0 = matrixVtV0;
matrixVtV0 = matrixCalc2;
matrixCalc2 = matrixTemp;
matrixTemp = matrixVt1AV1;
matrixVt1AV1 = matrixVtAV;
matrixVtAV = matrixTemp;
} /* end while */

/* Find matrix V1:V2 = B * (X-Y:V) */
{
int row;
for (row = matrixB.length-1; row >= 0; row--) {
matrixV1[row] = matrixV2[row] = 0;
}
for (row = matrixB.length-1; row >= 0; row--) {
rowMatrixB = matrixB[row];
for (index = rowMatrixB.length-1; index>=0; index--) {
matrixV1[rowMatrixB[index]] ^= matrixXmY[row];
matrixV2[rowMatrixB[index]] ^= matrixV[row];
}
}
}
maxdim=64;
dim=0;
while (dim<maxdim) {
for (i=dim; i<maxdim; i++) {
matr = (i>=32?matrixV1: matrixV2);
mask = 1 << (31-(i & 31));
vectorIndex[i] = -1;
for (j=0; j<matrixV1.length; j++) {
if ((matr[j] & mask) != 0) {
vectorIndex[i] = j;
break;
}
}
}
for (i=dim; i<maxdim; i++) {
if (vectorIndex[i] < 0) {
colexchange(matrixXmY, matrixV, matrixV1, matrixV2, dim, i);
vectorIndex[i] = vectorIndex[dim];
vectorIndex[dim]=-1;
dim++;
}
}
if (dim==maxdim) {break;}
min=vectorIndex[dim];
minind=dim;
for (i=dim+1; i<maxdim; i++) {
if (vectorIndex[i] < min) {
min=vectorIndex[i];
minind=i;
}
}
minanz=0;
for (i=dim; i<maxdim; i++) {
if (vectorIndex[i] == min) {
minanz++;
}
}
if (minanz>1) {
for (i=minind+1; i<maxdim; i++) {
if (vectorIndex[i]==min) {
coladd(matrixXmY, matrixV, matrixV1, matrixV2, minind, i);
}
}
}
else {
maxdim--;
colexchange(matrixXmY, matrixV, matrixV1, matrixV2, minind,
maxdim);
}
}
dim=0; /* find linear independent solutions */
while (dim<maxdim) {
for (i=dim; i<maxdim; i++) {
matr = (i>=32?matrixXmY: matrixV);
mask = 1 << (31-(i & 31));
vectorIndex[i] = -1;
for (j=0; j<matrixV1.length; j++) {
if ((matr[j] & mask) != 0) {
vectorIndex[i] = j;
break;
}
}
}
for (i=dim; i<maxdim; i++) {
if (vectorIndex[i] < 0) {
maxdim--;
colexchange(matrixXmY, matrixV, matrixV1, matrixV2, maxdim,
i);
vectorIndex[i] = vectorIndex[maxdim];
vectorIndex[maxdim]=-1;
}
}
if (dim==maxdim) {break;}
min=vectorIndex[dim];
minind=dim;
for (i=dim+1; i<maxdim; i++) {
if (vectorIndex[i] < min) {
min=vectorIndex[i];
minind=i;
}
}
minanz=0;
for (i=dim; i<maxdim; i++) {
if (vectorIndex[i] == min) {
minanz++;
}
}
if (minanz>1) {
for (i=minind+1; i<maxdim; i++) {
if (vectorIndex[i] == min) {
coladd(matrixXmY, matrixV, matrixV1, matrixV2, minind, i);
}
}
}
else {
colexchange(matrixXmY, matrixV, matrixV1, matrixV2, minind,
dim);
dim++;
}
}
return matrixV;
}

void colexchange(int [] XmY, int [] V, int [] V1, int [] V2,
int col1, int col2) {

int i;
int mask1, mask2;
int [] matr1, matr2;

if (col1==col2) {
return;
}
mask1 = 1 << (31-(col1 & 31));
mask2 = 1 << (31-(col2 & 31));
matr1 = (col1 >= 32? V1: V2);
matr2 = (col2 >= 32? V1: V2);
for (i=V.length-1; i>=0; i--) {
if (((matr1[i] & mask1) == 0) != ((matr2[i] & mask2) == 0)) {
matr1[i] ^= mask1;
matr2[i] ^= mask2;
}
}

matr1 = (col1 >= 32? XmY: V);
matr2 = (col2 >= 32? XmY: V);
for (i=V.length-1; i>=0; i--) {
if (((matr1[i] & mask1) == 0) != ((matr2[i] & mask2) == 0)) {
matr1[i] ^= mask1;
matr2[i] ^= mask2;
}
}
}

void coladd(int [] XmY, int [] V, int [] V1, int [] V2,
int col1, int col2) {

int i;
int mask1, mask2;
int [] matr1, matr2;

if (col1==col2) {
return;
}
mask1 = 1 << (31-(col1 & 31));
mask2 = 1 << (31-(col2 & 31));
matr1 = (col1 >= 32? V1: V2);
matr2 = (col2 >= 32? V1: V2);
for (i=V.length-1; i>=0; i--) {
if ((matr1[i] & mask1) != 0) {
matr2[i] ^= mask2;
}
}

matr1 = (col1 >= 32? XmY: V);
matr2 = (col2 >= 32? XmY: V);
for (i=V.length-1; i>=0; i--) {
if ((matr1[i] & mask1) != 0) {
matr2[i] ^= mask2;
}
}
}

/* Multiply the transpose of a binary matrix of length n x 32 by */
/* another binary matrix of length n x 32 */
/* The product matrix has size 32 x 32 */
void MatrTranspMult(int [] LeftMatr, int [] RightMatr, int []
ProdMatr) {
int prodMatr;
int matrLength = LeftMatr.length;
int row, col;
int iMask = 1;
for (col = 31; col >= 0; col--) {
prodMatr = 0;
for (row = 0; row < matrLength; row++) {
if ((LeftMatr[row] & iMask) != 0) {
prodMatr ^= RightMatr[row];
}
}
ProdMatr[col] = prodMatr;
iMask *= 2;
}
}

void MatrixAddition(int [] leftMatr, int [] rightMatr, int []
sumMatr) {
for (int row = leftMatr.length-1; row >= 0 ; row--) {
sumMatr[row] = leftMatr[row] ^ rightMatr[row];
}
}

/* Multiply binary matrices of length m x 32 by 32 x 32 */
/* The product matrix has size m x 32. Then add it to a m x 32
matrix. */
void MatrixMultAdd(int [] LeftMatr, int [] RightMatr, int []
ProdMatr) {
int leftMatr;
int matrLength = LeftMatr.length;
int prodMatr;
int row, col;
for (row = 0; row < matrLength; row++) {
prodMatr = ProdMatr[row];
leftMatr = LeftMatr[row];
col = 0;
while (leftMatr != 0) {
if (leftMatr < 0) {
prodMatr ^= RightMatr[col];
}
leftMatr *= 2;
col++;
}
ProdMatr[row] = prodMatr;
}
}

void MatrMultBySSt(int [] Matr, int diagS, int [] Prod) {
for (int row = Matr.length-1; row >= 0 ; row--) {
Prod[row] = diagS & Matr[row];
}
}

/* Multiply binary matrices of length m x 32 by 32 x 32 */
/* The product matrix has size m x 32 */
void MatrixMultiplication(int [] LeftMatr, int [] RightMatr, int []
ProdMatr) {
int leftMatr;
int matrLength = LeftMatr.length;
int prodMatr;
int row, col;
for (row = 0; row < matrLength; row++) {
prodMatr = 0;
leftMatr = LeftMatr[row];
col = 0;
while (leftMatr != 0) {
if (leftMatr < 0) {
prodMatr ^= RightMatr[col];
}
leftMatr *= 2;
col++;
}
ProdMatr[row] = prodMatr;
}
}

/* Compute Bt * B * input matrix where B is the matrix that holds
the */
/* factorization relations */
void MultiplyAByMatrix(int [][] matrixB, int [] Matr, int []
TempMatr, int [] ProdMatr) {
int index;
int prodMatr;
int row;
int [] rowMatrixB;

/* Compute TempMatr = B * Matr */
for (row = matrixB.length-1; row >= 0; row--) {
TempMatr[row] = 0;
}
for (row = matrixB.length-1; row >= 0; row--) {
rowMatrixB = matrixB[row];
for (index = rowMatrixB.length-1; index>=0; index--) {
TempMatr[rowMatrixB[index]] ^= Matr[row];
}
}

/* Compute ProdMatr = Bt * TempMatr */
for (row = matrixB.length-1; row >= 0; row--) {
prodMatr = 0;
rowMatrixB = matrixB[row];
for (index = rowMatrixB.length-1; index>=0; index--) {
prodMatr ^= TempMatr[rowMatrixB[index]];
}
ProdMatr[row] = prodMatr;
}
}

int EraseSingletons(int [][] matrixB, long [][] vectLeftHandSide,
int [] vectExpParity) {
int i,j,delta;
int [] rowMatrixB;
// Find singletons in matrixB.
for (i=matrixB.length-1; i>=0; i--) {
vectExpParity[i] = 0;
}
for (i=matrixB.length-1; i>=0; i--) {
rowMatrixB = matrixB[i];
for (j=rowMatrixB.length-1; j>=0; j--) {
vectExpParity[rowMatrixB[j]]++;
}
}
delta = 0;
// Erase singletons from matrixB.
for (i=0; i<matrixB.length; i++) {
rowMatrixB = matrixB[i];
for (j=rowMatrixB.length-1; j>=0; j--) {
if (vectExpParity[rowMatrixB[j]] == 1) { // Singleton found
delta++;
break;
}
}
if (j<0) { // Singleton not
found
matrixB[i-delta] = matrixB[i];
vectLeftHandSide[i-delta] = vectLeftHandSide[i];
}
}
return delta;
}

boolean InsertNewRelation(long [] biR, long [] biT, long [] biU,
int nbrColumns, int [][] matrixB, int [] rowMatrixB, int pp,
long [][] vectLeftHandSide) {
/* Convert negative numbers to the range 0 <= n < TestNbr */

int i,k,index;

if ((TestNbr[0] & 1) == 0) {
DivBigNbrByLong(TestNbr,2,TestNbr);
for (i=0; i<NumberLength; i++) {
biT[i] = 0;
}
AddBigNbrModN(biR, biT, biR);
NumberLength--;
ModInvBigNbr(biR, biT, TestNbr);
NumberLength++;
MultBigNbrByLong(TestNbr,2,TestNbr);
}
else {
ModInvBigNbr(biR, biT, TestNbr);
}
if ((biU[NumberLength-1] & 0x80000000L) != 0) {
AddBigNbr(biU, TestNbr, biU);
}

// Compute biU / biR (mod TestNbr)

NumberLength--;
MultBigNbrModN(biU, biT, biR);
NumberLength++;

// Insert it only if it is different from previous relations.

if (pp == 0) { // First relation. Insert it always.
matrixB[pp] = new int [nbrColumns]; /* Add relation to matrix B
*/
for (k=0; k<nbrColumns; k++) {
matrixB[0][k] = rowMatrixB[k];
}
vectLeftHandSide[0] = new long[NumberLength];
for (index=0; index<NumberLength; index++) {
vectLeftHandSide[0][index] = biR[index];
}
return true;
}
for (i=0; i<pp; i++) {
if (nbrColumns > matrixB[i].length) {
continue;
}
if (nbrColumns == matrixB[i].length) {
for (k=0; k<nbrColumns; k++) {
if (rowMatrixB[k] != matrixB[i][k]) {
break;
}
}
if (k==nbrColumns) {
return false; // Do not insert same relation.
}
if (rowMatrixB[k] > matrixB[i][k]) {
continue;
}
}
for (k=pp-1; k>=i; k--) {
matrixB[k+1] = matrixB[k];
vectLeftHandSide[k+1] = vectLeftHandSide[k];
}
break;
}
matrixB[i] = new int [nbrColumns]; // Add relation to matrix B.
for (k=0; k<nbrColumns; k++) {
matrixB[i][k] = rowMatrixB[k];
}
vectLeftHandSide[i] = new long[NumberLength];
for (index=0; index<NumberLength; index++) {
vectLeftHandSide[i][index] = biR[index];
}
return true;
}

/
***********************************************************************/
/* 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;
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 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);
}

long modInv(long NbrMod, long currentPrime) {
long QQ, T1, T3;
long U1 = 1;
long U3 = NbrMod;
long V1 = 0;
long V3 = currentPrime;
while (V3 != 0) {
QQ = U3/V3;
T1 = U1 - V1*QQ;
T3 = U3 - V3*QQ;
U1 = V1; U3 = V3;
V1 = T1; V3 = T3;
}
return (U1+currentPrime)%currentPrime;
}

long modPow(long NbrMod, long Expon, long currentPrime) {
long Power = 1;
long Square = NbrMod;
while (Expon != 0) {
if ((Expon & 1) == 1) {
Power = (Power * Square) % currentPrime;
}
Square = (Square * Square) % currentPrime;
Expon /= 2;
}
return Power;
}

BigInteger factorizeecm(BigInteger N) {

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

int I,J,Pass,Qaux;
long L1,L2,LS,P,Q,IP;
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 [] 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 [] 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 [] sieve = new byte[23100];
byte [] sieve2310 = new byte[2310];
int [] sieveidx = new int[480];
int i, j, u;

// added JGW
if (N.bitLength()>32000)
return N;

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

// added JGW
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');

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

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 && EC<5000);

return BigNbrToBigInt(GD);
}

/**********************/
/* Auxiliary routines */
/**********************/
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);
}

/* 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 */
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);
}

// Perform JS <- JS * JW
void JS_JW(int PK, int PL, int PM, int P) {
int I,J,K;
for (I=0; I<PL; I++) {
for (J=0; J<PL; J++) {
K = (I+J)%PK;
MontgomeryMult(aiJS[I], aiJW[J], biTmp);
AddBigNbrModN(aiJX[K], biTmp, aiJX[K]);
}
}
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJX[I][J];
aiJX[I][J] = 0;
}
}
NormalizeJS(PK, PL, PM, P);
}

// Normalize coefficient of JS
void NormalizeJS(int PK, int PL, int PM, int P) {
int I,J;
for (I=PL; I<PK; I++) {
if (BigNbrIsZERO(aiJS[I]) == false) {
for (J=0; J<NumberLength; J++) {
biT[J] = aiJS[I][J];
}
for (J=1; J<P; J++) {
SubtractBigNbrModN(aiJS[I-J*PM], biT, aiJS[I-J*PM]);
}
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = 0;
}
}
}
}

// Perform JS <- JS ^ 2
void JS_2(int PK, int PL, int PM, int P) {
int I,J,K;
for (I=0; I<PL; I++) {
K = 2*I%PK;
MontgomeryMult(aiJS[I], aiJS[I], biTmp);
AddBigNbrModN(aiJX[K], biTmp, aiJX[K]);
AddBigNbrModN(aiJS[I], aiJS[I], biT);
for (J=I+1; J<PL; J++) {
K = (I+J)%PK;
MontgomeryMult(biT, aiJS[J], biTmp);
AddBigNbrModN(aiJX[K], biTmp, aiJX[K]);
}
}
for (I=0; I<PK; I++) {
for (J=0; J<NumberLength; J++) {
aiJS[I][J] = aiJX[I][J];
aiJX[I][J] = 0;
}
}
NormalizeJS(PK, PL, PM, P);
}

void JacobiSum(int A, int B, int P, int PK, int PL, int PM, int Q) {
int I,J,K;

for (I=0; I<PL; I++) {
for (J=0; J<NumberLength; J++) {
aiJ0[I][J] = 0;
}
}
for (I=1; I<=Q-2; I++) {
J = (A*I+B*aiF[I])%PK;
if (J<PL) {
AddBigNbrModN(aiJ0[J], MontgomeryMultR1, aiJ0[J]);
}
else {
for (K=1; K<P; K++) {
SubtractBigNbrModN(aiJ0[J-K*PM], MontgomeryMultR1, aiJ0[J-
K*PM]);
}
}
}
}

// Perform JS <- JS ^ E
void JS_E(int PK, int PL, int PM, int P) {
int I, J, K;
long Mask;

for (I=NumberLength-1; I>0; I--) {
if (biExp[I] != 0) {break;}
}
if (I == 0 && biExp[0] == 1) {return;} // Return if E == 1
for (K=0; K<PL; K++) {
for (J=0; J<NumberLength; J++) {
aiJW[K][J] = aiJS[K][J];
}
}
Mask = DosALa31;
while (true) {
if ((biExp[I] & Mask) != 0) {break;}
Mask /= 2;
}
do {
JS_2(PK, PL, PM, P);
Mask /= 2;
if (Mask == 0) {
Mask = DosALa31;
I--;
}
if ((biExp[I] & Mask) != 0) {
JS_JW(PK, PL, PM, P);
}
} while (I>0 || Mask != 1);
}

// Normalize coefficient of JW
void NormalizeJW(int PK, int PL, int PM, int P) {
int I,J;
for (I=PL; I<PK; I++) {
if (BigNbrIsZERO(aiJW[I]) == false) {
for (J=0; J<NumberLength; J++) {
biT[J] = aiJW[I][J];
}
for (J=1; J<P; J++) {
SubtractBigNbrModN(aiJW[I-J*PM], biT, aiJW[I-J*PM]);
}
for (J=0; J<NumberLength; J++) {
aiJW[I][J] = 0;
}
}
}
}

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

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);
}

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 MontDig, Prd;
int MontgomeryMultN = (int)this.MontgomeryMultN;
long TestNbr[] = this.TestNbr;
int NumberLength = this.NumberLength;

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 MultBigNbrByLong(long Nbr1[], long Nbr2, long Prod[]) {
int NumberLength = this.NumberLength;
long Cy, Pr;
int i;
Cy = 0;
for (i=0; i<NumberLength; i++) {
Pr = (Cy & MaxUInt) + Nbr2*Nbr1[i];
Cy = (Cy >>> 32) + (Pr >>> 32);
Prod[i] = Pr & MaxUInt;
}
}

void BigNbrModN(long Nbr[], int Length, long Mod[]) {
int i,j;
for (i=0; i<NumberLength; i++) {
Mod[i] = Nbr[i + Length - NumberLength];
}
Mod[i] = 0;
AdjustModN(Mod);
for (i=Length - NumberLength - 1; i>=0; i--) {
for (j=NumberLength; j>0; j--) {
Mod[j] = Mod[j-1];
}
Mod[0] = Nbr[i];
AdjustModN(Mod);
}
}

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 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;
}
}
}

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;
}
}

// Compare Nbr1^2 vs. Nbr2
int CompareSquare(long Nbr1[], long Nbr2[]) {
int I,k;

for (I=NumberLength-1; I>0; I--) {
if (Nbr1[I] != 0) {break;}
}
k = NumberLength/2;
if (NumberLength%2 == 0) {
if (I >= k) {return 1;} // Nbr1^2 >
Nbr2
if (I < k - 1 || biS[k - 1] < 65536) {return -1;} // Nbr1^2 <
Nbr2
}
else {
if (I < k) {return -1;} // Nbr1^2 <
Nbr2
if (I > k || biS[k] >= 65536) {return 1;} // Nbr1^2 >
Nbr2
}
MultBigNbr(biS, biS, biTmp);
SubtractBigNbr(biTmp, TestNbr, biTmp);
if (BigNbrIsZERO(biTmp) == true) {return 0;} // Nbr1^2 ==
Nbr2
if (biTmp[NumberLength-1] >= 0) {return 1;} // Nbr1^2 >
Nbr2
return -1; // Nbr1^2 <
Nbr2
}

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 MultBigNbr(long Nbr1[], long Nbr2[], long Prod[]) {
int NumberLength = this.NumberLength;
long Cy, Pr;
int i,j;
Cy = Pr = 0;
for (i=0; i<NumberLength; i++) {
Pr = Cy & MaxUInt;
Cy >>>= 32;
for (j=0; j<=i; j++) {
Pr += Nbr1[j]*Nbr2[i-j];
Cy += (Pr >>> 32);
Pr &= MaxUInt;
}
Prod[i] = Pr;
}
}

long BigNbrModLong(long Nbr1[], long Nbr2) {
int i;
long Rem = 0;

for (i=NumberLength-1; i>=0; i--) {
Rem = ((Rem << 32) + Nbr1[i]) % Nbr2;
}
return Rem;
}

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 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 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);
}
}

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

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);
}
}

long RemDivBigNbrByLong(long Dividend[], long Divisor) {
int i;
long Divid, Rem = 0;

if (Divisor < 0) {
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;
}
return Rem;
}

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;
}
}

BigInteger factorizewe(BigInteger n) {
BigInteger T = new BigInteger("1");
BigInteger a = new BigInteger("2");
BigInteger b = new BigInteger("1");
BigInteger wanless = new BigInteger("2");

// workaround - apparent java bug in modPow - JGW
if (n.compareTo(TWO) < 0)
return ONE;

if (n.remainder(TWO).compareTo(ZERO) == 0)
return TWO;
// end workaround

while (wanless.compareTo(n) < 0)
wanless = wanless.multiply(TWO);

while (a.multiply(a).compareTo(n) <= 0 && (T.compareTo(ONE) == 0 ||
T.compareTo(n) == 0)) {
b = a.modPow(wanless, n);
T = n.gcd(b.modPow(n, n).subtract(b));
a = a.add(TWO);
}
return T;
}

}

[bearnol]

... in case it's of (more) general interest.
Two points to note:
a) The automatic, fast finding of algebraic factors
b) The input format (roughly) as follows:
"1) You can either enter just a (denary) number in the normal way:
eg, 123 or 7654321 or

2) you can enter a formula using the following operators by prefixing
with f or F:
+ for addition
- for subtraction
* for multiplication
/ for integer division
% for modulo division
^ or # for exponentiation
[ or { or ( for open parenthesis (nestable)
] or } or ) for close parenthesis (nestable)
eg, f2#700+1 or F2#(2#12)+1 or

3) you can enter a randomly generated number of approximate size 2^n,
by using the syntax rn or Rn, where n is an integer:
eg, r200 or R1000
"

[bearnol]

Hmmm..., it would appear that just posting superfac9.java as a text
message, mangles the code slightly, so I've also uploaded it as a file
to this group. I'd advise people to try using that by preference
(though hopefully it's nice to have the text, for the record, too)
J