64)
return euclidModInverse(k);
int t = inverseMod32(value[offset+intLen-1]);
if (k < 33) {
t = (k == 32 ? t : t & ((1 << k) - 1));
return new MutableBigInteger(t);
}
long pLong = (value[offset+intLen-1] & LONG_MASK);
if (intLen > 1)
pLong |= ((long)value[offset+intLen-2] << 32);
long tLong = t & LONG_MASK;
tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
MutableBigInteger result = new MutableBigInteger(new int[2]);
result.value[0] = (int)(tLong >>> 32);
result.value[1] = (int)tLong;
result.intLen = 2;
result.normalize();
return result;
}
/**
* Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
*/
static int inverseMod32(int val) {
// Newton's iteration!
int t = val;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
return t;
}
/**
* Returns the multiplicative inverse of val mod 2^64. Assumes val is odd.
*/
static long inverseMod64(long val) {
// Newton's iteration!
long t = val;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
assert(t * val == 1);
return t;
}
/**
* Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
*/
static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
// Copy the mod to protect original
return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
}
/**
* Calculate the multiplicative inverse of this mod mod, where mod is odd.
* This and mod are not changed by the calculation.
*
* This method implements an algorithm due to Richard Schroeppel, that uses
* the same intermediate representation as Montgomery Reduction
* ("Montgomery Form"). The algorithm is described in an unpublished
* manuscript entitled "Fast Modular Reciprocals."
*/
private MutableBigInteger modInverse(MutableBigInteger mod) {
MutableBigInteger p = new MutableBigInteger(mod);
MutableBigInteger f = new MutableBigInteger(this);
MutableBigInteger g = new MutableBigInteger(p);
SignedMutableBigInteger c = new SignedMutableBigInteger(1);
SignedMutableBigInteger d = new SignedMutableBigInteger();
MutableBigInteger temp = null;
SignedMutableBigInteger sTemp = null;
int k = 0;
// Right shift f k times until odd, left shift d k times
if (f.isEven()) {
int trailingZeros = f.getLowestSetBit();
f.rightShift(trailingZeros);
d.leftShift(trailingZeros);
k = trailingZeros;
}
// The Almost Inverse Algorithm
while (!f.isOne()) {
// If gcd(f, g) != 1, number is not invertible modulo mod
if (f.isZero())
throw new ArithmeticException("BigInteger not invertible.");
// If f < g exchange f, g and c, d
if (f.compare(g) < 0) {
temp = f; f = g; g = temp;
sTemp = d; d = c; c = sTemp;
}
// If f == g (mod 4)
if (((f.value[f.offset + f.intLen - 1] ^
g.value[g.offset + g.intLen - 1]) & 3) == 0) {
f.subtract(g);
c.signedSubtract(d);
} else { // If f != g (mod 4)
f.add(g);
c.signedAdd(d);
}
// Right shift f k times until odd, left shift d k times
int trailingZeros = f.getLowestSetBit();
f.rightShift(trailingZeros);
d.leftShift(trailingZeros);
k += trailingZeros;
}
while (c.sign < 0)
c.signedAdd(p);
return fixup(c, p, k);
}
/**
* The Fixup Algorithm
* Calculates X such that X = C * 2^(-k) (mod P)
* Assumes C> 5; i < numWords; i++) {
// V = R * c (mod 2^j)
int v = r * c.value[c.offset + c.intLen-1];
// c = c + (v * p)
p.mul(v, temp);
c.add(temp);
// c = c / 2^j
c.intLen--;
}
int numBits = k & 0x1f;
if (numBits != 0) {
// V = R * c (mod 2^j)
int v = r * c.value[c.offset + c.intLen-1];
v &= ((1<= 0)
c.subtract(p);
return c;
}
/**
* Uses the extended Euclidean algorithm to compute the modInverse of base
* mod a modulus that is a power of 2. The modulus is 2^k.
*/
MutableBigInteger euclidModInverse(int k) {
MutableBigInteger b = new MutableBigInteger(1);
b.leftShift(k);
MutableBigInteger mod = new MutableBigInteger(b);
MutableBigInteger a = new MutableBigInteger(this);
MutableBigInteger q = new MutableBigInteger();
MutableBigInteger r = b.divide(a, q);
MutableBigInteger swapper = b;
// swap b & r
b = r;
r = swapper;
MutableBigInteger t1 = new MutableBigInteger(q);
MutableBigInteger t0 = new MutableBigInteger(1);
MutableBigInteger temp = new MutableBigInteger();
while (!b.isOne()) {
r = a.divide(b, q);
if (r.intLen == 0)
throw new ArithmeticException("BigInteger not invertible.");
swapper = r;
a = swapper;
if (q.intLen == 1)
t1.mul(q.value[q.offset], temp);
else
q.multiply(t1, temp);
swapper = q;
q = temp;
temp = swapper;
t0.add(q);
if (a.isOne())
return t0;
r = b.divide(a, q);
if (r.intLen == 0)
throw new ArithmeticException("BigInteger not invertible.");
swapper = b;
b = r;
if (q.intLen == 1)
t0.mul(q.value[q.offset], temp);
else
q.multiply(t0, temp);
swapper = q; q = temp; temp = swapper;
t1.add(q);
}
mod.subtract(t1);
return mod;
}
}