All Downloads are FREE. Search and download functionalities are using the official Maven repository.

com.google.gwt.emul.java.math.Division Maven / Gradle / Ivy

/*
 * Copyright 2009 Google Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may not
 * use this file except in compliance with the License. You may obtain a copy of
 * the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations under
 * the License.
 */

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements. See the NOTICE file distributed with this
 * work for additional information regarding copyright ownership. The ASF
 * licenses this file to You under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations under
 * the License.
 *
 * INCLUDES MODIFICATIONS BY RICHARD ZSCHECH AS WELL AS GOOGLE.
 */
package java.math;

/**
 * Static library that provides all operations related with division and modular
 * arithmetic to {@link BigInteger}. Some methods are provided in both mutable
 * and immutable way. There are several variants provided listed below:
 *
 * 
  • Division
    • * {@link BigInteger} division and remainder by {@link BigInteger}.
    • * {@link BigInteger} division and remainder by {@code int}.
    • gcd * between {@link BigInteger} numbers.
  • Modular * arithmetic
    • Modular exponentiation between * {@link BigInteger} numbers.
    • Modular inverse of a {@link BigInteger} * numbers.
*/ class Division { /** * Divides the array 'a' by the array 'b' and gets the quotient and the * remainder. Implements the Knuth's division algorithm. See D. Knuth, The Art * of Computer Programming, vol. 2. Steps D1-D8 correspond the steps in the * algorithm description. * * @param quot the quotient * @param quotLength the quotient's length * @param a the dividend * @param aLength the dividend's length * @param b the divisor * @param bLength the divisor's length * @return the remainder */ static int[] divide(int quot[], int quotLength, int a[], int aLength, int b[], int bLength) { int normA[] = new int[aLength + 1]; // the normalized dividend // an extra byte is needed for correct shift int normB[] = new int[bLength + 1]; // the normalized divisor; int normBLength = bLength; /* * Step D1: normalize a and b and put the results to a1 and b1 the * normalized divisor's first digit must be >= 2^31 */ int divisorShift = Integer.numberOfLeadingZeros(b[bLength - 1]); if (divisorShift != 0) { BitLevel.shiftLeft(normB, b, 0, divisorShift); BitLevel.shiftLeft(normA, a, 0, divisorShift); } else { System.arraycopy(a, 0, normA, 0, aLength); System.arraycopy(b, 0, normB, 0, bLength); } int firstDivisorDigit = normB[normBLength - 1]; // Step D2: set the quotient index int i = quotLength - 1; int j = aLength; while (i >= 0) { // Step D3: calculate a guess digit guessDigit int guessDigit = 0; if (normA[j] == firstDivisorDigit) { // set guessDigit to the largest unsigned int value guessDigit = -1; } else { long product = (((normA[j] & 0xffffffffL) << 32) + (normA[j - 1] & 0xffffffffL)); long res = Division.divideLongByInt(product, firstDivisorDigit); guessDigit = (int) res; // the quotient of divideLongByInt int rem = (int) (res >> 32); // the remainder of // divideLongByInt // decrease guessDigit by 1 while leftHand > rightHand if (guessDigit != 0) { long leftHand = 0; long rightHand = 0; boolean rOverflowed = false; guessDigit++; // to have the proper value in the loop // below do { guessDigit--; if (rOverflowed) { break; } // leftHand always fits in an unsigned long leftHand = (guessDigit & 0xffffffffL) * (normB[normBLength - 2] & 0xffffffffL); /* * rightHand can overflow; in this case the loop condition will be * true in the next step of the loop */ rightHand = ((long) rem << 32) + (normA[j - 2] & 0xffffffffL); long longR = (rem & 0xffffffffL) + (firstDivisorDigit & 0xffffffffL); /* * checks that longR does not fit in an unsigned int; this ensures * that rightHand will overflow unsigned long in the next step */ if (Integer.numberOfLeadingZeros((int) (longR >>> 32)) < 32) { rOverflowed = true; } else { rem = (int) longR; } } while (((leftHand ^ 0x8000000000000000L) > (rightHand ^ 0x8000000000000000L))); } } // Step D4: multiply normB by guessDigit and subtract the production // from normA. if (guessDigit != 0) { int borrow = Division.multiplyAndSubtract(normA, j - normBLength, normB, normBLength, guessDigit); // Step D5: check the borrow if (borrow != 0) { // Step D6: compensating addition guessDigit--; long carry = 0; for (int k = 0; k < normBLength; k++) { carry += (normA[j - normBLength + k] & 0xffffffffL) + (normB[k] & 0xffffffffL); normA[j - normBLength + k] = (int) carry; carry >>>= 32; } } } if (quot != null) { quot[i] = guessDigit; } // Step D7 j--; i--; } /* * Step D8: we got the remainder in normA. Denormalize it id needed */ if (divisorShift != 0) { // reuse normB BitLevel.shiftRight(normB, normBLength, normA, 0, divisorShift); return normB; } System.arraycopy(normA, 0, normB, 0, bLength); return normA; } /** * Computes the quotient and the remainder after a division by an {@code int} * number. * * @return an array of the form {@code [quotient, remainder]}. */ static BigInteger[] divideAndRemainderByInteger(BigInteger val, int divisor, int divisorSign) { // res[0] is a quotient and res[1] is a remainder: int[] valDigits = val.digits; int valLen = val.numberLength; int valSign = val.sign; if (valLen == 1) { long a = (valDigits[0] & 0xffffffffL); long b = (divisor & 0xffffffffL); long quo = a / b; long rem = a % b; if (valSign != divisorSign) { quo = -quo; } if (valSign < 0) { rem = -rem; } return new BigInteger[] {BigInteger.valueOf(quo), BigInteger.valueOf(rem)}; } int quotientLength = valLen; int quotientSign = ((valSign == divisorSign) ? 1 : -1); int quotientDigits[] = new int[quotientLength]; int remainderDigits[]; remainderDigits = new int[] {Division.divideArrayByInt(quotientDigits, valDigits, valLen, divisor)}; BigInteger result0 = new BigInteger(quotientSign, quotientLength, quotientDigits); BigInteger result1 = new BigInteger(valSign, 1, remainderDigits); result0.cutOffLeadingZeroes(); result1.cutOffLeadingZeroes(); return new BigInteger[] {result0, result1}; } /** * Divides an array by an integer value. Implements the Knuth's division * algorithm. See D. Knuth, The Art of Computer Programming, vol. 2. * * @param dest the quotient * @param src the dividend * @param srcLength the length of the dividend * @param divisor the divisor * @return remainder */ static int divideArrayByInt(int dest[], int src[], final int srcLength, final int divisor) { long rem = 0; long bLong = divisor & 0xffffffffL; for (int i = srcLength - 1; i >= 0; i--) { long temp = (rem << 32) | (src[i] & 0xffffffffL); long quot; if (temp >= 0) { quot = (temp / bLong); rem = (temp % bLong); } else { /* * make the dividend positive shifting it right by 1 bit then get the * quotient an remainder and correct them properly */ long aPos = temp >>> 1; long bPos = divisor >>> 1; quot = aPos / bPos; rem = aPos % bPos; // double the remainder and add 1 if a is odd rem = (rem << 1) + (temp & 1); if ((divisor & 1) != 0) { // the divisor is odd if (quot <= rem) { rem -= quot; } else { if (quot - rem <= bLong) { rem += bLong - quot; quot -= 1; } else { rem += (bLong << 1) - quot; quot -= 2; } } } } dest[i] = (int) (quot & 0xffffffffL); } return (int) rem; } /** * Divides an unsigned long a by an unsigned int b. It is supposed that the * most significant bit of b is set to 1, i.e. b < 0 * * @param a the dividend * @param b the divisor * @return the long value containing the unsigned integer remainder in the * left half and the unsigned integer quotient in the right half */ static long divideLongByInt(long a, int b) { long quot; long rem; long bLong = b & 0xffffffffL; if (a >= 0) { quot = (a / bLong); rem = (a % bLong); } else { /* * Make the dividend positive shifting it right by 1 bit then get the * quotient an remainder and correct them properly */ long aPos = a >>> 1; long bPos = b >>> 1; quot = aPos / bPos; rem = aPos % bPos; // double the remainder and add 1 if a is odd rem = (rem << 1) + (a & 1); if ((b & 1) != 0) { // the divisor is odd if (quot <= rem) { rem -= quot; } else { if (quot - rem <= bLong) { rem += bLong - quot; quot -= 1; } else { rem += (bLong << 1) - quot; quot -= 2; } } } } return (rem << 32) | (quot & 0xffffffffL); } /** * Performs modular exponentiation using the Montgomery Reduction. It requires * that all parameters be positive and the modulus be even. Based The * square and multiply algorithm and the Montgomery Reduction C. K. Koc - * Montgomery Reduction with Even Modulus. The square and multiply * algorithm and the Montgomery Reduction. * * @ar.org.fitc.ref "C. K. Koc - Montgomery Reduction with Even Modulus" * @see BigInteger#modPow(BigInteger, BigInteger) */ static BigInteger evenModPow(BigInteger base, BigInteger exponent, BigInteger modulus) { // PRE: (base > 0), (exponent > 0), (modulus > 0) and (modulus even) // STEP 1: Obtain the factorization 'modulus'= q * 2^j. int j = modulus.getLowestSetBit(); BigInteger q = modulus.shiftRight(j); // STEP 2: Compute x1 := base^exponent (mod q). BigInteger x1 = oddModPow(base, exponent, q); // STEP 3: Compute x2 := base^exponent (mod 2^j). BigInteger x2 = pow2ModPow(base, exponent, j); // STEP 4: Compute q^(-1) (mod 2^j) and y := (x2-x1) * q^(-1) (mod 2^j) BigInteger qInv = modPow2Inverse(q, j); BigInteger y = (x2.subtract(x1)).multiply(qInv); inplaceModPow2(y, j); if (y.sign < 0) { y = y.add(BigInteger.getPowerOfTwo(j)); } // STEP 5: Compute and return: x1 + q * y return x1.add(q.multiply(y)); } /** * Performs the final reduction of the Montgomery algorithm. * * @see #monPro(BigInteger, BigInteger, BigInteger, long) * @see #monSquare(BigInteger, BigInteger, long) */ static BigInteger finalSubtraction(int res[], BigInteger modulus) { // skipping leading zeros int modulusLen = modulus.numberLength; boolean doSub = res[modulusLen] != 0; if (!doSub) { int modulusDigits[] = modulus.digits; doSub = true; for (int i = modulusLen - 1; i >= 0; i--) { if (res[i] != modulusDigits[i]) { doSub = (res[i] != 0) && ((res[i] & 0xFFFFFFFFL) > (modulusDigits[i] & 0xFFFFFFFFL)); break; } } } BigInteger result = new BigInteger(1, modulusLen + 1, res); // if (res >= modulusDigits) compute (res - modulusDigits) if (doSub) { Elementary.inplaceSubtract(result, modulus); } result.cutOffLeadingZeroes(); return result; } /** * @param m a positive modulus Return the greatest common divisor of op1 and * op2, * * @param op1 must be greater than zero * @param op2 must be greater than zero * @see BigInteger#gcd(BigInteger) * @return {@code GCD(op1, op2)} */ static BigInteger gcdBinary(BigInteger op1, BigInteger op2) { // PRE: (op1 > 0) and (op2 > 0) /* * Divide both number the maximal possible times by 2 without rounding * gcd(2*a, 2*b) = 2 * gcd(a,b) */ int lsb1 = op1.getLowestSetBit(); int lsb2 = op2.getLowestSetBit(); int pow2Count = Math.min(lsb1, lsb2); BitLevel.inplaceShiftRight(op1, lsb1); BitLevel.inplaceShiftRight(op2, lsb2); BigInteger swap; // I want op2 > op1 if (op1.compareTo(op2) == BigInteger.GREATER) { swap = op1; op1 = op2; op2 = swap; } do { // INV: op2 >= op1 && both are odd unless op1 = 0 // Optimization for small operands // (op2.bitLength() < 64) implies by INV (op1.bitLength() < 64) if ((op2.numberLength == 1) || ((op2.numberLength == 2) && (op2.digits[1] > 0))) { op2 = BigInteger.valueOf(Division.gcdBinary(op1.longValue(), op2.longValue())); break; } // Implements one step of the Euclidean algorithm // To reduce one operand if it's much smaller than the other one if (op2.numberLength > op1.numberLength * 1.2) { op2 = op2.remainder(op1); if (op2.signum() != 0) { BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit()); } } else { // Use Knuth's algorithm of successive subtract and shifting do { Elementary.inplaceSubtract(op2, op1); // both are odd BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit()); // op2 is even } while (op2.compareTo(op1) >= BigInteger.EQUALS); } // now op1 >= op2 swap = op2; op2 = op1; op1 = swap; } while (op1.sign != 0); return op2.shiftLeft(pow2Count); } /** * Performs the same as {@link #gcdBinary(BigInteger, BigInteger)}, but with * numbers of 63 bits, represented in positives values of {@code long} type. * * @param op1 a positive number * @param op2 a positive number * @see #gcdBinary(BigInteger, BigInteger) * @return GCD(op1, op2) */ static long gcdBinary(long op1, long op2) { // PRE: (op1 > 0) and (op2 > 0) int lsb1 = Long.numberOfTrailingZeros(op1); int lsb2 = Long.numberOfTrailingZeros(op2); int pow2Count = Math.min(lsb1, lsb2); if (lsb1 != 0) { op1 >>>= lsb1; } if (lsb2 != 0) { op2 >>>= lsb2; } do { if (op1 >= op2) { op1 -= op2; op1 >>>= Long.numberOfTrailingZeros(op1); } else { op2 -= op1; op2 >>>= Long.numberOfTrailingZeros(op2); } } while (op1 != 0); return (op2 << pow2Count); } /** * Performs {@code x = x mod (2n)}. * * @param x a positive number, it will store the result. * @param n a positive exponent of {@code 2}. */ static void inplaceModPow2(BigInteger x, int n) { // PRE: (x > 0) and (n >= 0) int fd = n >> 5; int leadingZeros; if ((x.numberLength < fd) || (x.bitLength() <= n)) { return; } leadingZeros = 32 - (n & 31); x.numberLength = fd + 1; x.digits[fd] &= (leadingZeros < 32) ? (-1 >>> leadingZeros) : 0; x.cutOffLeadingZeroes(); } /** * * Based on "New Algorithm for Classical Modular Inverse" Róbert Lórencz. LNCS * 2523 (2002) * * @return a^(-1) mod m */ static BigInteger modInverseLorencz(BigInteger a, BigInteger modulo) { // PRE: a is coprime with modulo, a < modulo int max = Math.max(a.numberLength, modulo.numberLength); int uDigits[] = new int[max + 1]; // enough place to make all the inplace // operation int vDigits[] = new int[max + 1]; System.arraycopy(modulo.digits, 0, uDigits, 0, modulo.numberLength); System.arraycopy(a.digits, 0, vDigits, 0, a.numberLength); BigInteger u = new BigInteger(modulo.sign, modulo.numberLength, uDigits); BigInteger v = new BigInteger(a.sign, a.numberLength, vDigits); BigInteger r = new BigInteger(0, 1, new int[max + 1]); // BigInteger.ZERO; BigInteger s = new BigInteger(1, 1, new int[max + 1]); s.digits[0] = 1; // r == 0 && s == 1, but with enough place int coefU = 0, coefV = 0; int n = modulo.bitLength(); int k; while (!isPowerOfTwo(u, coefU) && !isPowerOfTwo(v, coefV)) { // modification of original algorithm: I calculate how many times the // algorithm will enter in the same branch of if k = howManyIterations(u, n); if (k != 0) { BitLevel.inplaceShiftLeft(u, k); if (coefU >= coefV) { BitLevel.inplaceShiftLeft(r, k); } else { BitLevel.inplaceShiftRight(s, Math.min(coefV - coefU, k)); if (k - (coefV - coefU) > 0) { BitLevel.inplaceShiftLeft(r, k - coefV + coefU); } } coefU += k; } k = howManyIterations(v, n); if (k != 0) { BitLevel.inplaceShiftLeft(v, k); if (coefV >= coefU) { BitLevel.inplaceShiftLeft(s, k); } else { BitLevel.inplaceShiftRight(r, Math.min(coefU - coefV, k)); if (k - (coefU - coefV) > 0) { BitLevel.inplaceShiftLeft(s, k - coefU + coefV); } } coefV += k; } if (u.signum() == v.signum()) { if (coefU <= coefV) { Elementary.completeInPlaceSubtract(u, v); Elementary.completeInPlaceSubtract(r, s); } else { Elementary.completeInPlaceSubtract(v, u); Elementary.completeInPlaceSubtract(s, r); } } else { if (coefU <= coefV) { Elementary.completeInPlaceAdd(u, v); Elementary.completeInPlaceAdd(r, s); } else { Elementary.completeInPlaceAdd(v, u); Elementary.completeInPlaceAdd(s, r); } } if (v.signum() == 0 || u.signum() == 0) { // math.19: BigInteger not invertible throw new ArithmeticException("BigInteger not invertible."); } } if (isPowerOfTwo(v, coefV)) { r = s; if (v.signum() != u.signum()) { u = u.negate(); } } if (u.testBit(n)) { if (r.signum() < 0) { r = r.negate(); } else { r = modulo.subtract(r); } } if (r.signum() < 0) { r = r.add(modulo); } return r; } /** * Calculates a.modInverse(p) Based on: Savas, E; Koc, C "The Montgomery * Modular Inverse - Revised". */ static BigInteger modInverseMontgomery(BigInteger a, BigInteger p) { if (a.sign == 0) { // ZERO hasn't inverse // math.19: BigInteger not invertible throw new ArithmeticException("BigInteger not invertible."); } if (!p.testBit(0)) { // montgomery inverse require even modulo return modInverseLorencz(a, p); } int m = p.numberLength * 32; // PRE: a \in [1, p - 1] BigInteger u, v, r, s; u = p.copy(); // make copy to use inplace method v = a.copy(); int max = Math.max(v.numberLength, u.numberLength); r = new BigInteger(1, 1, new int[max + 1]); s = new BigInteger(1, 1, new int[max + 1]); s.digits[0] = 1; // s == 1 && v == 0 int k = 0; int lsbu = u.getLowestSetBit(); int lsbv = v.getLowestSetBit(); int toShift; if (lsbu > lsbv) { BitLevel.inplaceShiftRight(u, lsbu); BitLevel.inplaceShiftRight(v, lsbv); BitLevel.inplaceShiftLeft(r, lsbv); k += lsbu - lsbv; } else { BitLevel.inplaceShiftRight(u, lsbu); BitLevel.inplaceShiftRight(v, lsbv); BitLevel.inplaceShiftLeft(s, lsbu); k += lsbv - lsbu; } r.sign = 1; while (v.signum() > 0) { // INV v >= 0, u >= 0, v odd, u odd (except last iteration when v is even // (0)) while (u.compareTo(v) > BigInteger.EQUALS) { Elementary.inplaceSubtract(u, v); toShift = u.getLowestSetBit(); BitLevel.inplaceShiftRight(u, toShift); Elementary.inplaceAdd(r, s); BitLevel.inplaceShiftLeft(s, toShift); k += toShift; } while (u.compareTo(v) <= BigInteger.EQUALS) { Elementary.inplaceSubtract(v, u); if (v.signum() == 0) { break; } toShift = v.getLowestSetBit(); BitLevel.inplaceShiftRight(v, toShift); Elementary.inplaceAdd(s, r); BitLevel.inplaceShiftLeft(r, toShift); k += toShift; } } if (!u.isOne()) { // in u is stored the gcd // math.19: BigInteger not invertible. throw new ArithmeticException("BigInteger not invertible."); } if (r.compareTo(p) >= BigInteger.EQUALS) { Elementary.inplaceSubtract(r, p); } r = p.subtract(r); // Have pair: ((BigInteger)r, (Integer)k) where r == a^(-1) * 2^k mod // (module) int n1 = calcN(p); if (k > m) { r = monPro(r, BigInteger.ONE, p, n1); k = k - m; } r = monPro(r, BigInteger.getPowerOfTwo(m - k), p, n1); return r; } /** * @param x an odd positive number. * @param n the exponent by which 2 is raised. * @return {@code x-1 (mod 2n)}. */ static BigInteger modPow2Inverse(BigInteger x, int n) { // PRE: (x > 0), (x is odd), and (n > 0) BigInteger y = new BigInteger(1, new int[1 << n]); y.numberLength = 1; y.digits[0] = 1; y.sign = 1; for (int i = 1; i < n; i++) { if (BitLevel.testBit(x.multiply(y), i)) { // Adding 2^i to y (setting the i-th bit) y.digits[i >> 5] |= (1 << (i & 31)); } } return y; } /** * Implements the Montgomery Product of two integers represented by {@code * int} arrays. The arrays are supposed in little endian notation. * * @param a The first factor of the product. * @param b The second factor of the product. * @param modulus The modulus of the operations. Zmodulus. * @param n2 The digit modulus'[0]. * @ar.org.fitc.ref "C. K. Koc - Analyzing and Comparing Montgomery * Multiplication Algorithms" * @see #modPowOdd(BigInteger, BigInteger, BigInteger) */ static BigInteger monPro(BigInteger a, BigInteger b, BigInteger modulus, int n2) { int modulusLen = modulus.numberLength; int res[] = new int[(modulusLen << 1) + 1]; Multiplication.multArraysPAP(a.digits, Math.min(modulusLen, a.numberLength), b.digits, Math.min(modulusLen, b.numberLength), res); monReduction(res, modulus, n2); return finalSubtraction(res, modulus); } /** * Multiplies an array by int and subtracts it from a subarray of another * array. * * @param a the array to subtract from * @param start the start element of the subarray of a * @param b the array to be multiplied and subtracted * @param bLen the length of b * @param c the multiplier of b * @return the carry element of subtraction */ static int multiplyAndSubtract(int a[], int start, int b[], int bLen, int c) { long carry0 = 0; long carry1 = 0; for (int i = 0; i < bLen; i++) { carry0 = Multiplication.unsignedMultAddAdd(b[i], c, (int) carry0, 0); carry1 = (a[start + i] & 0xffffffffL) - (carry0 & 0xffffffffL) + carry1; a[start + i] = (int) carry1; carry1 >>= 32; // -1 or 0 carry0 >>>= 32; } carry1 = (a[start + bLen] & 0xffffffffL) - carry0 + carry1; a[start + bLen] = (int) carry1; return (int) (carry1 >> 32); // -1 or 0 } /** * Performs modular exponentiation using the Montgomery Reduction. It requires * that all parameters be positive and the modulus be odd. > * * @see BigInteger#modPow(BigInteger, BigInteger) * @see #monPro(BigInteger, BigInteger, BigInteger, int) * @see #slidingWindow(BigInteger, BigInteger, BigInteger, BigInteger, int) * @see #squareAndMultiply(BigInteger, BigInteger, BigInteger, BigInteger, * int) */ static BigInteger oddModPow(BigInteger base, BigInteger exponent, BigInteger modulus) { // PRE: (base > 0), (exponent > 0), (modulus > 0) and (odd modulus) int k = (modulus.numberLength << 5); // r = 2^k // n-residue of base [base * r (mod modulus)] BigInteger a2 = base.shiftLeft(k).mod(modulus); // n-residue of base [1 * r (mod modulus)] BigInteger x2 = BigInteger.getPowerOfTwo(k).mod(modulus); BigInteger res; // Compute (modulus[0]^(-1)) (mod 2^32) for odd modulus int n2 = calcN(modulus); if (modulus.numberLength == 1) { res = squareAndMultiply(x2, a2, exponent, modulus, n2); } else { res = slidingWindow(x2, a2, exponent, modulus, n2); } return monPro(res, BigInteger.ONE, modulus, n2); } /** * It requires that all parameters be positive. * * @return {@code baseexponent mod (2j)}. * @see BigInteger#modPow(BigInteger, BigInteger) */ static BigInteger pow2ModPow(BigInteger base, BigInteger exponent, int j) { // PRE: (base > 0), (exponent > 0) and (j > 0) BigInteger res = BigInteger.ONE; BigInteger e = exponent.copy(); BigInteger baseMod2toN = base.copy(); BigInteger res2; /* * If 'base' is odd then it's coprime with 2^j and phi(2^j) = 2^(j-1); so we * can reduce reduce the exponent (mod 2^(j-1)). */ if (base.testBit(0)) { inplaceModPow2(e, j - 1); } inplaceModPow2(baseMod2toN, j); for (int i = e.bitLength() - 1; i >= 0; i--) { res2 = res.copy(); inplaceModPow2(res2, j); res = res.multiply(res2); if (BitLevel.testBit(e, i)) { res = res.multiply(baseMod2toN); inplaceModPow2(res, j); } } inplaceModPow2(res, j); return res; } /** * Divides a BigInteger by a signed int and returns * the remainder. * * @param dividend the BigInteger to be divided. Must be non-negative. * @param divisor a signed int * @return divide % divisor */ static int remainder(BigInteger dividend, int divisor) { return remainderArrayByInt(dividend.digits, dividend.numberLength, divisor); } /** * Divides an array by an integer value. Implements the Knuth's division * algorithm. See D. Knuth, The Art of Computer Programming, vol. 2. * * @param src the dividend * @param srcLength the length of the dividend * @param divisor the divisor * @return remainder */ static int remainderArrayByInt(int src[], final int srcLength, final int divisor) { long result = 0; for (int i = srcLength - 1; i >= 0; i--) { long temp = (result << 32) + (src[i] & 0xffffffffL); long res = divideLongByInt(temp, divisor); result = (int) (res >> 32); } return (int) result; } /* * Implements the Montgomery modular exponentiation based in The sliding * windows algorithm and the MongomeryReduction. * * @ar.org.fitc.ref * "A. Menezes,P. van Oorschot, S. Vanstone - Handbook of Applied Cryptography" * ; * * @see #oddModPow(BigInteger, BigInteger, BigInteger) */ static BigInteger slidingWindow(BigInteger x2, BigInteger a2, BigInteger exponent, BigInteger modulus, int n2) { // fill odd low pows of a2 BigInteger pows[] = new BigInteger[8]; BigInteger res = x2; int lowexp; BigInteger x3; int acc3; pows[0] = a2; x3 = monPro(a2, a2, modulus, n2); for (int i = 1; i <= 7; i++) { pows[i] = monPro(pows[i - 1], x3, modulus, n2); } for (int i = exponent.bitLength() - 1; i >= 0; i--) { if (BitLevel.testBit(exponent, i)) { lowexp = 1; acc3 = i; for (int j = Math.max(i - 3, 0); j <= i - 1; j++) { if (BitLevel.testBit(exponent, j)) { if (j < acc3) { acc3 = j; lowexp = (lowexp << (i - j)) ^ 1; } else { lowexp = lowexp ^ (1 << (j - acc3)); } } } for (int j = acc3; j <= i; j++) { res = monPro(res, res, modulus, n2); } res = monPro(pows[(lowexp - 1) >> 1], res, modulus, n2); i = acc3; } else { res = monPro(res, res, modulus, n2); } } return res; } static BigInteger squareAndMultiply(BigInteger x2, BigInteger a2, BigInteger exponent, BigInteger modulus, int n2) { BigInteger res = x2; for (int i = exponent.bitLength() - 1; i >= 0; i--) { res = monPro(res, res, modulus, n2); if (BitLevel.testBit(exponent, i)) { res = monPro(res, a2, modulus, n2); } } return res; } /** * Calculate the first digit of the inverse. */ private static int calcN(BigInteger a) { long m0 = a.digits[0] & 0xFFFFFFFFL; long n2 = 1L; // this is a'[0] long powerOfTwo = 2L; do { if (((m0 * n2) & powerOfTwo) != 0) { n2 |= powerOfTwo; } powerOfTwo <<= 1; } while (powerOfTwo < 0x100000000L); n2 = -n2; return (int) (n2 & 0xFFFFFFFFL); } /** * Calculate how many iteration of Lorencz's algorithm would perform the same * operation. * * @param bi * @param n * @return */ private static int howManyIterations(BigInteger bi, int n) { int i = n - 1; if (bi.sign > 0) { while (!bi.testBit(i)) { i--; } return n - 1 - i; } else { while (bi.testBit(i)) { i--; } return n - 1 - Math.max(i, bi.getLowestSetBit()); } } /** * Returns {@code bi == abs(2^exp)}. */ private static boolean isPowerOfTwo(BigInteger bi, int exp) { boolean result = false; result = (exp >> 5 == bi.numberLength - 1) && (bi.digits[bi.numberLength - 1] == 1 << (exp & 31)); if (result) { for (int i = 0; result && i < bi.numberLength - 1; i++) { result = bi.digits[i] == 0; } } return result; } private static void monReduction(int[] res, BigInteger modulus, int n2) { /* res + m*modulus_digits */ int[] modulusDigits = modulus.digits; int modulusLen = modulus.numberLength; long outerCarry = 0; for (int i = 0; i < modulusLen; i++) { long innnerCarry = 0; int m = (int) Multiplication.unsignedMultAddAdd(res[i], n2, 0, 0); for (int j = 0; j < modulusLen; j++) { innnerCarry = Multiplication.unsignedMultAddAdd(m, modulusDigits[j], res[i + j], (int) innnerCarry); res[i + j] = (int) innnerCarry; innnerCarry >>>= 32; } outerCarry += (res[i + modulusLen] & 0xFFFFFFFFL) + innnerCarry; res[i + modulusLen] = (int) outerCarry; outerCarry >>>= 32; } res[modulusLen << 1] = (int) outerCarry; /* res / r */ for (int j = 0; j < modulusLen + 1; j++) { res[j] = res[j + modulusLen]; } } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy