0001 /*
0002 * Copyright 1999-2007 Sun Microsystems, Inc. All Rights Reserved.
0003 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
0004 *
0005 * This code is free software; you can redistribute it and/or modify it
0006 * under the terms of the GNU General Public License version 2 only, as
0007 * published by the Free Software Foundation. Sun designates this
0008 * particular file as subject to the "Classpath" exception as provided
0009 * by Sun in the LICENSE file that accompanied this code.
0010 *
0011 * This code is distributed in the hope that it will be useful, but WITHOUT
0012 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
0013 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
0014 * version 2 for more details (a copy is included in the LICENSE file that
0015 * accompanied this code).
0016 *
0017 * You should have received a copy of the GNU General Public License version
0018 * 2 along with this work; if not, write to the Free Software Foundation,
0019 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
0020 *
0021 * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
0022 * CA 95054 USA or visit www.sun.com if you need additional information or
0023 * have any questions.
0024 */
0025
0026 package java.math;
0027
0028 /**
0029 * A class used to represent multiprecision integers that makes efficient
0030 * use of allocated space by allowing a number to occupy only part of
0031 * an array so that the arrays do not have to be reallocated as often.
0032 * When performing an operation with many iterations the array used to
0033 * hold a number is only reallocated when necessary and does not have to
0034 * be the same size as the number it represents. A mutable number allows
0035 * calculations to occur on the same number without having to create
0036 * a new number for every step of the calculation as occurs with
0037 * BigIntegers.
0038 *
0039 * @see BigInteger
0040 * @author Michael McCloskey
0041 * @since 1.3
0042 */
0043
0044 class MutableBigInteger {
0045 /**
0046 * Holds the magnitude of this MutableBigInteger in big endian order.
0047 * The magnitude may start at an offset into the value array, and it may
0048 * end before the length of the value array.
0049 */
0050 int[] value;
0051
0052 /**
0053 * The number of ints of the value array that are currently used
0054 * to hold the magnitude of this MutableBigInteger. The magnitude starts
0055 * at an offset and offset + intLen may be less than value.length.
0056 */
0057 int intLen;
0058
0059 /**
0060 * The offset into the value array where the magnitude of this
0061 * MutableBigInteger begins.
0062 */
0063 int offset = 0;
0064
0065 /**
0066 * This mask is used to obtain the value of an int as if it were unsigned.
0067 */
0068 private final static long LONG_MASK = 0xffffffffL;
0069
0070 // Constructors
0071
0072 /**
0073 * The default constructor. An empty MutableBigInteger is created with
0074 * a one word capacity.
0075 */
0076 MutableBigInteger() {
0077 value = new int[1];
0078 intLen = 0;
0079 }
0080
0081 /**
0082 * Construct a new MutableBigInteger with a magnitude specified by
0083 * the int val.
0084 */
0085 MutableBigInteger(int val) {
0086 value = new int[1];
0087 intLen = 1;
0088 value[0] = val;
0089 }
0090
0091 /**
0092 * Construct a new MutableBigInteger with the specified value array
0093 * up to the specified length.
0094 */
0095 MutableBigInteger(int[] val, int len) {
0096 value = val;
0097 intLen = len;
0098 }
0099
0100 /**
0101 * Construct a new MutableBigInteger with the specified value array
0102 * up to the length of the array supplied.
0103 */
0104 MutableBigInteger(int[] val) {
0105 value = val;
0106 intLen = val.length;
0107 }
0108
0109 /**
0110 * Construct a new MutableBigInteger with a magnitude equal to the
0111 * specified BigInteger.
0112 */
0113 MutableBigInteger(BigInteger b) {
0114 value = (int[]) b.mag.clone();
0115 intLen = value.length;
0116 }
0117
0118 /**
0119 * Construct a new MutableBigInteger with a magnitude equal to the
0120 * specified MutableBigInteger.
0121 */
0122 MutableBigInteger(MutableBigInteger val) {
0123 intLen = val.intLen;
0124 value = new int[intLen];
0125
0126 for (int i = 0; i < intLen; i++)
0127 value[i] = val.value[val.offset + i];
0128 }
0129
0130 /**
0131 * Clear out a MutableBigInteger for reuse.
0132 */
0133 void clear() {
0134 offset = intLen = 0;
0135 for (int index = 0, n = value.length; index < n; index++)
0136 value[index] = 0;
0137 }
0138
0139 /**
0140 * Set a MutableBigInteger to zero, removing its offset.
0141 */
0142 void reset() {
0143 offset = intLen = 0;
0144 }
0145
0146 /**
0147 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
0148 * as this MutableBigInteger is numerically less than, equal to, or
0149 * greater than {@code b}.
0150 */
0151 final int compare(MutableBigInteger b) {
0152 if (intLen < b.intLen)
0153 return -1;
0154 if (intLen > b.intLen)
0155 return 1;
0156
0157 for (int i = 0; i < intLen; i++) {
0158 int b1 = value[offset + i] + 0x80000000;
0159 int b2 = b.value[b.offset + i] + 0x80000000;
0160 if (b1 < b2)
0161 return -1;
0162 if (b1 > b2)
0163 return 1;
0164 }
0165 return 0;
0166 }
0167
0168 /**
0169 * Return the index of the lowest set bit in this MutableBigInteger. If the
0170 * magnitude of this MutableBigInteger is zero, -1 is returned.
0171 */
0172 private final int getLowestSetBit() {
0173 if (intLen == 0)
0174 return -1;
0175 int j, b;
0176 for (j = intLen - 1; (j > 0) && (value[j + offset] == 0); j--)
0177 ;
0178 b = value[j + offset];
0179 if (b == 0)
0180 return -1;
0181 return ((intLen - 1 - j) << 5) + BigInteger.trailingZeroCnt(b);
0182 }
0183
0184 /**
0185 * Return the int in use in this MutableBigInteger at the specified
0186 * index. This method is not used because it is not inlined on all
0187 * platforms.
0188 */
0189 private final int getInt(int index) {
0190 return value[offset + index];
0191 }
0192
0193 /**
0194 * Return a long which is equal to the unsigned value of the int in
0195 * use in this MutableBigInteger at the specified index. This method is
0196 * not used because it is not inlined on all platforms.
0197 */
0198 private final long getLong(int index) {
0199 return value[offset + index] & LONG_MASK;
0200 }
0201
0202 /**
0203 * Ensure that the MutableBigInteger is in normal form, specifically
0204 * making sure that there are no leading zeros, and that if the
0205 * magnitude is zero, then intLen is zero.
0206 */
0207 final void normalize() {
0208 if (intLen == 0) {
0209 offset = 0;
0210 return;
0211 }
0212
0213 int index = offset;
0214 if (value[index] != 0)
0215 return;
0216
0217 int indexBound = index + intLen;
0218 do {
0219 index++;
0220 } while (index < indexBound && value[index] == 0);
0221
0222 int numZeros = index - offset;
0223 intLen -= numZeros;
0224 offset = (intLen == 0 ? 0 : offset + numZeros);
0225 }
0226
0227 /**
0228 * If this MutableBigInteger cannot hold len words, increase the size
0229 * of the value array to len words.
0230 */
0231 private final void ensureCapacity(int len) {
0232 if (value.length < len) {
0233 value = new int[len];
0234 offset = 0;
0235 intLen = len;
0236 }
0237 }
0238
0239 /**
0240 * Convert this MutableBigInteger into an int array with no leading
0241 * zeros, of a length that is equal to this MutableBigInteger's intLen.
0242 */
0243 int[] toIntArray() {
0244 int[] result = new int[intLen];
0245 for (int i = 0; i < intLen; i++)
0246 result[i] = value[offset + i];
0247 return result;
0248 }
0249
0250 /**
0251 * Sets the int at index+offset in this MutableBigInteger to val.
0252 * This does not get inlined on all platforms so it is not used
0253 * as often as originally intended.
0254 */
0255 void setInt(int index, int val) {
0256 value[offset + index] = val;
0257 }
0258
0259 /**
0260 * Sets this MutableBigInteger's value array to the specified array.
0261 * The intLen is set to the specified length.
0262 */
0263 void setValue(int[] val, int length) {
0264 value = val;
0265 intLen = length;
0266 offset = 0;
0267 }
0268
0269 /**
0270 * Sets this MutableBigInteger's value array to a copy of the specified
0271 * array. The intLen is set to the length of the new array.
0272 */
0273 void copyValue(MutableBigInteger val) {
0274 int len = val.intLen;
0275 if (value.length < len)
0276 value = new int[len];
0277
0278 for (int i = 0; i < len; i++)
0279 value[i] = val.value[val.offset + i];
0280 intLen = len;
0281 offset = 0;
0282 }
0283
0284 /**
0285 * Sets this MutableBigInteger's value array to a copy of the specified
0286 * array. The intLen is set to the length of the specified array.
0287 */
0288 void copyValue(int[] val) {
0289 int len = val.length;
0290 if (value.length < len)
0291 value = new int[len];
0292 for (int i = 0; i < len; i++)
0293 value[i] = val[i];
0294 intLen = len;
0295 offset = 0;
0296 }
0297
0298 /**
0299 * Returns true iff this MutableBigInteger has a value of one.
0300 */
0301 boolean isOne() {
0302 return (intLen == 1) && (value[offset] == 1);
0303 }
0304
0305 /**
0306 * Returns true iff this MutableBigInteger has a value of zero.
0307 */
0308 boolean isZero() {
0309 return (intLen == 0);
0310 }
0311
0312 /**
0313 * Returns true iff this MutableBigInteger is even.
0314 */
0315 boolean isEven() {
0316 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
0317 }
0318
0319 /**
0320 * Returns true iff this MutableBigInteger is odd.
0321 */
0322 boolean isOdd() {
0323 return ((value[offset + intLen - 1] & 1) == 1);
0324 }
0325
0326 /**
0327 * Returns true iff this MutableBigInteger is in normal form. A
0328 * MutableBigInteger is in normal form if it has no leading zeros
0329 * after the offset, and intLen + offset <= value.length.
0330 */
0331 boolean isNormal() {
0332 if (intLen + offset > value.length)
0333 return false;
0334 if (intLen == 0)
0335 return true;
0336 return (value[offset] != 0);
0337 }
0338
0339 /**
0340 * Returns a String representation of this MutableBigInteger in radix 10.
0341 */
0342 public String toString() {
0343 BigInteger b = new BigInteger(this , 1);
0344 return b.toString();
0345 }
0346
0347 /**
0348 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
0349 * in normal form.
0350 */
0351 void rightShift(int n) {
0352 if (intLen == 0)
0353 return;
0354 int nInts = n >>> 5;
0355 int nBits = n & 0x1F;
0356 this .intLen -= nInts;
0357 if (nBits == 0)
0358 return;
0359 int bitsInHighWord = BigInteger.bitLen(value[offset]);
0360 if (nBits >= bitsInHighWord) {
0361 this .primitiveLeftShift(32 - nBits);
0362 this .intLen--;
0363 } else {
0364 primitiveRightShift(nBits);
0365 }
0366 }
0367
0368 /**
0369 * Left shift this MutableBigInteger n bits.
0370 */
0371 void leftShift(int n) {
0372 /*
0373 * If there is enough storage space in this MutableBigInteger already
0374 * the available space will be used. Space to the right of the used
0375 * ints in the value array is faster to utilize, so the extra space
0376 * will be taken from the right if possible.
0377 */
0378 if (intLen == 0)
0379 return;
0380 int nInts = n >>> 5;
0381 int nBits = n & 0x1F;
0382 int bitsInHighWord = BigInteger.bitLen(value[offset]);
0383
0384 // If shift can be done without moving words, do so
0385 if (n <= (32 - bitsInHighWord)) {
0386 primitiveLeftShift(nBits);
0387 return;
0388 }
0389
0390 int newLen = intLen + nInts + 1;
0391 if (nBits <= (32 - bitsInHighWord))
0392 newLen--;
0393 if (value.length < newLen) {
0394 // The array must grow
0395 int[] result = new int[newLen];
0396 for (int i = 0; i < intLen; i++)
0397 result[i] = value[offset + i];
0398 setValue(result, newLen);
0399 } else if (value.length - offset >= newLen) {
0400 // Use space on right
0401 for (int i = 0; i < newLen - intLen; i++)
0402 value[offset + intLen + i] = 0;
0403 } else {
0404 // Must use space on left
0405 for (int i = 0; i < intLen; i++)
0406 value[i] = value[offset + i];
0407 for (int i = intLen; i < newLen; i++)
0408 value[i] = 0;
0409 offset = 0;
0410 }
0411 intLen = newLen;
0412 if (nBits == 0)
0413 return;
0414 if (nBits <= (32 - bitsInHighWord))
0415 primitiveLeftShift(nBits);
0416 else
0417 primitiveRightShift(32 - nBits);
0418 }
0419
0420 /**
0421 * A primitive used for division. This method adds in one multiple of the
0422 * divisor a back to the dividend result at a specified offset. It is used
0423 * when qhat was estimated too large, and must be adjusted.
0424 */
0425 private int divadd(int[] a, int[] result, int offset) {
0426 long carry = 0;
0427
0428 for (int j = a.length - 1; j >= 0; j--) {
0429 long sum = (a[j] & LONG_MASK)
0430 + (result[j + offset] & LONG_MASK) + carry;
0431 result[j + offset] = (int) sum;
0432 carry = sum >>> 32;
0433 }
0434 return (int) carry;
0435 }
0436
0437 /**
0438 * This method is used for division. It multiplies an n word input a by one
0439 * word input x, and subtracts the n word product from q. This is needed
0440 * when subtracting qhat*divisor from dividend.
0441 */
0442 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
0443 long xLong = x & LONG_MASK;
0444 long carry = 0;
0445 offset += len;
0446
0447 for (int j = len - 1; j >= 0; j--) {
0448 long product = (a[j] & LONG_MASK) * xLong + carry;
0449 long difference = q[offset] - product;
0450 q[offset--] = (int) difference;
0451 carry = (product >>> 32)
0452 + (((difference & LONG_MASK) > (((~(int) product) & LONG_MASK))) ? 1
0453 : 0);
0454 }
0455 return (int) carry;
0456 }
0457
0458 /**
0459 * Right shift this MutableBigInteger n bits, where n is
0460 * less than 32.
0461 * Assumes that intLen > 0, n > 0 for speed
0462 */
0463 private final void primitiveRightShift(int n) {
0464 int[] val = value;
0465 int n2 = 32 - n;
0466 for (int i = offset + intLen - 1, c = val[i]; i > offset; i--) {
0467 int b = c;
0468 c = val[i - 1];
0469 val[i] = (c << n2) | (b >>> n);
0470 }
0471 val[offset] >>>= n;
0472 }
0473
0474 /**
0475 * Left shift this MutableBigInteger n bits, where n is
0476 * less than 32.
0477 * Assumes that intLen > 0, n > 0 for speed
0478 */
0479 private final void primitiveLeftShift(int n) {
0480 int[] val = value;
0481 int n2 = 32 - n;
0482 for (int i = offset, c = val[i], m = i + intLen - 1; i < m; i++) {
0483 int b = c;
0484 c = val[i + 1];
0485 val[i] = (b << n) | (c >>> n2);
0486 }
0487 val[offset + intLen - 1] <<= n;
0488 }
0489
0490 /**
0491 * Adds the contents of two MutableBigInteger objects.The result
0492 * is placed within this MutableBigInteger.
0493 * The contents of the addend are not changed.
0494 */
0495 void add(MutableBigInteger addend) {
0496 int x = intLen;
0497 int y = addend.intLen;
0498 int resultLen = (intLen > addend.intLen ? intLen
0499 : addend.intLen);
0500 int[] result = (value.length < resultLen ? new int[resultLen]
0501 : value);
0502
0503 int rstart = result.length - 1;
0504 long sum = 0;
0505
0506 // Add common parts of both numbers
0507 while (x > 0 && y > 0) {
0508 x--;
0509 y--;
0510 sum = (value[x + offset] & LONG_MASK)
0511 + (addend.value[y + addend.offset] & LONG_MASK)
0512 + (sum >>> 32);
0513 result[rstart--] = (int) sum;
0514 }
0515
0516 // Add remainder of the longer number
0517 while (x > 0) {
0518 x--;
0519 sum = (value[x + offset] & LONG_MASK) + (sum >>> 32);
0520 result[rstart--] = (int) sum;
0521 }
0522 while (y > 0) {
0523 y--;
0524 sum = (addend.value[y + addend.offset] & LONG_MASK)
0525 + (sum >>> 32);
0526 result[rstart--] = (int) sum;
0527 }
0528
0529 if ((sum >>> 32) > 0) { // Result must grow in length
0530 resultLen++;
0531 if (result.length < resultLen) {
0532 int temp[] = new int[resultLen];
0533 for (int i = resultLen - 1; i > 0; i--)
0534 temp[i] = result[i - 1];
0535 temp[0] = 1;
0536 result = temp;
0537 } else {
0538 result[rstart--] = 1;
0539 }
0540 }
0541
0542 value = result;
0543 intLen = resultLen;
0544 offset = result.length - resultLen;
0545 }
0546
0547 /**
0548 * Subtracts the smaller of this and b from the larger and places the
0549 * result into this MutableBigInteger.
0550 */
0551 int subtract(MutableBigInteger b) {
0552 MutableBigInteger a = this ;
0553
0554 int[] result = value;
0555 int sign = a.compare(b);
0556
0557 if (sign == 0) {
0558 reset();
0559 return 0;
0560 }
0561 if (sign < 0) {
0562 MutableBigInteger tmp = a;
0563 a = b;
0564 b = tmp;
0565 }
0566
0567 int resultLen = a.intLen;
0568 if (result.length < resultLen)
0569 result = new int[resultLen];
0570
0571 long diff = 0;
0572 int x = a.intLen;
0573 int y = b.intLen;
0574 int rstart = result.length - 1;
0575
0576 // Subtract common parts of both numbers
0577 while (y > 0) {
0578 x--;
0579 y--;
0580
0581 diff = (a.value[x + a.offset] & LONG_MASK)
0582 - (b.value[y + b.offset] & LONG_MASK)
0583 - ((int) -(diff >> 32));
0584 result[rstart--] = (int) diff;
0585 }
0586 // Subtract remainder of longer number
0587 while (x > 0) {
0588 x--;
0589 diff = (a.value[x + a.offset] & LONG_MASK)
0590 - ((int) -(diff >> 32));
0591 result[rstart--] = (int) diff;
0592 }
0593
0594 value = result;
0595 intLen = resultLen;
0596 offset = value.length - resultLen;
0597 normalize();
0598 return sign;
0599 }
0600
0601 /**
0602 * Subtracts the smaller of a and b from the larger and places the result
0603 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
0604 * operation was performed.
0605 */
0606 private int difference(MutableBigInteger b) {
0607 MutableBigInteger a = this ;
0608 int sign = a.compare(b);
0609 if (sign == 0)
0610 return 0;
0611 if (sign < 0) {
0612 MutableBigInteger tmp = a;
0613 a = b;
0614 b = tmp;
0615 }
0616
0617 long diff = 0;
0618 int x = a.intLen;
0619 int y = b.intLen;
0620
0621 // Subtract common parts of both numbers
0622 while (y > 0) {
0623 x--;
0624 y--;
0625 diff = (a.value[a.offset + x] & LONG_MASK)
0626 - (b.value[b.offset + y] & LONG_MASK)
0627 - ((int) -(diff >> 32));
0628 a.value[a.offset + x] = (int) diff;
0629 }
0630 // Subtract remainder of longer number
0631 while (x > 0) {
0632 x--;
0633 diff = (a.value[a.offset + x] & LONG_MASK)
0634 - ((int) -(diff >> 32));
0635 a.value[a.offset + x] = (int) diff;
0636 }
0637
0638 a.normalize();
0639 return sign;
0640 }
0641
0642 /**
0643 * Multiply the contents of two MutableBigInteger objects. The result is
0644 * placed into MutableBigInteger z. The contents of y are not changed.
0645 */
0646 void multiply(MutableBigInteger y, MutableBigInteger z) {
0647 int xLen = intLen;
0648 int yLen = y.intLen;
0649 int newLen = xLen + yLen;
0650
0651 // Put z into an appropriate state to receive product
0652 if (z.value.length < newLen)
0653 z.value = new int[newLen];
0654 z.offset = 0;
0655 z.intLen = newLen;
0656
0657 // The first iteration is hoisted out of the loop to avoid extra add
0658 long carry = 0;
0659 for (int j = yLen - 1, k = yLen + xLen - 1; j >= 0; j--, k--) {
0660 long product = (y.value[j + y.offset] & LONG_MASK)
0661 * (value[xLen - 1 + offset] & LONG_MASK) + carry;
0662 z.value[k] = (int) product;
0663 carry = product >>> 32;
0664 }
0665 z.value[xLen - 1] = (int) carry;
0666
0667 // Perform the multiplication word by word
0668 for (int i = xLen - 2; i >= 0; i--) {
0669 carry = 0;
0670 for (int j = yLen - 1, k = yLen + i; j >= 0; j--, k--) {
0671 long product = (y.value[j + y.offset] & LONG_MASK)
0672 * (value[i + offset] & LONG_MASK)
0673 + (z.value[k] & LONG_MASK) + carry;
0674 z.value[k] = (int) product;
0675 carry = product >>> 32;
0676 }
0677 z.value[i] = (int) carry;
0678 }
0679
0680 // Remove leading zeros from product
0681 z.normalize();
0682 }
0683
0684 /**
0685 * Multiply the contents of this MutableBigInteger by the word y. The
0686 * result is placed into z.
0687 */
0688 void mul(int y, MutableBigInteger z) {
0689 if (y == 1) {
0690 z.copyValue(this );
0691 return;
0692 }
0693
0694 if (y == 0) {
0695 z.clear();
0696 return;
0697 }
0698
0699 // Perform the multiplication word by word
0700 long ylong = y & LONG_MASK;
0701 int[] zval = (z.value.length < intLen + 1 ? new int[intLen + 1]
0702 : z.value);
0703 long carry = 0;
0704 for (int i = intLen - 1; i >= 0; i--) {
0705 long product = ylong * (value[i + offset] & LONG_MASK)
0706 + carry;
0707 zval[i + 1] = (int) product;
0708 carry = product >>> 32;
0709 }
0710
0711 if (carry == 0) {
0712 z.offset = 1;
0713 z.intLen = intLen;
0714 } else {
0715 z.offset = 0;
0716 z.intLen = intLen + 1;
0717 zval[0] = (int) carry;
0718 }
0719 z.value = zval;
0720 }
0721
0722 /**
0723 * This method is used for division of an n word dividend by a one word
0724 * divisor. The quotient is placed into quotient. The one word divisor is
0725 * specified by divisor. The value of this MutableBigInteger is the
0726 * dividend at invocation but is replaced by the remainder.
0727 *
0728 * NOTE: The value of this MutableBigInteger is modified by this method.
0729 */
0730 void divideOneWord(int divisor, MutableBigInteger quotient) {
0731 long divLong = divisor & LONG_MASK;
0732
0733 // Special case of one word dividend
0734 if (intLen == 1) {
0735 long remValue = value[offset] & LONG_MASK;
0736 quotient.value[0] = (int) (remValue / divLong);
0737 quotient.intLen = (quotient.value[0] == 0) ? 0 : 1;
0738 quotient.offset = 0;
0739
0740 value[0] = (int) (remValue - (quotient.value[0] * divLong));
0741 offset = 0;
0742 intLen = (value[0] == 0) ? 0 : 1;
0743
0744 return;
0745 }
0746
0747 if (quotient.value.length < intLen)
0748 quotient.value = new int[intLen];
0749 quotient.offset = 0;
0750 quotient.intLen = intLen;
0751
0752 // Normalize the divisor
0753 int shift = 32 - BigInteger.bitLen(divisor);
0754
0755 int rem = value[offset];
0756 long remLong = rem & LONG_MASK;
0757 if (remLong < divLong) {
0758 quotient.value[0] = 0;
0759 } else {
0760 quotient.value[0] = (int) (remLong / divLong);
0761 rem = (int) (remLong - (quotient.value[0] * divLong));
0762 remLong = rem & LONG_MASK;
0763 }
0764
0765 int xlen = intLen;
0766 int[] qWord = new int[2];
0767 while (--xlen > 0) {
0768 long dividendEstimate = (remLong << 32)
0769 | (value[offset + intLen - xlen] & LONG_MASK);
0770 if (dividendEstimate >= 0) {
0771 qWord[0] = (int) (dividendEstimate / divLong);
0772 qWord[1] = (int) (dividendEstimate - (qWord[0] * divLong));
0773 } else {
0774 divWord(qWord, dividendEstimate, divisor);
0775 }
0776 quotient.value[intLen - xlen] = (int) qWord[0];
0777 rem = (int) qWord[1];
0778 remLong = rem & LONG_MASK;
0779 }
0780
0781 // Unnormalize
0782 if (shift > 0)
0783 value[0] = rem %= divisor;
0784 else
0785 value[0] = rem;
0786 intLen = (value[0] == 0) ? 0 : 1;
0787
0788 quotient.normalize();
0789 }
0790
0791 /**
0792 * Calculates the quotient and remainder of this div b and places them
0793 * in the MutableBigInteger objects provided.
0794 *
0795 * Uses Algorithm D in Knuth section 4.3.1.
0796 * Many optimizations to that algorithm have been adapted from the Colin
0797 * Plumb C library.
0798 * It special cases one word divisors for speed.
0799 * The contents of a and b are not changed.
0800 *
0801 */
0802 void divide(MutableBigInteger b, MutableBigInteger quotient,
0803 MutableBigInteger rem) {
0804 if (b.intLen == 0)
0805 throw new ArithmeticException("BigInteger divide by zero");
0806
0807 // Dividend is zero
0808 if (intLen == 0) {
0809 quotient.intLen = quotient.offset = rem.intLen = rem.offset = 0;
0810 return;
0811 }
0812
0813 int cmp = compare(b);
0814
0815 // Dividend less than divisor
0816 if (cmp < 0) {
0817 quotient.intLen = quotient.offset = 0;
0818 rem.copyValue(this );
0819 return;
0820 }
0821 // Dividend equal to divisor
0822 if (cmp == 0) {
0823 quotient.value[0] = quotient.intLen = 1;
0824 quotient.offset = rem.intLen = rem.offset = 0;
0825 return;
0826 }
0827
0828 quotient.clear();
0829
0830 // Special case one word divisor
0831 if (b.intLen == 1) {
0832 rem.copyValue(this );
0833 rem.divideOneWord(b.value[b.offset], quotient);
0834 return;
0835 }
0836
0837 // Copy divisor value to protect divisor
0838 int[] d = new int[b.intLen];
0839 for (int i = 0; i < b.intLen; i++)
0840 d[i] = b.value[b.offset + i];
0841 int dlen = b.intLen;
0842
0843 // Remainder starts as dividend with space for a leading zero
0844 if (rem.value.length < intLen + 1)
0845 rem.value = new int[intLen + 1];
0846
0847 for (int i = 0; i < intLen; i++)
0848 rem.value[i + 1] = value[i + offset];
0849 rem.intLen = intLen;
0850 rem.offset = 1;
0851
0852 int nlen = rem.intLen;
0853
0854 // Set the quotient size
0855 int limit = nlen - dlen + 1;
0856 if (quotient.value.length < limit) {
0857 quotient.value = new int[limit];
0858 quotient.offset = 0;
0859 }
0860 quotient.intLen = limit;
0861 int[] q = quotient.value;
0862
0863 // D1 normalize the divisor
0864 int shift = 32 - BigInteger.bitLen(d[0]);
0865 if (shift > 0) {
0866 // First shift will not grow array
0867 BigInteger.primitiveLeftShift(d, dlen, shift);
0868 // But this one might
0869 rem.leftShift(shift);
0870 }
0871
0872 // Must insert leading 0 in rem if its length did not change
0873 if (rem.intLen == nlen) {
0874 rem.offset = 0;
0875 rem.value[0] = 0;
0876 rem.intLen++;
0877 }
0878
0879 int dh = d[0];
0880 long dhLong = dh & LONG_MASK;
0881 int dl = d[1];
0882 int[] qWord = new int[2];
0883
0884 // D2 Initialize j
0885 for (int j = 0; j < limit; j++) {
0886 // D3 Calculate qhat
0887 // estimate qhat
0888 int qhat = 0;
0889 int qrem = 0;
0890 boolean skipCorrection = false;
0891 int nh = rem.value[j + rem.offset];
0892 int nh2 = nh + 0x80000000;
0893 int nm = rem.value[j + 1 + rem.offset];
0894
0895 if (nh == dh) {
0896 qhat = ~0;
0897 qrem = nh + nm;
0898 skipCorrection = qrem + 0x80000000 < nh2;
0899 } else {
0900 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
0901 if (nChunk >= 0) {
0902 qhat = (int) (nChunk / dhLong);
0903 qrem = (int) (nChunk - (qhat * dhLong));
0904 } else {
0905 divWord(qWord, nChunk, dh);
0906 qhat = qWord[0];
0907 qrem = qWord[1];
0908 }
0909 }
0910
0911 if (qhat == 0)
0912 continue;
0913
0914 if (!skipCorrection) { // Correct qhat
0915 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
0916 long rs = ((qrem & LONG_MASK) << 32) | nl;
0917 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
0918
0919 if (unsignedLongCompare(estProduct, rs)) {
0920 qhat--;
0921 qrem = (int) ((qrem & LONG_MASK) + dhLong);
0922 if ((qrem & LONG_MASK) >= dhLong) {
0923 estProduct = (dl & LONG_MASK)
0924 * (qhat & LONG_MASK);
0925 rs = ((qrem & LONG_MASK) << 32) | nl;
0926 if (unsignedLongCompare(estProduct, rs))
0927 qhat--;
0928 }
0929 }
0930 }
0931
0932 // D4 Multiply and subtract
0933 rem.value[j + rem.offset] = 0;
0934 int borrow = mulsub(rem.value, d, qhat, dlen, j
0935 + rem.offset);
0936
0937 // D5 Test remainder
0938 if (borrow + 0x80000000 > nh2) {
0939 // D6 Add back
0940 divadd(d, rem.value, j + 1 + rem.offset);
0941 qhat--;
0942 }
0943
0944 // Store the quotient digit
0945 q[j] = qhat;
0946 } // D7 loop on j
0947
0948 // D8 Unnormalize
0949 if (shift > 0)
0950 rem.rightShift(shift);
0951
0952 rem.normalize();
0953 quotient.normalize();
0954 }
0955
0956 /**
0957 * Compare two longs as if they were unsigned.
0958 * Returns true iff one is bigger than two.
0959 */
0960 private boolean unsignedLongCompare(long one, long two) {
0961 return (one + Long.MIN_VALUE) > (two + Long.MIN_VALUE);
0962 }
0963
0964 /**
0965 * This method divides a long quantity by an int to estimate
0966 * qhat for two multi precision numbers. It is used when
0967 * the signed value of n is less than zero.
0968 */
0969 private void divWord(int[] result, long n, int d) {
0970 long dLong = d & LONG_MASK;
0971
0972 if (dLong == 1) {
0973 result[0] = (int) n;
0974 result[1] = 0;
0975 return;
0976 }
0977
0978 // Approximate the quotient and remainder
0979 long q = (n >>> 1) / (dLong >>> 1);
0980 long r = n - q * dLong;
0981
0982 // Correct the approximation
0983 while (r < 0) {
0984 r += dLong;
0985 q--;
0986 }
0987 while (r >= dLong) {
0988 r -= dLong;
0989 q++;
0990 }
0991
0992 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
0993 result[0] = (int) q;
0994 result[1] = (int) r;
0995 }
0996
0997 /**
0998 * Calculate GCD of this and b. This and b are changed by the computation.
0999 */
1000 MutableBigInteger hybridGCD(MutableBigInteger b) {
1001 // Use Euclid's algorithm until the numbers are approximately the
1002 // same length, then use the binary GCD algorithm to find the GCD.
1003 MutableBigInteger a = this ;
1004 MutableBigInteger q = new MutableBigInteger(), r = new MutableBigInteger();
1005
1006 while (b.intLen != 0) {
1007 if (Math.abs(a.intLen - b.intLen) < 2)
1008 return a.binaryGCD(b);
1009
1010 a.divide(b, q, r);
1011 MutableBigInteger swapper = a;
1012 a = b;
1013 b = r;
1014 r = swapper;
1015 }
1016 return a;
1017 }
1018
1019 /**
1020 * Calculate GCD of this and v.
1021 * Assumes that this and v are not zero.
1022 */
1023 private MutableBigInteger binaryGCD(MutableBigInteger v) {
1024 // Algorithm B from Knuth section 4.5.2
1025 MutableBigInteger u = this ;
1026 MutableBigInteger r = new MutableBigInteger();
1027
1028 // step B1
1029 int s1 = u.getLowestSetBit();
1030 int s2 = v.getLowestSetBit();
1031 int k = (s1 < s2) ? s1 : s2;
1032 if (k != 0) {
1033 u.rightShift(k);
1034 v.rightShift(k);
1035 }
1036
1037 // step B2
1038 boolean uOdd = (k == s1);
1039 MutableBigInteger t = uOdd ? v : u;
1040 int tsign = uOdd ? -1 : 1;
1041
1042 int lb;
1043 while ((lb = t.getLowestSetBit()) >= 0) {
1044 // steps B3 and B4
1045 t.rightShift(lb);
1046 // step B5
1047 if (tsign > 0)
1048 u = t;
1049 else
1050 v = t;
1051
1052 // Special case one word numbers
1053 if (u.intLen < 2 && v.intLen < 2) {
1054 int x = u.value[u.offset];
1055 int y = v.value[v.offset];
1056 x = binaryGcd(x, y);
1057 r.value[0] = x;
1058 r.intLen = 1;
1059 r.offset = 0;
1060 if (k > 0)
1061 r.leftShift(k);
1062 return r;
1063 }
1064
1065 // step B6
1066 if ((tsign = u.difference(v)) == 0)
1067 break;
1068 t = (tsign >= 0) ? u : v;
1069 }
1070
1071 if (k > 0)
1072 u.leftShift(k);
1073 return u;
1074 }
1075
1076 /**
1077 * Calculate GCD of a and b interpreted as unsigned integers.
1078 */
1079 static int binaryGcd(int a, int b) {
1080 if (b == 0)
1081 return a;
1082 if (a == 0)
1083 return b;
1084
1085 int x;
1086 int aZeros = 0;
1087 while ((x = (int) a & 0xff) == 0) {
1088 a >>>= 8;
1089 aZeros += 8;
1090 }
1091 int y = BigInteger.trailingZeroTable[x];
1092 aZeros += y;
1093 a >>>= y;
1094
1095 int bZeros = 0;
1096 while ((x = (int) b & 0xff) == 0) {
1097 b >>>= 8;
1098 bZeros += 8;
1099 }
1100 y = BigInteger.trailingZeroTable[x];
1101 bZeros += y;
1102 b >>>= y;
1103
1104 int t = (aZeros < bZeros ? aZeros : bZeros);
1105
1106 while (a != b) {
1107 if ((a + 0x80000000) > (b + 0x80000000)) { // a > b as unsigned
1108 a -= b;
1109
1110 while ((x = (int) a & 0xff) == 0)
1111 a >>>= 8;
1112 a >>>= BigInteger.trailingZeroTable[x];
1113 } else {
1114 b -= a;
1115
1116 while ((x = (int) b & 0xff) == 0)
1117 b >>>= 8;
1118 b >>>= BigInteger.trailingZeroTable[x];
1119 }
1120 }
1121 return a << t;
1122 }
1123
1124 /**
1125 * Returns the modInverse of this mod p. This and p are not affected by
1126 * the operation.
1127 */
1128 MutableBigInteger mutableModInverse(MutableBigInteger p) {
1129 // Modulus is odd, use Schroeppel's algorithm
1130 if (p.isOdd())
1131 return modInverse(p);
1132
1133 // Base and modulus are even, throw exception
1134 if (isEven())
1135 throw new ArithmeticException("BigInteger not invertible.");
1136
1137 // Get even part of modulus expressed as a power of 2
1138 int powersOf2 = p.getLowestSetBit();
1139
1140 // Construct odd part of modulus
1141 MutableBigInteger oddMod = new MutableBigInteger(p);
1142 oddMod.rightShift(powersOf2);
1143
1144 if (oddMod.isOne())
1145 return modInverseMP2(powersOf2);
1146
1147 // Calculate 1/a mod oddMod
1148 MutableBigInteger oddPart = modInverse(oddMod);
1149
1150 // Calculate 1/a mod evenMod
1151 MutableBigInteger evenPart = modInverseMP2(powersOf2);
1152
1153 // Combine the results using Chinese Remainder Theorem
1154 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1155 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1156
1157 MutableBigInteger temp1 = new MutableBigInteger();
1158 MutableBigInteger temp2 = new MutableBigInteger();
1159 MutableBigInteger result = new MutableBigInteger();
1160
1161 oddPart.leftShift(powersOf2);
1162 oddPart.multiply(y1, result);
1163
1164 evenPart.multiply(oddMod, temp1);
1165 temp1.multiply(y2, temp2);
1166
1167 result.add(temp2);
1168 result.divide(p, temp1, temp2);
1169 return temp2;
1170 }
1171
1172 /*
1173 * Calculate the multiplicative inverse of this mod 2^k.
1174 */
1175 MutableBigInteger modInverseMP2(int k) {
1176 if (isEven())
1177 throw new ArithmeticException("Non-invertible. (GCD != 1)");
1178
1179 if (k > 64)
1180 return euclidModInverse(k);
1181
1182 int t = inverseMod32(value[offset + intLen - 1]);
1183
1184 if (k < 33) {
1185 t = (k == 32 ? t : t & ((1 << k) - 1));
1186 return new MutableBigInteger(t);
1187 }
1188
1189 long pLong = (value[offset + intLen - 1] & LONG_MASK);
1190 if (intLen > 1)
1191 pLong |= ((long) value[offset + intLen - 2] << 32);
1192 long tLong = t & LONG_MASK;
1193 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
1194 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
1195
1196 MutableBigInteger result = new MutableBigInteger(new int[2]);
1197 result.value[0] = (int) (tLong >>> 32);
1198 result.value[1] = (int) tLong;
1199 result.intLen = 2;
1200 result.normalize();
1201 return result;
1202 }
1203
1204 /*
1205 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
1206 */
1207 static int inverseMod32(int val) {
1208 // Newton's iteration!
1209 int t = val;
1210 t *= 2 - val * t;
1211 t *= 2 - val * t;
1212 t *= 2 - val * t;
1213 t *= 2 - val * t;
1214 return t;
1215 }
1216
1217 /*
1218 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
1219 */
1220 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
1221 // Copy the mod to protect original
1222 return fixup(new MutableBigInteger(1), new MutableBigInteger(
1223 mod), k);
1224 }
1225
1226 /**
1227 * Calculate the multiplicative inverse of this mod mod, where mod is odd.
1228 * This and mod are not changed by the calculation.
1229 *
1230 * This method implements an algorithm due to Richard Schroeppel, that uses
1231 * the same intermediate representation as Montgomery Reduction
1232 * ("Montgomery Form"). The algorithm is described in an unpublished
1233 * manuscript entitled "Fast Modular Reciprocals."
1234 */
1235 private MutableBigInteger modInverse(MutableBigInteger mod) {
1236 MutableBigInteger p = new MutableBigInteger(mod);
1237 MutableBigInteger f = new MutableBigInteger(this );
1238 MutableBigInteger g = new MutableBigInteger(p);
1239 SignedMutableBigInteger c = new SignedMutableBigInteger(1);
1240 SignedMutableBigInteger d = new SignedMutableBigInteger();
1241 MutableBigInteger temp = null;
1242 SignedMutableBigInteger sTemp = null;
1243
1244 int k = 0;
1245 // Right shift f k times until odd, left shift d k times
1246 if (f.isEven()) {
1247 int trailingZeros = f.getLowestSetBit();
1248 f.rightShift(trailingZeros);
1249 d.leftShift(trailingZeros);
1250 k = trailingZeros;
1251 }
1252
1253 // The Almost Inverse Algorithm
1254 while (!f.isOne()) {
1255 // If gcd(f, g) != 1, number is not invertible modulo mod
1256 if (f.isZero())
1257 throw new ArithmeticException(
1258 "BigInteger not invertible.");
1259
1260 // If f < g exchange f, g and c, d
1261 if (f.compare(g) < 0) {
1262 temp = f;
1263 f = g;
1264 g = temp;
1265 sTemp = d;
1266 d = c;
1267 c = sTemp;
1268 }
1269
1270 // If f == g (mod 4)
1271 if (((f.value[f.offset + f.intLen - 1] ^ g.value[g.offset
1272 + g.intLen - 1]) & 3) == 0) {
1273 f.subtract(g);
1274 c.signedSubtract(d);
1275 } else { // If f != g (mod 4)
1276 f.add(g);
1277 c.signedAdd(d);
1278 }
1279
1280 // Right shift f k times until odd, left shift d k times
1281 int trailingZeros = f.getLowestSetBit();
1282 f.rightShift(trailingZeros);
1283 d.leftShift(trailingZeros);
1284 k += trailingZeros;
1285 }
1286
1287 while (c.sign < 0)
1288 c.signedAdd(p);
1289
1290 return fixup(c, p, k);
1291 }
1292
1293 /*
1294 * The Fixup Algorithm
1295 * Calculates X such that X = C * 2^(-k) (mod P)
1296 * Assumes C<P and P is odd.
1297 */
1298 static MutableBigInteger fixup(MutableBigInteger c,
1299 MutableBigInteger p, int k) {
1300 MutableBigInteger temp = new MutableBigInteger();
1301 // Set r to the multiplicative inverse of p mod 2^32
1302 int r = -inverseMod32(p.value[p.offset + p.intLen - 1]);
1303
1304 for (int i = 0, numWords = k >> 5; i < numWords; i++) {
1305 // V = R * c (mod 2^j)
1306 int v = r * c.value[c.offset + c.intLen - 1];
1307 // c = c + (v * p)
1308 p.mul(v, temp);
1309 c.add(temp);
1310 // c = c / 2^j
1311 c.intLen--;
1312 }
1313 int numBits = k & 0x1f;
1314 if (numBits != 0) {
1315 // V = R * c (mod 2^j)
1316 int v = r * c.value[c.offset + c.intLen - 1];
1317 v &= ((1 << numBits) - 1);
1318 // c = c + (v * p)
1319 p.mul(v, temp);
1320 c.add(temp);
1321 // c = c / 2^j
1322 c.rightShift(numBits);
1323 }
1324
1325 // In theory, c may be greater than p at this point (Very rare!)
1326 while (c.compare(p) >= 0)
1327 c.subtract(p);
1328
1329 return c;
1330 }
1331
1332 /**
1333 * Uses the extended Euclidean algorithm to compute the modInverse of base
1334 * mod a modulus that is a power of 2. The modulus is 2^k.
1335 */
1336 MutableBigInteger euclidModInverse(int k) {
1337 MutableBigInteger b = new MutableBigInteger(1);
1338 b.leftShift(k);
1339 MutableBigInteger mod = new MutableBigInteger(b);
1340
1341 MutableBigInteger a = new MutableBigInteger(this );
1342 MutableBigInteger q = new MutableBigInteger();
1343 MutableBigInteger r = new MutableBigInteger();
1344
1345 b.divide(a, q, r);
1346 MutableBigInteger swapper = b;
1347 b = r;
1348 r = swapper;
1349
1350 MutableBigInteger t1 = new MutableBigInteger(q);
1351 MutableBigInteger t0 = new MutableBigInteger(1);
1352 MutableBigInteger temp = new MutableBigInteger();
1353
1354 while (!b.isOne()) {
1355 a.divide(b, q, r);
1356
1357 if (r.intLen == 0)
1358 throw new ArithmeticException(
1359 "BigInteger not invertible.");
1360
1361 swapper = r;
1362 r = a;
1363 a = swapper;
1364
1365 if (q.intLen == 1)
1366 t1.mul(q.value[q.offset], temp);
1367 else
1368 q.multiply(t1, temp);
1369 swapper = q;
1370 q = temp;
1371 temp = swapper;
1372
1373 t0.add(q);
1374
1375 if (a.isOne())
1376 return t0;
1377
1378 b.divide(a, q, r);
1379
1380 if (r.intLen == 0)
1381 throw new ArithmeticException(
1382 "BigInteger not invertible.");
1383
1384 swapper = b;
1385 b = r;
1386 r = swapper;
1387
1388 if (q.intLen == 1)
1389 t0.mul(q.value[q.offset], temp);
1390 else
1391 q.multiply(t0, temp);
1392 swapper = q;
1393 q = temp;
1394 temp = swapper;
1395
1396 t1.add(q);
1397 }
1398 mod.subtract(t1);
1399 return mod;
1400 }
1401
1402 }
|