| /* |
| * 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: |
| * |
| * <ul type="circle"> <li><b>Division</b> <ul type="circle"> <li> |
| * {@link BigInteger} division and remainder by {@link BigInteger}.</li> <li> |
| * {@link BigInteger} division and remainder by {@code int}.</li> <li><i>gcd</i> |
| * between {@link BigInteger} numbers.</li> </ul> </li> <li><b>Modular |
| * arithmetic </b> <ul type="circle"> <li>Modular exponentiation between |
| * {@link BigInteger} numbers.</li> <li>Modular inverse of a {@link BigInteger} |
| * numbers.</li> </ul> </li> </ul> |
| */ |
| 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 <i>The |
| * square and multiply algorithm and the Montgomery Reduction C. K. Koc - |
| * Montgomery Reduction with Even Modulus</i>. 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 <code>GCD(op1, op2)</code> |
| */ |
| 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 (2<sup>n</sup>)}. |
| * |
| * @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<sup>-1</sup> (mod 2<sup>n</sup>)}. |
| */ |
| 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 <i>little endian</i> notation. |
| * |
| * @param a The first factor of the product. |
| * @param b The second factor of the product. |
| * @param modulus The modulus of the operations. Z<sub>modulus</sub>. |
| * @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 base<sup>exponent</sup> mod (2<sup>j</sup>)}. |
| * @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 <code>BigInteger</code> by a signed <code>int</code> 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 <i>The sliding |
| * windows algorithm and the MongomeryReduction</i>. |
| * |
| * @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()); |
| } |
| } |
| |
| /** |
| * @return 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]; |
| } |
| } |
| } |
