0001: /*
0002: * Licensed to the Apache Software Foundation (ASF) under one or more
0003: * contributor license agreements. See the NOTICE file distributed with
0004: * this work for additional information regarding copyright ownership.
0005: * The ASF licenses this file to You under the Apache License, Version 2.0
0006: * (the "License"); you may not use this file except in compliance with
0007: * the License. You may obtain a copy of the License at
0008: *
0009: * http://www.apache.org/licenses/LICENSE-2.0
0010: *
0011: * Unless required by applicable law or agreed to in writing, software
0012: * distributed under the License is distributed on an "AS IS" BASIS,
0013: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0014: * See the License for the specific language governing permissions and
0015: * limitations under the License.
0016: */
0017:
0018: package java.math;
0019:
0020: import org.apache.harmony.math.internal.nls.Messages;
0021:
0022: /**
0023: * Static library that provides all operations related with division and modular
0024: * arithmetic to {@link BigInteger}. Some methods are provided in both mutable
0025: * and immutable way. There are several variants provided listed below:
0026: *
0027: * <ul type="circle">
0028: * <li> <b>Division</b>
0029: * <ul type="circle">
0030: * <li>{@link BigInteger} division and remainder by {@link BigInteger}.</li>
0031: * <li>{@link BigInteger} division and remainder by {@code int}.</li>
0032: * <li><i>gcd</i> between {@link BigInteger} numbers.</li>
0033: * </ul>
0034: * </li>
0035: * <li> <b>Modular arithmetic </b>
0036: * <ul type="circle">
0037: * <li>Modular exponentiation between {@link BigInteger} numbers.</li>
0038: * <li>Modular inverse of a {@link BigInteger} numbers.</li>
0039: * </ul>
0040: * </li>
0041: *</ul>
0042: *
0043: * @author Intel Middleware Product Division
0044: * @author Instituto Tecnologico de Cordoba
0045: */
0046: class Division {
0047:
0048: /**
0049: * Divides the array 'a' by the array 'b' and gets the quotient and the
0050: * remainder. Implements the Knuth's division algorithm. See D. Knuth, The
0051: * Art of Computer Programming, vol. 2. Steps D1-D8 correspond the steps in
0052: * the algorithm description.
0053: *
0054: * @param quot the quotient
0055: * @param quotLength the quotient's length
0056: * @param a the dividend
0057: * @param aLength the dividend's length
0058: * @param b the divisor
0059: * @param bLength the divisor's length
0060: * @return the remainder
0061: */
0062: static int[] divide(int quot[], int quotLength, int a[],
0063: int aLength, int b[], int bLength) {
0064:
0065: int normA[] = new int[aLength + 1]; // the normalized dividend
0066: // an extra byte is needed for correct shift
0067: int normB[] = new int[bLength + 1]; // the normalized divisor;
0068: int normBLength = bLength;
0069: /*
0070: * Step D1: normalize a and b and put the results to a1 and b1 the
0071: * normalized divisor's first digit must be >= 2^31
0072: */
0073: int divisorShift = Integer.numberOfLeadingZeros(b[bLength - 1]);
0074: if (divisorShift != 0) {
0075: BitLevel.shiftLeft(normB, b, 0, divisorShift);
0076: BitLevel.shiftLeft(normA, a, 0, divisorShift);
0077: } else {
0078: System.arraycopy(a, 0, normA, 0, aLength);
0079: System.arraycopy(b, 0, normB, 0, bLength);
0080: }
0081: int firstDivisorDigit = normB[normBLength - 1];
0082: // Step D2: set the quotient index
0083: int i = quotLength - 1;
0084: int j = aLength;
0085:
0086: while (i >= 0) {
0087: // Step D3: calculate a guess digit guessDigit
0088: int guessDigit = 0;
0089: if (normA[j] == firstDivisorDigit) {
0090: // set guessDigit to the largest unsigned int value
0091: guessDigit = -1;
0092: } else {
0093: long product = (((normA[j] & 0xffffffffL) << 32) + (normA[j - 1] & 0xffffffffL));
0094: long res = Division.divideLongByInt(product,
0095: firstDivisorDigit);
0096: guessDigit = (int) res; // the quotient of divideLongByInt
0097: int rem = (int) (res >> 32); // the remainder of
0098: // divideLongByInt
0099: // decrease guessDigit by 1 while leftHand > rightHand
0100: if (guessDigit != 0) {
0101: long leftHand = 0;
0102: long rightHand = 0;
0103: boolean rOverflowed = false;
0104: guessDigit++; // to have the proper value in the loop
0105: // below
0106: do {
0107: guessDigit--;
0108: if (rOverflowed) {
0109: break;
0110: }
0111: // leftHand always fits in an unsigned long
0112: leftHand = (guessDigit & 0xffffffffL)
0113: * (normB[normBLength - 2] & 0xffffffffL);
0114: /*
0115: * rightHand can overflow; in this case the loop
0116: * condition will be true in the next step of the loop
0117: */
0118: rightHand = ((long) rem << 32)
0119: + (normA[j - 2] & 0xffffffffL);
0120: long longR = (rem & 0xffffffffL)
0121: + (firstDivisorDigit & 0xffffffffL);
0122: /*
0123: * checks that longR does not fit in an unsigned int;
0124: * this ensures that rightHand will overflow unsigned
0125: * long in the next step
0126: */
0127: if (Integer
0128: .numberOfLeadingZeros((int) (longR >>> 32)) < 32) {
0129: rOverflowed = true;
0130: } else {
0131: rem = (int) longR;
0132: }
0133: } while (((leftHand ^ 0x8000000000000000L) > (rightHand ^ 0x8000000000000000L)));
0134: }
0135: }
0136: // Step D4: multiply normB by guessDigit and subtract the production
0137: // from normA.
0138: if (guessDigit != 0) {
0139: int borrow = Division.multiplyAndSubtract(normA, j
0140: - normBLength, normB, normBLength,
0141: guessDigit & 0xffffffffL);
0142: // Step D5: check the borrow
0143: if (borrow != 0) {
0144: // Step D6: compensating addition
0145: guessDigit--;
0146: long carry = 0;
0147: for (int k = 0; k < normBLength; k++) {
0148: carry += (normA[j - normBLength + k] & 0xffffffffL)
0149: + (normB[k] & 0xffffffffL);
0150: normA[j - normBLength + k] = (int) carry;
0151: carry >>>= 32;
0152: }
0153: }
0154: }
0155: if (quot != null) {
0156: quot[i] = guessDigit;
0157: }
0158: // Step D7
0159: j--;
0160: i--;
0161: }
0162: /*
0163: * Step D8: we got the remainder in normA. Denormalize it id needed
0164: */
0165: if (divisorShift != 0) {
0166: // reuse normB
0167: BitLevel.shiftRight(normB, normBLength, normA, 0,
0168: divisorShift);
0169: return normB;
0170: }
0171: System.arraycopy(normA, 0, normB, 0, bLength);
0172: return normA;
0173: }
0174:
0175: /**
0176: * Divides an array by an integer value. Implements the Knuth's division
0177: * algorithm. See D. Knuth, The Art of Computer Programming, vol. 2.
0178: *
0179: * @param dest the quotient
0180: * @param src the dividend
0181: * @param srcLength the length of the dividend
0182: * @param divisor the divisor
0183: * @return remainder
0184: */
0185: static int divideArrayByInt(int dest[], int src[],
0186: final int srcLength, final int divisor) {
0187:
0188: long rem = 0;
0189: long bLong = divisor & 0xffffffffL;
0190:
0191: for (int i = srcLength - 1; i >= 0; i--) {
0192: long temp = (rem << 32) | (src[i] & 0xffffffffL);
0193: long quot;
0194: if (temp >= 0) {
0195: quot = (temp / bLong);
0196: rem = (temp % bLong);
0197: } else {
0198: /*
0199: * make the dividend positive shifting it right by 1 bit then
0200: * get the quotient an remainder and correct them properly
0201: */
0202: long aPos = temp >>> 1;
0203: long bPos = divisor >>> 1;
0204: quot = aPos / bPos;
0205: rem = aPos % bPos;
0206: // double the remainder and add 1 if a is odd
0207: rem = (rem << 1) + (temp & 1);
0208: if ((divisor & 1) != 0) {
0209: // the divisor is odd
0210: if (quot <= rem) {
0211: rem -= quot;
0212: } else {
0213: if (quot - rem <= bLong) {
0214: rem += bLong - quot;
0215: quot -= 1;
0216: } else {
0217: rem += (bLong << 1) - quot;
0218: quot -= 2;
0219: }
0220: }
0221: }
0222: }
0223: dest[i] = (int) (quot & 0xffffffffL);
0224: }
0225: return (int) rem;
0226: }
0227:
0228: /**
0229: * Divides an array by an integer value. Implements the Knuth's division
0230: * algorithm. See D. Knuth, The Art of Computer Programming, vol. 2.
0231: *
0232: * @param src the dividend
0233: * @param srcLength the length of the dividend
0234: * @param divisor the divisor
0235: * @return remainder
0236: */
0237: static int remainderArrayByInt(int src[], final int srcLength,
0238: final int divisor) {
0239:
0240: long result = 0;
0241:
0242: for (int i = srcLength - 1; i >= 0; i--) {
0243: long temp = (result << 32) + (src[i] & 0xffffffffL);
0244: long res = divideLongByInt(temp, divisor);
0245: result = (int) (res >> 32);
0246: }
0247: return (int) result;
0248: }
0249:
0250: /**
0251: * Divides a <code>BigInteger</code> by a signed <code>int</code> and
0252: * returns the remainder.
0253: *
0254: * @param dividend the BigInteger to be divided. Must be non-negative.
0255: * @param divisor a signed int
0256: * @return divide % divisor
0257: */
0258: static int remainder(BigInteger dividend, int divisor) {
0259: return remainderArrayByInt(dividend.digits,
0260: dividend.numberLength, divisor);
0261: }
0262:
0263: /**
0264: * Divides an unsigned long a by an unsigned int b. It is supposed that the
0265: * most significant bit of b is set to 1, i.e. b < 0
0266: *
0267: * @param a the dividend
0268: * @param b the divisor
0269: * @return the long value containing the unsigned integer remainder in the
0270: * left half and the unsigned integer quotient in the right half
0271: */
0272: static long divideLongByInt(long a, int b) {
0273: long quot;
0274: long rem;
0275: long bLong = b & 0xffffffffL;
0276:
0277: if (a >= 0) {
0278: quot = (a / bLong);
0279: rem = (a % bLong);
0280: } else {
0281: /*
0282: * Make the dividend positive shifting it right by 1 bit then get
0283: * the quotient an remainder and correct them properly
0284: */
0285: long aPos = a >>> 1;
0286: long bPos = b >>> 1;
0287: quot = aPos / bPos;
0288: rem = aPos % bPos;
0289: // double the remainder and add 1 if a is odd
0290: rem = (rem << 1) + (a & 1);
0291: if ((b & 1) != 0) { // the divisor is odd
0292: if (quot <= rem) {
0293: rem -= quot;
0294: } else {
0295: if (quot - rem <= bLong) {
0296: rem += bLong - quot;
0297: quot -= 1;
0298: } else {
0299: rem += (bLong << 1) - quot;
0300: quot -= 2;
0301: }
0302: }
0303: }
0304: }
0305: return (rem << 32) | (quot & 0xffffffffL);
0306: }
0307:
0308: /**
0309: * Computes the quotient and the remainder after a division by an {@code int}
0310: * number.
0311: *
0312: * @return an array of the form {@code [quotient, remainder]}.
0313: */
0314: static BigInteger[] divideAndRemainderByInteger(BigInteger val,
0315: int divisor, int divisorSign) {
0316: // res[0] is a quotient and res[1] is a remainder:
0317: int[] valDigits = val.digits;
0318: int valLen = val.numberLength;
0319: int valSign = val.sign;
0320: if (valLen == 1) {
0321: long a = (valDigits[0] & 0xffffffffL);
0322: long b = (divisor & 0xffffffffL);
0323: long quo = a / b;
0324: long rem = a % b;
0325: if (valSign != divisorSign) {
0326: quo = -quo;
0327: }
0328: if (valSign < 0) {
0329: rem = -rem;
0330: }
0331: return new BigInteger[] { BigInteger.valueOf(quo),
0332: BigInteger.valueOf(rem) };
0333: }
0334: int quotientLength = valLen;
0335: int quotientSign = ((valSign == divisorSign) ? 1 : -1);
0336: int quotientDigits[] = new int[quotientLength];
0337: int remainderDigits[];
0338: remainderDigits = new int[] { Division.divideArrayByInt(
0339: quotientDigits, valDigits, valLen, divisor) };
0340: BigInteger result0 = new BigInteger(quotientSign,
0341: quotientLength, quotientDigits);
0342: BigInteger result1 = new BigInteger(valSign, 1, remainderDigits);
0343: result0.cutOffLeadingZeroes();
0344: result1.cutOffLeadingZeroes();
0345: return new BigInteger[] { result0, result1 };
0346: }
0347:
0348: /**
0349: * Multiplies an array by int and subtracts it from a subarray of another
0350: * array.
0351: *
0352: * @param a the array to subtract from
0353: * @param start the start element of the subarray of a
0354: * @param b the array to be multiplied and subtracted
0355: * @param bLen the length of b
0356: * @param c the multiplier of b
0357: * @return the carry element of subtraction
0358: */
0359: static int multiplyAndSubtract(int a[], int start, int b[],
0360: int bLen, long c) {
0361: int i;
0362: int carry = 0;
0363: long product;
0364: int productInt;
0365:
0366: for (i = 0; i < bLen; i++) {
0367: product = c * (b[i] & 0xffffffffL);
0368: productInt = (int) product;
0369: productInt += carry;
0370: carry = (int) (product >> 32)
0371: + ((productInt ^ 0x80000000) < (carry ^ 0x80000000) ? 1
0372: : 0);
0373: productInt = a[start + i] - productInt;
0374: if ((productInt ^ 0x80000000) > (a[start + i] ^ 0x80000000)) {
0375: carry++;
0376: }
0377: a[start + i] = productInt;
0378: }
0379: product = (a[start + i] & 0xffffffffL) - (carry & 0xffffffffL);
0380: a[start + i] = (int) product;
0381: carry = (int) (product >> 32); // -1 or 0
0382: return carry;
0383: }
0384:
0385: /**
0386: * @param m a positive modulus
0387: * Return the greatest common divisor of op1 and op2,
0388: *
0389: * @param op1
0390: * must be greater than zero
0391: * @param op2
0392: * must be greater than zero
0393: * @see BigInteger#gcd(BigInteger)
0394: * @return {@code GCD(op1, op2)}
0395: */
0396: static BigInteger gcdBinary(BigInteger op1, BigInteger op2) {
0397: // PRE: (op1 > 0) and (op2 > 0)
0398:
0399: /*
0400: * Divide both number the maximal possible times by 2 without rounding
0401: * gcd(2*a, 2*b) = 2 * gcd(a,b)
0402: */
0403: int lsb1 = op1.getLowestSetBit();
0404: int lsb2 = op2.getLowestSetBit();
0405: int pow2Count = Math.min(lsb1, lsb2);
0406:
0407: BitLevel.inplaceShiftRight(op1, lsb1);
0408: BitLevel.inplaceShiftRight(op2, lsb2);
0409:
0410: BigInteger swap;
0411: // I want op2 > op1
0412: if (op1.compareTo(op2) == BigInteger.GREATER) {
0413: swap = op1;
0414: op1 = op2;
0415: op2 = swap;
0416: }
0417:
0418: do { // INV: op2 >= op1 && both are odd unless op1 = 0
0419:
0420: // Optimization for small operands
0421: // (op2.bitLength() < 64) implies by INV (op1.bitLength() < 64)
0422: if ((op2.numberLength == 1)
0423: || ((op2.numberLength == 2) && (op2.digits[1] > 0))) {
0424: op2 = BigInteger.valueOf(Division.gcdBinary(op1
0425: .longValue(), op2.longValue()));
0426: break;
0427: }
0428:
0429: // Implements one step of the Euclidean algorithm
0430: // To reduce one operand if it's much smaller than the other one
0431: if (op2.numberLength > op1.numberLength * 1.2) {
0432: op2 = op2.remainder(op1);
0433: if (op2.signum() != 0) {
0434: BitLevel.inplaceShiftRight(op2, op2
0435: .getLowestSetBit());
0436: }
0437: } else {
0438:
0439: // Use Knuth's algorithm of successive subtract and shifting
0440: do {
0441: Elementary.inplaceSubtract(op2, op1); // both are odd
0442: BitLevel.inplaceShiftRight(op2, op2
0443: .getLowestSetBit()); // op2 is even
0444: } while (op2.compareTo(op1) >= BigInteger.EQUALS);
0445: }
0446: // now op1 >= op2
0447: swap = op2;
0448: op2 = op1;
0449: op1 = swap;
0450: } while (op1.sign != 0);
0451: return op2.shiftLeft(pow2Count);
0452: }
0453:
0454: /**
0455: * Performs the same as {@link #gcdBinary(BigInteger, BigInteger)}, but
0456: * with numbers of 63 bits, represented in positives values of {@code long}
0457: * type.
0458: *
0459: * @param op1
0460: * a positive number
0461: * @param op2
0462: * a positive number
0463: * @see #gcdBinary(BigInteger, BigInteger)
0464: * @return <code>GCD(op1, op2)</code>
0465: */
0466: static long gcdBinary(long op1, long op2) {
0467: // PRE: (op1 > 0) and (op2 > 0)
0468: int lsb1 = Long.numberOfTrailingZeros(op1);
0469: int lsb2 = Long.numberOfTrailingZeros(op2);
0470: int pow2Count = Math.min(lsb1, lsb2);
0471:
0472: if (lsb1 != 0) {
0473: op1 >>>= lsb1;
0474: }
0475: if (lsb2 != 0) {
0476: op2 >>>= lsb2;
0477: }
0478: do {
0479: if (op1 >= op2) {
0480: op1 -= op2;
0481: op1 >>>= Long.numberOfTrailingZeros(op1);
0482: } else {
0483: op2 -= op1;
0484: op2 >>>= Long.numberOfTrailingZeros(op2);
0485: }
0486: } while (op1 != 0);
0487: return (op2 << pow2Count);
0488: }
0489:
0490: /**
0491: * Calculates a.modInverse(p) Based on: Savas, E; Koc, C "The Montgomery Modular
0492: * Inverse - Revised"
0493: */
0494: static BigInteger modInverseMontgomery(BigInteger a, BigInteger p) {
0495:
0496: if (a.sign == 0) {
0497: // ZERO hasn't inverse
0498: // math.19: BigInteger not invertible
0499: throw new ArithmeticException(Messages.getString("math.19"));
0500: }
0501:
0502: if (!p.testBit(0)) {
0503: // montgomery inverse require even modulo
0504: return modInverseLorencz(a, p);
0505: }
0506:
0507: int m = p.numberLength * 32;
0508:
0509: Object[] save = almostMonInv(a, p);
0510: BigInteger r = ((BigInteger) save[0]);
0511: int k = ((Integer) save[1]).intValue();
0512: // assert ( k >= n && k <= m + n );
0513:
0514: long n1 = calcN(p);
0515:
0516: if (k > m) {
0517: r = monPro(r, BigInteger.ONE, p, n1);
0518: k = k - m;
0519: // assert k < m;
0520: }
0521: r = monPro(r, BigInteger.ONE.shiftLeft(m - k), p, n1);
0522: return r;
0523: }
0524:
0525: /**
0526: * Calculate the first digit of the inverse
0527: */
0528: private static long calcN(BigInteger a) {
0529: long m0 = a.digits[0] & 0xFFFFFFFFL;
0530: long n2 = 1L; // this is a'[0]
0531: long powerOfTwo = 2L;
0532:
0533: do {
0534: if (((m0 * n2) & powerOfTwo) != 0) {
0535: n2 |= powerOfTwo;
0536: }
0537: powerOfTwo <<= 1;
0538: } while (powerOfTwo < 0x100000000L);
0539: n2 = -n2;
0540: return n2;
0541: }
0542:
0543: /**
0544: * Used for an intermediate result of the modInverse algorithm
0545: * @return the pair: ((BigInteger)r, (Integer)k) where r == a^(-1) * 2^k mod (module)
0546: */
0547: private static Object[] almostMonInv(BigInteger a, BigInteger module) {
0548: // PRE: a \in [1, p - 1]
0549: BigInteger u, v, r, s;
0550: // make copy to use inplace method
0551: u = module.copy();
0552: v = a.copy();
0553:
0554: int max = Math.max(v.numberLength, u.numberLength);
0555: r = new BigInteger(1, 1, new int[max + 1]);
0556: s = new BigInteger(1, 1, new int[max + 1]);
0557: s.digits[0] = 1;
0558: // s == 1 && v == 0
0559:
0560: int k = 0;
0561:
0562: int lsbu = u.getLowestSetBit();
0563: int lsbv = v.getLowestSetBit();
0564: int toShift;
0565:
0566: if (lsbu > lsbv) {
0567: BitLevel.inplaceShiftRight(u, lsbu);
0568: BitLevel.inplaceShiftRight(v, lsbv);
0569: BitLevel.inplaceShiftLeft(r, lsbv);
0570: k += lsbu - lsbv;
0571: } else {
0572: BitLevel.inplaceShiftRight(u, lsbu);
0573: BitLevel.inplaceShiftRight(v, lsbv);
0574: BitLevel.inplaceShiftLeft(s, lsbu);
0575: k += lsbv - lsbu;
0576:
0577: }
0578:
0579: r.sign = 1;
0580: while (v.signum() > 0) {
0581: // INV v >= 0, u >= 0, v odd, u odd (except last iteration when v is even (0))
0582:
0583: while (u.compareTo(v) > BigInteger.EQUALS) {
0584: Elementary.inplaceSubtract(u, v);
0585: toShift = u.getLowestSetBit();
0586: BitLevel.inplaceShiftRight(u, toShift);
0587: Elementary.inplaceAdd(r, s);
0588: BitLevel.inplaceShiftLeft(s, toShift);
0589: k += toShift;
0590:
0591: }
0592:
0593: while (u.compareTo(v) <= BigInteger.EQUALS) {
0594: Elementary.inplaceSubtract(v, u);
0595: if (v.signum() == 0)
0596: break;
0597: toShift = v.getLowestSetBit();
0598: BitLevel.inplaceShiftRight(v, toShift);
0599: Elementary.inplaceAdd(s, r);
0600: BitLevel.inplaceShiftLeft(r, toShift);
0601: k += toShift;
0602:
0603: }
0604:
0605: }
0606: if (!u.isOne()) {
0607: // in u is stored the gcd
0608: // math.19: BigInteger not invertible.
0609: throw new ArithmeticException(Messages.getString("math.19"));
0610: }
0611:
0612: if (r.compareTo(module) >= BigInteger.EQUALS) {
0613: Elementary.inplaceSubtract(r, module);
0614: }
0615: r = module.subtract(r);
0616:
0617: return new Object[] { r, k };
0618: }
0619:
0620: /**
0621: * @return bi == abs(2^exp)
0622: */
0623: private static boolean isPowerOfTwo(BigInteger bi, int exp) {
0624: boolean result = false;
0625: result = (exp >> 5 == bi.numberLength - 1)
0626: && (bi.digits[bi.numberLength - 1] == 1 << (exp & 31));
0627: if (result) {
0628: for (int i = 0; result && i < bi.numberLength - 1; i++) {
0629: result = bi.digits[i] == 0;
0630: }
0631: }
0632: return result;
0633: }
0634:
0635: /**
0636: * Calculate how many iteration of Lorencz's algorithm would perform the
0637: * same operation
0638: *
0639: * @param bi
0640: * @param n
0641: * @return
0642: */
0643: private static int howManyIterations(BigInteger bi, int n) {
0644: int i = n - 1;
0645: if (bi.sign > 0) {
0646: while (!bi.testBit(i))
0647: i--;
0648: return n - 1 - i;
0649: } else {
0650: while (bi.testBit(i))
0651: i--;
0652: return n - 1 - Math.max(i, bi.getLowestSetBit());
0653: }
0654:
0655: }
0656:
0657: /**
0658: *
0659: * Based on "New Algorithm for Classical Modular Inverse" Róbert Lórencz.
0660: * LNCS 2523 (2002)
0661: *
0662: * @return a^(-1) mod m
0663: */
0664: static BigInteger modInverseLorencz(BigInteger a, BigInteger modulo) {
0665: // PRE: a is coprime with modulo, a < modulo
0666:
0667: int max = Math.max(a.numberLength, modulo.numberLength);
0668: int uDigits[] = new int[max + 1]; // enough place to make all the inplace operation
0669: int vDigits[] = new int[max + 1];
0670: System.arraycopy(modulo.digits, 0, uDigits, 0,
0671: modulo.numberLength);
0672: System.arraycopy(a.digits, 0, vDigits, 0, a.numberLength);
0673: BigInteger u = new BigInteger(modulo.sign, modulo.numberLength,
0674: uDigits);
0675: BigInteger v = new BigInteger(a.sign, a.numberLength, vDigits);
0676:
0677: BigInteger r = new BigInteger(0, 1, new int[max + 1]); // BigInteger.ZERO;
0678: BigInteger s = new BigInteger(1, 1, new int[max + 1]);
0679: s.digits[0] = 1;
0680: // r == 0 && s == 1, but with enough place
0681:
0682: int coefU = 0, coefV = 0;
0683: int n = modulo.bitLength();
0684: int k;
0685: while (!isPowerOfTwo(u, coefU) && !isPowerOfTwo(v, coefV)) {
0686:
0687: // modification of original algorithm: I calculate how many times the algorithm will enter in the same branch of if
0688: k = howManyIterations(u, n);
0689:
0690: if (k != 0) {
0691: BitLevel.inplaceShiftLeft(u, k);
0692: if (coefU >= coefV) {
0693: BitLevel.inplaceShiftLeft(r, k);
0694: } else {
0695: BitLevel.inplaceShiftRight(s, Math.min(coefV
0696: - coefU, k));
0697: if (k - (coefV - coefU) > 0) {
0698: BitLevel.inplaceShiftLeft(r, k - coefV + coefU);
0699: }
0700: }
0701: coefU += k;
0702: }
0703:
0704: k = howManyIterations(v, n);
0705: if (k != 0) {
0706: BitLevel.inplaceShiftLeft(v, k);
0707: if (coefV >= coefU) {
0708: BitLevel.inplaceShiftLeft(s, k);
0709: } else {
0710: BitLevel.inplaceShiftRight(r, Math.min(coefU
0711: - coefV, k));
0712: if (k - (coefU - coefV) > 0) {
0713: BitLevel.inplaceShiftLeft(s, k - coefU + coefV);
0714: }
0715: }
0716: coefV += k;
0717:
0718: }
0719:
0720: if (u.signum() == v.signum()) {
0721: if (coefU <= coefV) {
0722: Elementary.completeInPlaceSubtract(u, v);
0723: Elementary.completeInPlaceSubtract(r, s);
0724: } else {
0725: Elementary.completeInPlaceSubtract(v, u);
0726: Elementary.completeInPlaceSubtract(s, r);
0727: }
0728: } else {
0729: if (coefU <= coefV) {
0730: Elementary.completeInPlaceAdd(u, v);
0731: Elementary.completeInPlaceAdd(r, s);
0732: } else {
0733: Elementary.completeInPlaceAdd(v, u);
0734: Elementary.completeInPlaceAdd(s, r);
0735: }
0736: }
0737: if (v.signum() == 0 || u.signum() == 0) {
0738: // math.19: BigInteger not invertible
0739: throw new ArithmeticException(Messages
0740: .getString("math.19"));
0741: }
0742: }
0743:
0744: if (isPowerOfTwo(v, coefV)) {
0745: r = s;
0746: if (v.signum() != u.signum())
0747: u = u.negate();
0748: }
0749: if (u.testBit(n)) {
0750: if (r.signum() < 0) {
0751: r = r.negate();
0752: } else {
0753: r = modulo.subtract(r);
0754: }
0755: }
0756: if (r.signum() < 0) {
0757: r = r.add(modulo);
0758: }
0759:
0760: return r;
0761: }
0762:
0763: static BigInteger squareAndMultiply(BigInteger x2, BigInteger a2,
0764: BigInteger exponent, BigInteger modulus, long n2) {
0765: BigInteger res = x2;
0766: for (int i = exponent.bitLength() - 1; i >= 0; i--) {
0767: res = monPro(res, res, modulus, n2);
0768: if (BitLevel.testBit(exponent, i)) {
0769: res = monPro(res, a2, modulus, n2);
0770: }
0771: }
0772: return res;
0773: }
0774:
0775: /*Implements the Montgomery modular exponentiation based in <i>The sliding windows algorithm and the Mongomery
0776: *Reduction</i>.
0777: *@ar.org.fitc.ref "A. Menezes,P. van Oorschot, S. Vanstone - Handbook of Applied Cryptography";
0778: *@see #oddModPow(BigInteger, BigInteger,
0779: * BigInteger)
0780: */
0781: static BigInteger slidingWindow(BigInteger x2, BigInteger a2,
0782: BigInteger exponent, BigInteger modulus, long n2) {
0783: // fill odd low pows of a2
0784: BigInteger pows[] = new BigInteger[8];
0785: BigInteger res = x2;
0786: int lowexp;
0787: BigInteger x3;
0788: int acc3;
0789: pows[0] = a2;
0790:
0791: x3 = monSquare(a2, modulus, n2);
0792: for (int i = 1; i <= 7; i++) {
0793: pows[i] = monPro(pows[i - 1], x3, modulus, n2);
0794: }
0795:
0796: for (int i = exponent.bitLength() - 1; i >= 0; i--) {
0797: if (BitLevel.testBit(exponent, i)) {
0798: lowexp = 1;
0799: acc3 = i;
0800:
0801: for (int j = Math.max(i - 3, 0); j <= i - 1; j++) {
0802: if (BitLevel.testBit(exponent, j)) {
0803: if (j < acc3) {
0804: acc3 = j;
0805: lowexp = (lowexp << (i - j)) ^ 1;
0806: } else {
0807: lowexp = lowexp ^ (1 << (j - acc3));
0808: }
0809: }
0810: }
0811:
0812: for (int j = acc3; j <= i; j++) {
0813: res = monSquare(res, modulus, n2);
0814: }
0815: res = monPro(pows[(lowexp - 1) >> 1], res, modulus, n2);
0816: i = acc3;
0817: } else {
0818: res = monSquare(res, modulus, n2);
0819: }
0820: }
0821: return res;
0822: }
0823:
0824: /**
0825: * Performs modular exponentiation using the Montgomery Reduction. It
0826: * requires that all parameters be positive and the modulus be odd. >
0827: *
0828: * @see BigInteger#modPow(BigInteger, BigInteger)
0829: * @see #monPro(BigInteger, BigInteger, BigInteger, long)
0830: * @see #slidingWindow(BigInteger, BigInteger, BigInteger, BigInteger,
0831: * long)
0832: * @see #squareAndMultiply(BigInteger, BigInteger, BigInteger, BigInteger,
0833: * long)
0834: */
0835: static BigInteger oddModPow(BigInteger base, BigInteger exponent,
0836: BigInteger modulus) {
0837: // PRE: (base > 0), (exponent > 0), (modulus > 0) and (odd modulus)
0838: int k = (modulus.numberLength << 5); // r = 2^k
0839: // n-residue of base [base * r (mod modulus)]
0840: BigInteger a2 = base.shiftLeft(k).mod(modulus);
0841: // n-residue of base [1 * r (mod modulus)]
0842: BigInteger x2 = BigInteger.ZERO.setBit(k).mod(modulus);
0843: BigInteger res;
0844: // Compute (modulus[0]^(-1)) (mod 2^32) for odd modulus
0845:
0846: long m0 = modulus.digits[0] & 0xFFFFFFFFL;
0847: long n2 = 1L; // this is n'[0]
0848: long powerOfTwo = 2L;
0849: // compute n2
0850: do {
0851: if (((m0 * n2) & powerOfTwo) != 0) {
0852: n2 |= powerOfTwo;
0853: }
0854: powerOfTwo <<= 1;
0855: } while (powerOfTwo < 0x100000000L);
0856: n2 = -n2;
0857:
0858: if (modulus.numberLength == 1) {
0859: res = squareAndMultiply(x2, a2, exponent, modulus, n2);
0860: } else {
0861: res = slidingWindow(x2, a2, exponent, modulus, n2);
0862: }
0863:
0864: return monPro(res, BigInteger.ONE, modulus, n2);
0865: }
0866:
0867: /**
0868: * Performs modular exponentiation using the Montgomery Reduction. It
0869: * requires that all parameters be positive and the modulus be even. Based
0870: * <i>The square and multiply algorithm and the Montgomery Reduction C. K.
0871: * Koc - Montgomery Reduction with Even Modulus</i>. The square and
0872: * multiply algorithm and the Montgomery Reduction.
0873: *
0874: * @ar.org.fitc.ref "C. K. Koc - Montgomery Reduction with Even Modulus"
0875: * @see BigInteger#modPow(BigInteger, BigInteger)
0876: */
0877: static BigInteger evenModPow(BigInteger base, BigInteger exponent,
0878: BigInteger modulus) {
0879: // PRE: (base > 0), (exponent > 0), (modulus > 0) and (modulus even)
0880: // STEP 1: Obtain the factorization 'modulus'= q * 2^j.
0881: int j = modulus.getLowestSetBit();
0882: BigInteger q = modulus.shiftRight(j);
0883:
0884: // STEP 2: Compute x1 := base^exponent (mod q).
0885: BigInteger x1 = oddModPow(base, exponent, q);
0886:
0887: // STEP 3: Compute x2 := base^exponent (mod 2^j).
0888: BigInteger x2 = pow2ModPow(base, exponent, j);
0889:
0890: // STEP 4: Compute q^(-1) (mod 2^j) and y := (x2-x1) * q^(-1) (mod 2^j)
0891: BigInteger qInv = modPow2Inverse(q, j);
0892: BigInteger y = (x2.subtract(x1)).multiply(qInv);
0893: inplaceModPow2(y, j);
0894: if (y.sign < 0) {
0895: y = y.add(BigInteger.ZERO.setBit(j));
0896: }
0897: // STEP 5: Compute and return: x1 + q * y
0898: return x1.add(q.multiply(y));
0899: }
0900:
0901: /**
0902: * It requires that all parameters be positive.
0903: *
0904: * @return {@code base<sup>exponent</sup> mod (2<sup>j</sup>)}.
0905: * @see BigInteger#modPow(BigInteger, BigInteger)
0906: */
0907: static BigInteger pow2ModPow(BigInteger base, BigInteger exponent,
0908: int j) {
0909: // PRE: (base > 0), (exponent > 0) and (j > 0)
0910: BigInteger res = BigInteger.ONE;
0911: BigInteger e = exponent.copy();
0912: BigInteger baseMod2toN = base.copy();
0913: BigInteger res2;
0914: /*
0915: * If 'base' is odd then it's coprime with 2^j and phi(2^j) = 2^(j-1);
0916: * so we can reduce reduce the exponent (mod 2^(j-1)).
0917: */
0918: if (base.testBit(0)) {
0919: inplaceModPow2(e, j - 1);
0920: }
0921: inplaceModPow2(baseMod2toN, j);
0922:
0923: for (int i = e.bitLength() - 1; i >= 0; i--) {
0924: res2 = res.copy();
0925: inplaceModPow2(res2, j);
0926: res = res.multiply(res2);
0927: if (BitLevel.testBit(e, i)) {
0928: res = res.multiply(baseMod2toN);
0929: inplaceModPow2(res, j);
0930: }
0931: }
0932: inplaceModPow2(res, j);
0933: return res;
0934: }
0935:
0936: /** Implements the Montgomery Square of a BigInteger.
0937: * @see #monPro(BigInteger, BigInteger, BigInteger,
0938: * long)
0939: */
0940: static BigInteger monSquare(BigInteger aBig, BigInteger modulus,
0941: long n2) {
0942: if (modulus.numberLength == 1) {
0943: return monPro(aBig, aBig, modulus, n2);
0944: }
0945: //Squaring
0946: int[] a = aBig.digits;
0947: int[] n = modulus.digits;
0948: int s = modulus.numberLength;
0949:
0950: //Multiplying...
0951: int[] t = new int[(s << 1) + 1];
0952: long cs;
0953:
0954: int limit = Math.min(s, aBig.numberLength);
0955: for (int i = 0; i < limit; i++) {
0956: cs = 0;
0957: for (int j = i + 1; j < limit; j++) {
0958: cs += (0xFFFFFFFFL & t[i + j]) + (0xFFFFFFFFL & a[i])
0959: * (0xFFFFFFFFL & a[j]);
0960: t[i + j] = (int) cs;
0961: cs >>>= 32;
0962: }
0963:
0964: t[i + limit] = (int) cs;
0965: }
0966: BitLevel.shiftLeft(t, t, 0, 1);
0967:
0968: cs = 0;
0969: long carry = 0;
0970: for (int i = 0, index = 0; i < s; i++, index++) {
0971: cs += (0xFFFFFFFFL & a[i]) * (0xFFFFFFFFL & a[i])
0972: + (t[index] & 0xFFFFFFFFL);
0973: t[index] = (int) cs;
0974: cs >>>= 32;
0975: index++;
0976: cs += t[index] & 0xFFFFFFFFL;
0977: t[index] = (int) cs;
0978: cs >>>= 32;
0979: }
0980:
0981: //Reducing...
0982: /* t + m*n */
0983: int m = 0;
0984: int i, j;
0985: cs = carry = 0;
0986: for (i = 0; i < s; i++) {
0987: cs = 0;
0988: m = (int) ((t[i] & 0xFFFFFFFFL) * (n2 & 0xFFFFFFFFL));
0989: for (j = 0; j < s; j++) {
0990: cs = (t[i + j] & 0xFFFFFFFFL) + (m & 0xFFFFFFFFL)
0991: * (n[j] & 0xFFFFFFFFL) + (cs >>> 32);
0992: t[i + j] = (int) cs;
0993: }
0994: //Adding C to the result
0995: carry += (t[i + s] & 0xFFFFFFFFL)
0996: + ((cs >>> 32) & 0xFFFFFFFFL);
0997: t[i + s] = (int) carry;
0998: carry >>>= 32;
0999: }
1000:
1001: t[s << 1] = (int) carry;
1002:
1003: /* t / r */
1004: for (j = 0; j < s + 1; j++) {
1005: t[j] = t[j + s];
1006: }
1007: /*step 3*/
1008: return finalSubtraction(t, s, s, modulus);
1009: }
1010:
1011: /**
1012: * Implements the Montgomery Product of two integers represented by
1013: * {@code int} arrays. The arrays are supposed in <i>little
1014: * endian</i> notation.
1015: *
1016: * @param a The first factor of the product.
1017: * @param b The second factor of the product.
1018: * @param modulus The modulus of the operations. Z<sub>modulus</sub>.
1019: * @param n2 The digit modulus'[0].
1020: * @ar.org.fitc.ref "C. K. Koc - Analyzing and Comparing Montgomery
1021: * Multiplication Algorithms"
1022: * @see #modPowOdd(BigInteger, BigInteger, BigInteger)
1023: */
1024: static BigInteger monPro(BigInteger a, BigInteger b,
1025: BigInteger modulus, long n2) {
1026: int aFirst = a.numberLength - 1;
1027: int bFirst = b.numberLength - 1;
1028: int s = modulus.numberLength;
1029: int n[] = modulus.digits;
1030:
1031: int m;
1032: int i, j;
1033: int t[] = new int[(s << 1) + 1];
1034: long product;
1035: long C;
1036: long aI;
1037:
1038: for (i = 0; i < s; i++) {
1039: C = 0;
1040: aI = (i > aFirst) ? 0 : (a.digits[i] & 0xFFFFFFFFL);
1041: for (j = 0; j < s; j++) {
1042: product = (t[j] & 0xFFFFFFFFL)
1043: + C
1044: + ((j > bFirst) ? 0
1045: : (b.digits[j] & 0xFFFFFFFFL) * aI);
1046: C = product >>> 32;
1047: t[j] = (int) product;
1048: }
1049: product = (t[s] & 0xFFFFFFFFL) + C;
1050: C = product >>> 32;
1051: t[s] = (int) product;
1052: t[s + 1] = (int) C;
1053:
1054: m = (int) ((t[0] & 0xFFFFFFFFL) * n2);
1055: product = (t[0] & 0xFFFFFFFFL) + (m & 0xFFFFFFFFL)
1056: * (n[0] & 0xFFFFFFFFL);
1057: C = (int) (product >>> 32);
1058: for (j = 1; j < s; j++) {
1059: product = (t[j] & 0xFFFFFFFFL) + (C & 0xFFFFFFFFL)
1060: + (m & 0xFFFFFFFFL) * (n[j] & 0xFFFFFFFFL);
1061: C = product >>> 32;
1062: t[j - 1] = (int) product;
1063: }
1064: product = (t[s] & 0xFFFFFFFFL) + (C & 0xFFFFFFFFL);
1065: C = product >>> 32;
1066: t[s - 1] = (int) product;
1067: t[s] = t[s + 1] + (int) C;
1068: }
1069:
1070: return finalSubtraction(t, t.length - 1, s, modulus);
1071:
1072: }
1073:
1074: /*Performs the final reduction of the Montgomery algorithm.
1075: *@see monPro(BigInteger, BigInteger, BigInteger,
1076: *long )
1077: *@see monSquare(BigInteger, BigInteger ,
1078: *long)
1079: */
1080: static BigInteger finalSubtraction(int t[], int tLength, int s,
1081: BigInteger modulus) {
1082: // skipping leading zeros
1083: int i;
1084: int n[] = modulus.digits;
1085: boolean lower = false;
1086:
1087: for (i = tLength; (i > 0) && (t[i] == 0); i--)
1088: ;
1089:
1090: if (i == s - 1) {
1091: for (; (i >= 0) && (t[i] == n[i]); i--)
1092: ;
1093: lower = (i >= 0)
1094: && (t[i] & 0xFFFFFFFFL) < (n[i] & 0xFFFFFFFFL);
1095: } else {
1096: lower = (i < s - 1);
1097: }
1098: BigInteger res = new BigInteger(1, s + 1, t);
1099: // if (t >= n) compute (t - n)
1100: if (!lower) {
1101: Elementary.inplaceSubtract(res, modulus);
1102: }
1103: res.cutOffLeadingZeroes();
1104: return res;
1105: }
1106:
1107: /**
1108: * @param x an odd positive number.
1109: * @param n the exponent by which 2 is raised.
1110: * @return {@code x<sup>-1</sup> (mod 2<sup>n</sup>)}.
1111: */
1112: static BigInteger modPow2Inverse(BigInteger x, int n) {
1113: // PRE: (x > 0), (x is odd), and (n > 0)
1114: BigInteger y = new BigInteger(1, new int[1 << n]);
1115: y.numberLength = 1;
1116: y.digits[0] = 1;
1117: y.sign = 1;
1118:
1119: for (int i = 1; i < n; i++) {
1120: if (BitLevel.testBit(x.multiply(y), i)) {
1121: // Adding 2^i to y (setting the i-th bit)
1122: y.digits[i >> 5] |= (1 << (i & 31));
1123: }
1124: }
1125: return y;
1126: }
1127:
1128: /**
1129: * Performs {@code x = x mod (2<sup>n</sup>)}.
1130: *
1131: * @param x a positive number, it will store the result.
1132: * @param n a positive exponent of {@code 2}.
1133: */
1134: static void inplaceModPow2(BigInteger x, int n) {
1135: // PRE: (x > 0) and (n >= 0)
1136: int fd = n >> 5;
1137: int leadingZeros;
1138:
1139: if ((x.numberLength < fd) || (x.bitLength() <= n)) {
1140: return;
1141: }
1142: leadingZeros = 32 - (n & 31);
1143: x.numberLength = fd + 1;
1144: x.digits[fd] &= (leadingZeros < 32) ? (-1 >>> leadingZeros) : 0;
1145: x.cutOffLeadingZeroes();
1146: }
1147:
1148: }
|