📄 gcd.java
字号:
package javax.math.factorization;
/**
* <p>Title: Factorization Library</p>
* <p>Description: </p>
* <p>Copyright: Copyright (c) 2004</p>
* <p>Company: </p>
* @author Vladimir Silva
* @version 1.0
*/
/**
* I am truly sorry about all this awful spagethi code
* I got it from: Dario Alejandro Alpern (Buenos Aires - Argentina)
* Last updated October 16th, 2003. See http://www.alpertron.com.ar/ECM.HTM
* I am not to blame for its crappy design and implementation, but....
* what the hell. it works!
* I modified it though, to make it less horrible (extracted the pieces i needed)
* If you are a brave soul, try to port this stuff to an object oriented, well-designed
* structure. I tried, but just gave up....
*/
public class GCD {
static final long CalcAuxGcdU[] = new long[Constants.MP_INITIAL_SIZE]; /* Used inside GCD calculations in */
static final long CalcAuxGcdV[] = new long[Constants.MP_INITIAL_SIZE]; /* multiple precision numbers */
static final long CalcAuxGcdT[] = new long[Constants.MP_INITIAL_SIZE];
/***********************************************************************/
/* NAME: ModInvBigNbr */
/* */
/* PURPOSE: Find the inverse multiplicative modulo v. */
/* */
/* The algorithm terminates with u1 = u^(-1) MOD v. */
/***********************************************************************/
static void ModInvBigNbr(long[] a, long[] inv, long[] b) {
int i, j;
int NumberLength = Factorizer.NumberLength;
int Dif, E;
int st1, st2;
long Yaa, Yab; // 2^E * A' = Yaa A + Yab B
long Yba, Ybb; // 2^E * B' = Yba A + Ybb B
long Ygb0; // 2^E * Mu' = Yaa Mu + Yab Gamma + Ymb0 B0
long Ymb0; // 2^E * Gamma' = Yba Mu + Ybb Gamma + Ygb0 B0
int Iaa, Iab, Iba, Ibb;
long Tmp1, Tmp2, Tmp3, Tmp4, Tmp5;
int B0l = (int) b[0];
int invB0l;
int Al, Bl, T1, Gl, Ml, P;
long T;
long Cy1, Cy2, Cy3, Cy4;
int Yaah, Yabh, Ybah, Ybbh;
int Ymb0h, Ygb0h;
long Pr1, Pr2, Pr3, Pr4, Pr5, Pr6, Pr7;
long[] CalcAuxModInvA = CalcAuxGcdU;
long[] CalcAuxModInvB = CalcAuxGcdV;
long[] CalcAuxModInvMu = inv;
long[] CalcAuxModInvGamma = 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) {
BigNumOps.SubtractBigNbr(b, CalcAuxModInvMu, CalcAuxModInvMu);
}
if ( (int) CalcAuxModInvMu[i = NumberLength - 1] < 0) {
BigNumOps.AddBigNbr(b, CalcAuxModInvMu, CalcAuxModInvMu);
}
for (; i >= 0; i--) {
if (b[i] != CalcAuxModInvMu[i])
break;
}
if (i < 0 || b[i] < CalcAuxModInvMu[i]) { // If B < Mu
BigNumOps.SubtractBigNbr(CalcAuxModInvMu, b, CalcAuxModInvMu); // Mu <- Mu - B
}
System.arraycopy(CalcAuxModInvMu, 0, inv, 0, NumberLength);
}
// 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.
static void GcdBigNbr(long Nbr1[], long Nbr2[], long Gcd[]) {
int i,k;
int NumberLength = Factorizer.NumberLength;
System.arraycopy(Nbr1,0, GCD.CalcAuxGcdU,0,NumberLength);
System.arraycopy(Nbr2,0, GCD.CalcAuxGcdV,0,NumberLength);
for (i=0; i<NumberLength; i++) {
if (GCD.CalcAuxGcdU[i] != 0) {break;}
}
if (i == NumberLength) {
System.arraycopy(GCD.CalcAuxGcdV,0,Gcd,0,NumberLength);
return;
}
for (i=0; i<NumberLength; i++) {
if (GCD.CalcAuxGcdV[i] != 0) {break;}
}
if (i == NumberLength) {
System.arraycopy(GCD.CalcAuxGcdU,0,Gcd,0,NumberLength);
return;
}
if (GCD.CalcAuxGcdU[NumberLength-1] >= Constants.DosALa31) {
BigNumOps.ChSignBigNbr(GCD.CalcAuxGcdU);
}
if (GCD.CalcAuxGcdV[NumberLength-1] >= Constants.DosALa31) {
BigNumOps.ChSignBigNbr(GCD.CalcAuxGcdV);
}
k=0;
while ((GCD.CalcAuxGcdU[0] & 1)==0 && (GCD.CalcAuxGcdV[0] & 1)==0) { // Step 1
k++;
BigNumOps.DivBigNbrByLong(GCD.CalcAuxGcdU,2L,GCD.CalcAuxGcdU);
BigNumOps.DivBigNbrByLong(GCD.CalcAuxGcdV,2L,GCD.CalcAuxGcdV);
}
if ((CalcAuxGcdU[0] & 1)==1) { // Step 2
System.arraycopy(GCD.CalcAuxGcdV,0,GCD.CalcAuxGcdT,0,NumberLength);
BigNumOps.ChSignBigNbr(GCD.CalcAuxGcdT);
}
else {
System.arraycopy(GCD.CalcAuxGcdU,0,GCD.CalcAuxGcdT,0,NumberLength);
}
do {
while ((GCD.CalcAuxGcdT[0] & 1)==0) { // Step 4
BigNumOps.DivBigNbrByLong(GCD.CalcAuxGcdT,2L,GCD.CalcAuxGcdT); // Step 3
}
if (CalcAuxGcdT[NumberLength-1] < Constants.DosALa31) { // Step 5
System.arraycopy(GCD.CalcAuxGcdT,0,GCD.CalcAuxGcdU,0,NumberLength);
}
else {
System.arraycopy(GCD.CalcAuxGcdT,0,GCD.CalcAuxGcdV,0,NumberLength);
BigNumOps.ChSignBigNbr(GCD.CalcAuxGcdV);
}
BigNumOps.SubtractBigNbr(GCD.CalcAuxGcdU,GCD.CalcAuxGcdV,GCD.CalcAuxGcdT); // Step 6
for (i=0; i<NumberLength; i++) {
if (GCD.CalcAuxGcdT[i] != 0) {break;}
}
} while (i != NumberLength);
System.arraycopy(GCD.CalcAuxGcdU,0,Gcd,0,NumberLength); // Step 7
while (k>0) {
BigNumOps.AddBigNbr(Gcd,Gcd,Gcd);
k--;
}
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -