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