1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/emul/compact/src/main/java/java/math/MutableBigInteger.java Sat Sep 07 13:51:24 2013 +0200
1.3 @@ -0,0 +1,1477 @@
1.4 +/*
1.5 + * Copyright (c) 1999, 2007, Oracle and/or its affiliates. All rights reserved.
1.6 + * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
1.7 + *
1.8 + * This code is free software; you can redistribute it and/or modify it
1.9 + * under the terms of the GNU General Public License version 2 only, as
1.10 + * published by the Free Software Foundation. Oracle designates this
1.11 + * particular file as subject to the "Classpath" exception as provided
1.12 + * by Oracle in the LICENSE file that accompanied this code.
1.13 + *
1.14 + * This code is distributed in the hope that it will be useful, but WITHOUT
1.15 + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
1.16 + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
1.17 + * version 2 for more details (a copy is included in the LICENSE file that
1.18 + * accompanied this code).
1.19 + *
1.20 + * You should have received a copy of the GNU General Public License version
1.21 + * 2 along with this work; if not, write to the Free Software Foundation,
1.22 + * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
1.23 + *
1.24 + * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
1.25 + * or visit www.oracle.com if you need additional information or have any
1.26 + * questions.
1.27 + */
1.28 +
1.29 +package java.math;
1.30 +
1.31 +/**
1.32 + * A class used to represent multiprecision integers that makes efficient
1.33 + * use of allocated space by allowing a number to occupy only part of
1.34 + * an array so that the arrays do not have to be reallocated as often.
1.35 + * When performing an operation with many iterations the array used to
1.36 + * hold a number is only reallocated when necessary and does not have to
1.37 + * be the same size as the number it represents. A mutable number allows
1.38 + * calculations to occur on the same number without having to create
1.39 + * a new number for every step of the calculation as occurs with
1.40 + * BigIntegers.
1.41 + *
1.42 + * @see BigInteger
1.43 + * @author Michael McCloskey
1.44 + * @since 1.3
1.45 + */
1.46 +
1.47 +import java.util.Arrays;
1.48 +
1.49 +import static java.math.BigInteger.LONG_MASK;
1.50 +import static java.math.BigDecimal.INFLATED;
1.51 +
1.52 +class MutableBigInteger {
1.53 + /**
1.54 + * Holds the magnitude of this MutableBigInteger in big endian order.
1.55 + * The magnitude may start at an offset into the value array, and it may
1.56 + * end before the length of the value array.
1.57 + */
1.58 + int[] value;
1.59 +
1.60 + /**
1.61 + * The number of ints of the value array that are currently used
1.62 + * to hold the magnitude of this MutableBigInteger. The magnitude starts
1.63 + * at an offset and offset + intLen may be less than value.length.
1.64 + */
1.65 + int intLen;
1.66 +
1.67 + /**
1.68 + * The offset into the value array where the magnitude of this
1.69 + * MutableBigInteger begins.
1.70 + */
1.71 + int offset = 0;
1.72 +
1.73 + // Constants
1.74 + /**
1.75 + * MutableBigInteger with one element value array with the value 1. Used by
1.76 + * BigDecimal divideAndRound to increment the quotient. Use this constant
1.77 + * only when the method is not going to modify this object.
1.78 + */
1.79 + static final MutableBigInteger ONE = new MutableBigInteger(1);
1.80 +
1.81 + // Constructors
1.82 +
1.83 + /**
1.84 + * The default constructor. An empty MutableBigInteger is created with
1.85 + * a one word capacity.
1.86 + */
1.87 + MutableBigInteger() {
1.88 + value = new int[1];
1.89 + intLen = 0;
1.90 + }
1.91 +
1.92 + /**
1.93 + * Construct a new MutableBigInteger with a magnitude specified by
1.94 + * the int val.
1.95 + */
1.96 + MutableBigInteger(int val) {
1.97 + value = new int[1];
1.98 + intLen = 1;
1.99 + value[0] = val;
1.100 + }
1.101 +
1.102 + /**
1.103 + * Construct a new MutableBigInteger with the specified value array
1.104 + * up to the length of the array supplied.
1.105 + */
1.106 + MutableBigInteger(int[] val) {
1.107 + value = val;
1.108 + intLen = val.length;
1.109 + }
1.110 +
1.111 + /**
1.112 + * Construct a new MutableBigInteger with a magnitude equal to the
1.113 + * specified BigInteger.
1.114 + */
1.115 + MutableBigInteger(BigInteger b) {
1.116 + intLen = b.mag.length;
1.117 + value = Arrays.copyOf(b.mag, intLen);
1.118 + }
1.119 +
1.120 + /**
1.121 + * Construct a new MutableBigInteger with a magnitude equal to the
1.122 + * specified MutableBigInteger.
1.123 + */
1.124 + MutableBigInteger(MutableBigInteger val) {
1.125 + intLen = val.intLen;
1.126 + value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
1.127 + }
1.128 +
1.129 + /**
1.130 + * Internal helper method to return the magnitude array. The caller is not
1.131 + * supposed to modify the returned array.
1.132 + */
1.133 + private int[] getMagnitudeArray() {
1.134 + if (offset > 0 || value.length != intLen)
1.135 + return Arrays.copyOfRange(value, offset, offset + intLen);
1.136 + return value;
1.137 + }
1.138 +
1.139 + /**
1.140 + * Convert this MutableBigInteger to a long value. The caller has to make
1.141 + * sure this MutableBigInteger can be fit into long.
1.142 + */
1.143 + private long toLong() {
1.144 + assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
1.145 + if (intLen == 0)
1.146 + return 0;
1.147 + long d = value[offset] & LONG_MASK;
1.148 + return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
1.149 + }
1.150 +
1.151 + /**
1.152 + * Convert this MutableBigInteger to a BigInteger object.
1.153 + */
1.154 + BigInteger toBigInteger(int sign) {
1.155 + if (intLen == 0 || sign == 0)
1.156 + return BigInteger.ZERO;
1.157 + return new BigInteger(getMagnitudeArray(), sign);
1.158 + }
1.159 +
1.160 + /**
1.161 + * Convert this MutableBigInteger to BigDecimal object with the specified sign
1.162 + * and scale.
1.163 + */
1.164 + BigDecimal toBigDecimal(int sign, int scale) {
1.165 + if (intLen == 0 || sign == 0)
1.166 + return BigDecimal.valueOf(0, scale);
1.167 + int[] mag = getMagnitudeArray();
1.168 + int len = mag.length;
1.169 + int d = mag[0];
1.170 + // If this MutableBigInteger can't be fit into long, we need to
1.171 + // make a BigInteger object for the resultant BigDecimal object.
1.172 + if (len > 2 || (d < 0 && len == 2))
1.173 + return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
1.174 + long v = (len == 2) ?
1.175 + ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
1.176 + d & LONG_MASK;
1.177 + return new BigDecimal(null, sign == -1 ? -v : v, scale, 0);
1.178 + }
1.179 +
1.180 + /**
1.181 + * Clear out a MutableBigInteger for reuse.
1.182 + */
1.183 + void clear() {
1.184 + offset = intLen = 0;
1.185 + for (int index=0, n=value.length; index < n; index++)
1.186 + value[index] = 0;
1.187 + }
1.188 +
1.189 + /**
1.190 + * Set a MutableBigInteger to zero, removing its offset.
1.191 + */
1.192 + void reset() {
1.193 + offset = intLen = 0;
1.194 + }
1.195 +
1.196 + /**
1.197 + * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
1.198 + * as this MutableBigInteger is numerically less than, equal to, or
1.199 + * greater than <tt>b</tt>.
1.200 + */
1.201 + final int compare(MutableBigInteger b) {
1.202 + int blen = b.intLen;
1.203 + if (intLen < blen)
1.204 + return -1;
1.205 + if (intLen > blen)
1.206 + return 1;
1.207 +
1.208 + // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
1.209 + // comparison.
1.210 + int[] bval = b.value;
1.211 + for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
1.212 + int b1 = value[i] + 0x80000000;
1.213 + int b2 = bval[j] + 0x80000000;
1.214 + if (b1 < b2)
1.215 + return -1;
1.216 + if (b1 > b2)
1.217 + return 1;
1.218 + }
1.219 + return 0;
1.220 + }
1.221 +
1.222 + /**
1.223 + * Compare this against half of a MutableBigInteger object (Needed for
1.224 + * remainder tests).
1.225 + * Assumes no leading unnecessary zeros, which holds for results
1.226 + * from divide().
1.227 + */
1.228 + final int compareHalf(MutableBigInteger b) {
1.229 + int blen = b.intLen;
1.230 + int len = intLen;
1.231 + if (len <= 0)
1.232 + return blen <=0 ? 0 : -1;
1.233 + if (len > blen)
1.234 + return 1;
1.235 + if (len < blen - 1)
1.236 + return -1;
1.237 + int[] bval = b.value;
1.238 + int bstart = 0;
1.239 + int carry = 0;
1.240 + // Only 2 cases left:len == blen or len == blen - 1
1.241 + if (len != blen) { // len == blen - 1
1.242 + if (bval[bstart] == 1) {
1.243 + ++bstart;
1.244 + carry = 0x80000000;
1.245 + } else
1.246 + return -1;
1.247 + }
1.248 + // compare values with right-shifted values of b,
1.249 + // carrying shifted-out bits across words
1.250 + int[] val = value;
1.251 + for (int i = offset, j = bstart; i < len + offset;) {
1.252 + int bv = bval[j++];
1.253 + long hb = ((bv >>> 1) + carry) & LONG_MASK;
1.254 + long v = val[i++] & LONG_MASK;
1.255 + if (v != hb)
1.256 + return v < hb ? -1 : 1;
1.257 + carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
1.258 + }
1.259 + return carry == 0? 0 : -1;
1.260 + }
1.261 +
1.262 + /**
1.263 + * Return the index of the lowest set bit in this MutableBigInteger. If the
1.264 + * magnitude of this MutableBigInteger is zero, -1 is returned.
1.265 + */
1.266 + private final int getLowestSetBit() {
1.267 + if (intLen == 0)
1.268 + return -1;
1.269 + int j, b;
1.270 + for (j=intLen-1; (j>0) && (value[j+offset]==0); j--)
1.271 + ;
1.272 + b = value[j+offset];
1.273 + if (b==0)
1.274 + return -1;
1.275 + return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
1.276 + }
1.277 +
1.278 + /**
1.279 + * Return the int in use in this MutableBigInteger at the specified
1.280 + * index. This method is not used because it is not inlined on all
1.281 + * platforms.
1.282 + */
1.283 + private final int getInt(int index) {
1.284 + return value[offset+index];
1.285 + }
1.286 +
1.287 + /**
1.288 + * Return a long which is equal to the unsigned value of the int in
1.289 + * use in this MutableBigInteger at the specified index. This method is
1.290 + * not used because it is not inlined on all platforms.
1.291 + */
1.292 + private final long getLong(int index) {
1.293 + return value[offset+index] & LONG_MASK;
1.294 + }
1.295 +
1.296 + /**
1.297 + * Ensure that the MutableBigInteger is in normal form, specifically
1.298 + * making sure that there are no leading zeros, and that if the
1.299 + * magnitude is zero, then intLen is zero.
1.300 + */
1.301 + final void normalize() {
1.302 + if (intLen == 0) {
1.303 + offset = 0;
1.304 + return;
1.305 + }
1.306 +
1.307 + int index = offset;
1.308 + if (value[index] != 0)
1.309 + return;
1.310 +
1.311 + int indexBound = index+intLen;
1.312 + do {
1.313 + index++;
1.314 + } while(index < indexBound && value[index]==0);
1.315 +
1.316 + int numZeros = index - offset;
1.317 + intLen -= numZeros;
1.318 + offset = (intLen==0 ? 0 : offset+numZeros);
1.319 + }
1.320 +
1.321 + /**
1.322 + * If this MutableBigInteger cannot hold len words, increase the size
1.323 + * of the value array to len words.
1.324 + */
1.325 + private final void ensureCapacity(int len) {
1.326 + if (value.length < len) {
1.327 + value = new int[len];
1.328 + offset = 0;
1.329 + intLen = len;
1.330 + }
1.331 + }
1.332 +
1.333 + /**
1.334 + * Convert this MutableBigInteger into an int array with no leading
1.335 + * zeros, of a length that is equal to this MutableBigInteger's intLen.
1.336 + */
1.337 + int[] toIntArray() {
1.338 + int[] result = new int[intLen];
1.339 + for(int i=0; i<intLen; i++)
1.340 + result[i] = value[offset+i];
1.341 + return result;
1.342 + }
1.343 +
1.344 + /**
1.345 + * Sets the int at index+offset in this MutableBigInteger to val.
1.346 + * This does not get inlined on all platforms so it is not used
1.347 + * as often as originally intended.
1.348 + */
1.349 + void setInt(int index, int val) {
1.350 + value[offset + index] = val;
1.351 + }
1.352 +
1.353 + /**
1.354 + * Sets this MutableBigInteger's value array to the specified array.
1.355 + * The intLen is set to the specified length.
1.356 + */
1.357 + void setValue(int[] val, int length) {
1.358 + value = val;
1.359 + intLen = length;
1.360 + offset = 0;
1.361 + }
1.362 +
1.363 + /**
1.364 + * Sets this MutableBigInteger's value array to a copy of the specified
1.365 + * array. The intLen is set to the length of the new array.
1.366 + */
1.367 + void copyValue(MutableBigInteger src) {
1.368 + int len = src.intLen;
1.369 + if (value.length < len)
1.370 + value = new int[len];
1.371 + System.arraycopy(src.value, src.offset, value, 0, len);
1.372 + intLen = len;
1.373 + offset = 0;
1.374 + }
1.375 +
1.376 + /**
1.377 + * Sets this MutableBigInteger's value array to a copy of the specified
1.378 + * array. The intLen is set to the length of the specified array.
1.379 + */
1.380 + void copyValue(int[] val) {
1.381 + int len = val.length;
1.382 + if (value.length < len)
1.383 + value = new int[len];
1.384 + System.arraycopy(val, 0, value, 0, len);
1.385 + intLen = len;
1.386 + offset = 0;
1.387 + }
1.388 +
1.389 + /**
1.390 + * Returns true iff this MutableBigInteger has a value of one.
1.391 + */
1.392 + boolean isOne() {
1.393 + return (intLen == 1) && (value[offset] == 1);
1.394 + }
1.395 +
1.396 + /**
1.397 + * Returns true iff this MutableBigInteger has a value of zero.
1.398 + */
1.399 + boolean isZero() {
1.400 + return (intLen == 0);
1.401 + }
1.402 +
1.403 + /**
1.404 + * Returns true iff this MutableBigInteger is even.
1.405 + */
1.406 + boolean isEven() {
1.407 + return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
1.408 + }
1.409 +
1.410 + /**
1.411 + * Returns true iff this MutableBigInteger is odd.
1.412 + */
1.413 + boolean isOdd() {
1.414 + return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
1.415 + }
1.416 +
1.417 + /**
1.418 + * Returns true iff this MutableBigInteger is in normal form. A
1.419 + * MutableBigInteger is in normal form if it has no leading zeros
1.420 + * after the offset, and intLen + offset <= value.length.
1.421 + */
1.422 + boolean isNormal() {
1.423 + if (intLen + offset > value.length)
1.424 + return false;
1.425 + if (intLen ==0)
1.426 + return true;
1.427 + return (value[offset] != 0);
1.428 + }
1.429 +
1.430 + /**
1.431 + * Returns a String representation of this MutableBigInteger in radix 10.
1.432 + */
1.433 + public String toString() {
1.434 + BigInteger b = toBigInteger(1);
1.435 + return b.toString();
1.436 + }
1.437 +
1.438 + /**
1.439 + * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
1.440 + * in normal form.
1.441 + */
1.442 + void rightShift(int n) {
1.443 + if (intLen == 0)
1.444 + return;
1.445 + int nInts = n >>> 5;
1.446 + int nBits = n & 0x1F;
1.447 + this.intLen -= nInts;
1.448 + if (nBits == 0)
1.449 + return;
1.450 + int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
1.451 + if (nBits >= bitsInHighWord) {
1.452 + this.primitiveLeftShift(32 - nBits);
1.453 + this.intLen--;
1.454 + } else {
1.455 + primitiveRightShift(nBits);
1.456 + }
1.457 + }
1.458 +
1.459 + /**
1.460 + * Left shift this MutableBigInteger n bits.
1.461 + */
1.462 + void leftShift(int n) {
1.463 + /*
1.464 + * If there is enough storage space in this MutableBigInteger already
1.465 + * the available space will be used. Space to the right of the used
1.466 + * ints in the value array is faster to utilize, so the extra space
1.467 + * will be taken from the right if possible.
1.468 + */
1.469 + if (intLen == 0)
1.470 + return;
1.471 + int nInts = n >>> 5;
1.472 + int nBits = n&0x1F;
1.473 + int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
1.474 +
1.475 + // If shift can be done without moving words, do so
1.476 + if (n <= (32-bitsInHighWord)) {
1.477 + primitiveLeftShift(nBits);
1.478 + return;
1.479 + }
1.480 +
1.481 + int newLen = intLen + nInts +1;
1.482 + if (nBits <= (32-bitsInHighWord))
1.483 + newLen--;
1.484 + if (value.length < newLen) {
1.485 + // The array must grow
1.486 + int[] result = new int[newLen];
1.487 + for (int i=0; i<intLen; i++)
1.488 + result[i] = value[offset+i];
1.489 + setValue(result, newLen);
1.490 + } else if (value.length - offset >= newLen) {
1.491 + // Use space on right
1.492 + for(int i=0; i<newLen - intLen; i++)
1.493 + value[offset+intLen+i] = 0;
1.494 + } else {
1.495 + // Must use space on left
1.496 + for (int i=0; i<intLen; i++)
1.497 + value[i] = value[offset+i];
1.498 + for (int i=intLen; i<newLen; i++)
1.499 + value[i] = 0;
1.500 + offset = 0;
1.501 + }
1.502 + intLen = newLen;
1.503 + if (nBits == 0)
1.504 + return;
1.505 + if (nBits <= (32-bitsInHighWord))
1.506 + primitiveLeftShift(nBits);
1.507 + else
1.508 + primitiveRightShift(32 -nBits);
1.509 + }
1.510 +
1.511 + /**
1.512 + * A primitive used for division. This method adds in one multiple of the
1.513 + * divisor a back to the dividend result at a specified offset. It is used
1.514 + * when qhat was estimated too large, and must be adjusted.
1.515 + */
1.516 + private int divadd(int[] a, int[] result, int offset) {
1.517 + long carry = 0;
1.518 +
1.519 + for (int j=a.length-1; j >= 0; j--) {
1.520 + long sum = (a[j] & LONG_MASK) +
1.521 + (result[j+offset] & LONG_MASK) + carry;
1.522 + result[j+offset] = (int)sum;
1.523 + carry = sum >>> 32;
1.524 + }
1.525 + return (int)carry;
1.526 + }
1.527 +
1.528 + /**
1.529 + * This method is used for division. It multiplies an n word input a by one
1.530 + * word input x, and subtracts the n word product from q. This is needed
1.531 + * when subtracting qhat*divisor from dividend.
1.532 + */
1.533 + private int mulsub(int[] q, int[] a, int x, int len, int offset) {
1.534 + long xLong = x & LONG_MASK;
1.535 + long carry = 0;
1.536 + offset += len;
1.537 +
1.538 + for (int j=len-1; j >= 0; j--) {
1.539 + long product = (a[j] & LONG_MASK) * xLong + carry;
1.540 + long difference = q[offset] - product;
1.541 + q[offset--] = (int)difference;
1.542 + carry = (product >>> 32)
1.543 + + (((difference & LONG_MASK) >
1.544 + (((~(int)product) & LONG_MASK))) ? 1:0);
1.545 + }
1.546 + return (int)carry;
1.547 + }
1.548 +
1.549 + /**
1.550 + * Right shift this MutableBigInteger n bits, where n is
1.551 + * less than 32.
1.552 + * Assumes that intLen > 0, n > 0 for speed
1.553 + */
1.554 + private final void primitiveRightShift(int n) {
1.555 + int[] val = value;
1.556 + int n2 = 32 - n;
1.557 + for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
1.558 + int b = c;
1.559 + c = val[i-1];
1.560 + val[i] = (c << n2) | (b >>> n);
1.561 + }
1.562 + val[offset] >>>= n;
1.563 + }
1.564 +
1.565 + /**
1.566 + * Left shift this MutableBigInteger n bits, where n is
1.567 + * less than 32.
1.568 + * Assumes that intLen > 0, n > 0 for speed
1.569 + */
1.570 + private final void primitiveLeftShift(int n) {
1.571 + int[] val = value;
1.572 + int n2 = 32 - n;
1.573 + for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
1.574 + int b = c;
1.575 + c = val[i+1];
1.576 + val[i] = (b << n) | (c >>> n2);
1.577 + }
1.578 + val[offset+intLen-1] <<= n;
1.579 + }
1.580 +
1.581 + /**
1.582 + * Adds the contents of two MutableBigInteger objects.The result
1.583 + * is placed within this MutableBigInteger.
1.584 + * The contents of the addend are not changed.
1.585 + */
1.586 + void add(MutableBigInteger addend) {
1.587 + int x = intLen;
1.588 + int y = addend.intLen;
1.589 + int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
1.590 + int[] result = (value.length < resultLen ? new int[resultLen] : value);
1.591 +
1.592 + int rstart = result.length-1;
1.593 + long sum;
1.594 + long carry = 0;
1.595 +
1.596 + // Add common parts of both numbers
1.597 + while(x>0 && y>0) {
1.598 + x--; y--;
1.599 + sum = (value[x+offset] & LONG_MASK) +
1.600 + (addend.value[y+addend.offset] & LONG_MASK) + carry;
1.601 + result[rstart--] = (int)sum;
1.602 + carry = sum >>> 32;
1.603 + }
1.604 +
1.605 + // Add remainder of the longer number
1.606 + while(x>0) {
1.607 + x--;
1.608 + if (carry == 0 && result == value && rstart == (x + offset))
1.609 + return;
1.610 + sum = (value[x+offset] & LONG_MASK) + carry;
1.611 + result[rstart--] = (int)sum;
1.612 + carry = sum >>> 32;
1.613 + }
1.614 + while(y>0) {
1.615 + y--;
1.616 + sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
1.617 + result[rstart--] = (int)sum;
1.618 + carry = sum >>> 32;
1.619 + }
1.620 +
1.621 + if (carry > 0) { // Result must grow in length
1.622 + resultLen++;
1.623 + if (result.length < resultLen) {
1.624 + int temp[] = new int[resultLen];
1.625 + // Result one word longer from carry-out; copy low-order
1.626 + // bits into new result.
1.627 + System.arraycopy(result, 0, temp, 1, result.length);
1.628 + temp[0] = 1;
1.629 + result = temp;
1.630 + } else {
1.631 + result[rstart--] = 1;
1.632 + }
1.633 + }
1.634 +
1.635 + value = result;
1.636 + intLen = resultLen;
1.637 + offset = result.length - resultLen;
1.638 + }
1.639 +
1.640 +
1.641 + /**
1.642 + * Subtracts the smaller of this and b from the larger and places the
1.643 + * result into this MutableBigInteger.
1.644 + */
1.645 + int subtract(MutableBigInteger b) {
1.646 + MutableBigInteger a = this;
1.647 +
1.648 + int[] result = value;
1.649 + int sign = a.compare(b);
1.650 +
1.651 + if (sign == 0) {
1.652 + reset();
1.653 + return 0;
1.654 + }
1.655 + if (sign < 0) {
1.656 + MutableBigInteger tmp = a;
1.657 + a = b;
1.658 + b = tmp;
1.659 + }
1.660 +
1.661 + int resultLen = a.intLen;
1.662 + if (result.length < resultLen)
1.663 + result = new int[resultLen];
1.664 +
1.665 + long diff = 0;
1.666 + int x = a.intLen;
1.667 + int y = b.intLen;
1.668 + int rstart = result.length - 1;
1.669 +
1.670 + // Subtract common parts of both numbers
1.671 + while (y>0) {
1.672 + x--; y--;
1.673 +
1.674 + diff = (a.value[x+a.offset] & LONG_MASK) -
1.675 + (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
1.676 + result[rstart--] = (int)diff;
1.677 + }
1.678 + // Subtract remainder of longer number
1.679 + while (x>0) {
1.680 + x--;
1.681 + diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
1.682 + result[rstart--] = (int)diff;
1.683 + }
1.684 +
1.685 + value = result;
1.686 + intLen = resultLen;
1.687 + offset = value.length - resultLen;
1.688 + normalize();
1.689 + return sign;
1.690 + }
1.691 +
1.692 + /**
1.693 + * Subtracts the smaller of a and b from the larger and places the result
1.694 + * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
1.695 + * operation was performed.
1.696 + */
1.697 + private int difference(MutableBigInteger b) {
1.698 + MutableBigInteger a = this;
1.699 + int sign = a.compare(b);
1.700 + if (sign ==0)
1.701 + return 0;
1.702 + if (sign < 0) {
1.703 + MutableBigInteger tmp = a;
1.704 + a = b;
1.705 + b = tmp;
1.706 + }
1.707 +
1.708 + long diff = 0;
1.709 + int x = a.intLen;
1.710 + int y = b.intLen;
1.711 +
1.712 + // Subtract common parts of both numbers
1.713 + while (y>0) {
1.714 + x--; y--;
1.715 + diff = (a.value[a.offset+ x] & LONG_MASK) -
1.716 + (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
1.717 + a.value[a.offset+x] = (int)diff;
1.718 + }
1.719 + // Subtract remainder of longer number
1.720 + while (x>0) {
1.721 + x--;
1.722 + diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
1.723 + a.value[a.offset+x] = (int)diff;
1.724 + }
1.725 +
1.726 + a.normalize();
1.727 + return sign;
1.728 + }
1.729 +
1.730 + /**
1.731 + * Multiply the contents of two MutableBigInteger objects. The result is
1.732 + * placed into MutableBigInteger z. The contents of y are not changed.
1.733 + */
1.734 + void multiply(MutableBigInteger y, MutableBigInteger z) {
1.735 + int xLen = intLen;
1.736 + int yLen = y.intLen;
1.737 + int newLen = xLen + yLen;
1.738 +
1.739 + // Put z into an appropriate state to receive product
1.740 + if (z.value.length < newLen)
1.741 + z.value = new int[newLen];
1.742 + z.offset = 0;
1.743 + z.intLen = newLen;
1.744 +
1.745 + // The first iteration is hoisted out of the loop to avoid extra add
1.746 + long carry = 0;
1.747 + for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
1.748 + long product = (y.value[j+y.offset] & LONG_MASK) *
1.749 + (value[xLen-1+offset] & LONG_MASK) + carry;
1.750 + z.value[k] = (int)product;
1.751 + carry = product >>> 32;
1.752 + }
1.753 + z.value[xLen-1] = (int)carry;
1.754 +
1.755 + // Perform the multiplication word by word
1.756 + for (int i = xLen-2; i >= 0; i--) {
1.757 + carry = 0;
1.758 + for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
1.759 + long product = (y.value[j+y.offset] & LONG_MASK) *
1.760 + (value[i+offset] & LONG_MASK) +
1.761 + (z.value[k] & LONG_MASK) + carry;
1.762 + z.value[k] = (int)product;
1.763 + carry = product >>> 32;
1.764 + }
1.765 + z.value[i] = (int)carry;
1.766 + }
1.767 +
1.768 + // Remove leading zeros from product
1.769 + z.normalize();
1.770 + }
1.771 +
1.772 + /**
1.773 + * Multiply the contents of this MutableBigInteger by the word y. The
1.774 + * result is placed into z.
1.775 + */
1.776 + void mul(int y, MutableBigInteger z) {
1.777 + if (y == 1) {
1.778 + z.copyValue(this);
1.779 + return;
1.780 + }
1.781 +
1.782 + if (y == 0) {
1.783 + z.clear();
1.784 + return;
1.785 + }
1.786 +
1.787 + // Perform the multiplication word by word
1.788 + long ylong = y & LONG_MASK;
1.789 + int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
1.790 + : z.value);
1.791 + long carry = 0;
1.792 + for (int i = intLen-1; i >= 0; i--) {
1.793 + long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1.794 + zval[i+1] = (int)product;
1.795 + carry = product >>> 32;
1.796 + }
1.797 +
1.798 + if (carry == 0) {
1.799 + z.offset = 1;
1.800 + z.intLen = intLen;
1.801 + } else {
1.802 + z.offset = 0;
1.803 + z.intLen = intLen + 1;
1.804 + zval[0] = (int)carry;
1.805 + }
1.806 + z.value = zval;
1.807 + }
1.808 +
1.809 + /**
1.810 + * This method is used for division of an n word dividend by a one word
1.811 + * divisor. The quotient is placed into quotient. The one word divisor is
1.812 + * specified by divisor.
1.813 + *
1.814 + * @return the remainder of the division is returned.
1.815 + *
1.816 + */
1.817 + int divideOneWord(int divisor, MutableBigInteger quotient) {
1.818 + long divisorLong = divisor & LONG_MASK;
1.819 +
1.820 + // Special case of one word dividend
1.821 + if (intLen == 1) {
1.822 + long dividendValue = value[offset] & LONG_MASK;
1.823 + int q = (int) (dividendValue / divisorLong);
1.824 + int r = (int) (dividendValue - q * divisorLong);
1.825 + quotient.value[0] = q;
1.826 + quotient.intLen = (q == 0) ? 0 : 1;
1.827 + quotient.offset = 0;
1.828 + return r;
1.829 + }
1.830 +
1.831 + if (quotient.value.length < intLen)
1.832 + quotient.value = new int[intLen];
1.833 + quotient.offset = 0;
1.834 + quotient.intLen = intLen;
1.835 +
1.836 + // Normalize the divisor
1.837 + int shift = Integer.numberOfLeadingZeros(divisor);
1.838 +
1.839 + int rem = value[offset];
1.840 + long remLong = rem & LONG_MASK;
1.841 + if (remLong < divisorLong) {
1.842 + quotient.value[0] = 0;
1.843 + } else {
1.844 + quotient.value[0] = (int)(remLong / divisorLong);
1.845 + rem = (int) (remLong - (quotient.value[0] * divisorLong));
1.846 + remLong = rem & LONG_MASK;
1.847 + }
1.848 +
1.849 + int xlen = intLen;
1.850 + int[] qWord = new int[2];
1.851 + while (--xlen > 0) {
1.852 + long dividendEstimate = (remLong<<32) |
1.853 + (value[offset + intLen - xlen] & LONG_MASK);
1.854 + if (dividendEstimate >= 0) {
1.855 + qWord[0] = (int) (dividendEstimate / divisorLong);
1.856 + qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong);
1.857 + } else {
1.858 + divWord(qWord, dividendEstimate, divisor);
1.859 + }
1.860 + quotient.value[intLen - xlen] = qWord[0];
1.861 + rem = qWord[1];
1.862 + remLong = rem & LONG_MASK;
1.863 + }
1.864 +
1.865 + quotient.normalize();
1.866 + // Unnormalize
1.867 + if (shift > 0)
1.868 + return rem % divisor;
1.869 + else
1.870 + return rem;
1.871 + }
1.872 +
1.873 + /**
1.874 + * Calculates the quotient of this div b and places the quotient in the
1.875 + * provided MutableBigInteger objects and the remainder object is returned.
1.876 + *
1.877 + * Uses Algorithm D in Knuth section 4.3.1.
1.878 + * Many optimizations to that algorithm have been adapted from the Colin
1.879 + * Plumb C library.
1.880 + * It special cases one word divisors for speed. The content of b is not
1.881 + * changed.
1.882 + *
1.883 + */
1.884 + MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1.885 + if (b.intLen == 0)
1.886 + throw new ArithmeticException("BigInteger divide by zero");
1.887 +
1.888 + // Dividend is zero
1.889 + if (intLen == 0) {
1.890 + quotient.intLen = quotient.offset;
1.891 + return new MutableBigInteger();
1.892 + }
1.893 +
1.894 + int cmp = compare(b);
1.895 + // Dividend less than divisor
1.896 + if (cmp < 0) {
1.897 + quotient.intLen = quotient.offset = 0;
1.898 + return new MutableBigInteger(this);
1.899 + }
1.900 + // Dividend equal to divisor
1.901 + if (cmp == 0) {
1.902 + quotient.value[0] = quotient.intLen = 1;
1.903 + quotient.offset = 0;
1.904 + return new MutableBigInteger();
1.905 + }
1.906 +
1.907 + quotient.clear();
1.908 + // Special case one word divisor
1.909 + if (b.intLen == 1) {
1.910 + int r = divideOneWord(b.value[b.offset], quotient);
1.911 + if (r == 0)
1.912 + return new MutableBigInteger();
1.913 + return new MutableBigInteger(r);
1.914 + }
1.915 +
1.916 + // Copy divisor value to protect divisor
1.917 + int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen);
1.918 + return divideMagnitude(div, quotient);
1.919 + }
1.920 +
1.921 + /**
1.922 + * Internally used to calculate the quotient of this div v and places the
1.923 + * quotient in the provided MutableBigInteger object and the remainder is
1.924 + * returned.
1.925 + *
1.926 + * @return the remainder of the division will be returned.
1.927 + */
1.928 + long divide(long v, MutableBigInteger quotient) {
1.929 + if (v == 0)
1.930 + throw new ArithmeticException("BigInteger divide by zero");
1.931 +
1.932 + // Dividend is zero
1.933 + if (intLen == 0) {
1.934 + quotient.intLen = quotient.offset = 0;
1.935 + return 0;
1.936 + }
1.937 + if (v < 0)
1.938 + v = -v;
1.939 +
1.940 + int d = (int)(v >>> 32);
1.941 + quotient.clear();
1.942 + // Special case on word divisor
1.943 + if (d == 0)
1.944 + return divideOneWord((int)v, quotient) & LONG_MASK;
1.945 + else {
1.946 + int[] div = new int[]{ d, (int)(v & LONG_MASK) };
1.947 + return divideMagnitude(div, quotient).toLong();
1.948 + }
1.949 + }
1.950 +
1.951 + /**
1.952 + * Divide this MutableBigInteger by the divisor represented by its magnitude
1.953 + * array. The quotient will be placed into the provided quotient object &
1.954 + * the remainder object is returned.
1.955 + */
1.956 + private MutableBigInteger divideMagnitude(int[] divisor,
1.957 + MutableBigInteger quotient) {
1.958 +
1.959 + // Remainder starts as dividend with space for a leading zero
1.960 + MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1.961 + System.arraycopy(value, offset, rem.value, 1, intLen);
1.962 + rem.intLen = intLen;
1.963 + rem.offset = 1;
1.964 +
1.965 + int nlen = rem.intLen;
1.966 +
1.967 + // Set the quotient size
1.968 + int dlen = divisor.length;
1.969 + int limit = nlen - dlen + 1;
1.970 + if (quotient.value.length < limit) {
1.971 + quotient.value = new int[limit];
1.972 + quotient.offset = 0;
1.973 + }
1.974 + quotient.intLen = limit;
1.975 + int[] q = quotient.value;
1.976 +
1.977 + // D1 normalize the divisor
1.978 + int shift = Integer.numberOfLeadingZeros(divisor[0]);
1.979 + if (shift > 0) {
1.980 + // First shift will not grow array
1.981 + BigInteger.primitiveLeftShift(divisor, dlen, shift);
1.982 + // But this one might
1.983 + rem.leftShift(shift);
1.984 + }
1.985 +
1.986 + // Must insert leading 0 in rem if its length did not change
1.987 + if (rem.intLen == nlen) {
1.988 + rem.offset = 0;
1.989 + rem.value[0] = 0;
1.990 + rem.intLen++;
1.991 + }
1.992 +
1.993 + int dh = divisor[0];
1.994 + long dhLong = dh & LONG_MASK;
1.995 + int dl = divisor[1];
1.996 + int[] qWord = new int[2];
1.997 +
1.998 + // D2 Initialize j
1.999 + for(int j=0; j<limit; j++) {
1.1000 + // D3 Calculate qhat
1.1001 + // estimate qhat
1.1002 + int qhat = 0;
1.1003 + int qrem = 0;
1.1004 + boolean skipCorrection = false;
1.1005 + int nh = rem.value[j+rem.offset];
1.1006 + int nh2 = nh + 0x80000000;
1.1007 + int nm = rem.value[j+1+rem.offset];
1.1008 +
1.1009 + if (nh == dh) {
1.1010 + qhat = ~0;
1.1011 + qrem = nh + nm;
1.1012 + skipCorrection = qrem + 0x80000000 < nh2;
1.1013 + } else {
1.1014 + long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1.1015 + if (nChunk >= 0) {
1.1016 + qhat = (int) (nChunk / dhLong);
1.1017 + qrem = (int) (nChunk - (qhat * dhLong));
1.1018 + } else {
1.1019 + divWord(qWord, nChunk, dh);
1.1020 + qhat = qWord[0];
1.1021 + qrem = qWord[1];
1.1022 + }
1.1023 + }
1.1024 +
1.1025 + if (qhat == 0)
1.1026 + continue;
1.1027 +
1.1028 + if (!skipCorrection) { // Correct qhat
1.1029 + long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1.1030 + long rs = ((qrem & LONG_MASK) << 32) | nl;
1.1031 + long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1.1032 +
1.1033 + if (unsignedLongCompare(estProduct, rs)) {
1.1034 + qhat--;
1.1035 + qrem = (int)((qrem & LONG_MASK) + dhLong);
1.1036 + if ((qrem & LONG_MASK) >= dhLong) {
1.1037 + estProduct -= (dl & LONG_MASK);
1.1038 + rs = ((qrem & LONG_MASK) << 32) | nl;
1.1039 + if (unsignedLongCompare(estProduct, rs))
1.1040 + qhat--;
1.1041 + }
1.1042 + }
1.1043 + }
1.1044 +
1.1045 + // D4 Multiply and subtract
1.1046 + rem.value[j+rem.offset] = 0;
1.1047 + int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1.1048 +
1.1049 + // D5 Test remainder
1.1050 + if (borrow + 0x80000000 > nh2) {
1.1051 + // D6 Add back
1.1052 + divadd(divisor, rem.value, j+1+rem.offset);
1.1053 + qhat--;
1.1054 + }
1.1055 +
1.1056 + // Store the quotient digit
1.1057 + q[j] = qhat;
1.1058 + } // D7 loop on j
1.1059 +
1.1060 + // D8 Unnormalize
1.1061 + if (shift > 0)
1.1062 + rem.rightShift(shift);
1.1063 +
1.1064 + quotient.normalize();
1.1065 + rem.normalize();
1.1066 + return rem;
1.1067 + }
1.1068 +
1.1069 + /**
1.1070 + * Compare two longs as if they were unsigned.
1.1071 + * Returns true iff one is bigger than two.
1.1072 + */
1.1073 + private boolean unsignedLongCompare(long one, long two) {
1.1074 + return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1.1075 + }
1.1076 +
1.1077 + /**
1.1078 + * This method divides a long quantity by an int to estimate
1.1079 + * qhat for two multi precision numbers. It is used when
1.1080 + * the signed value of n is less than zero.
1.1081 + */
1.1082 + private void divWord(int[] result, long n, int d) {
1.1083 + long dLong = d & LONG_MASK;
1.1084 +
1.1085 + if (dLong == 1) {
1.1086 + result[0] = (int)n;
1.1087 + result[1] = 0;
1.1088 + return;
1.1089 + }
1.1090 +
1.1091 + // Approximate the quotient and remainder
1.1092 + long q = (n >>> 1) / (dLong >>> 1);
1.1093 + long r = n - q*dLong;
1.1094 +
1.1095 + // Correct the approximation
1.1096 + while (r < 0) {
1.1097 + r += dLong;
1.1098 + q--;
1.1099 + }
1.1100 + while (r >= dLong) {
1.1101 + r -= dLong;
1.1102 + q++;
1.1103 + }
1.1104 +
1.1105 + // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1.1106 + result[0] = (int)q;
1.1107 + result[1] = (int)r;
1.1108 + }
1.1109 +
1.1110 + /**
1.1111 + * Calculate GCD of this and b. This and b are changed by the computation.
1.1112 + */
1.1113 + MutableBigInteger hybridGCD(MutableBigInteger b) {
1.1114 + // Use Euclid's algorithm until the numbers are approximately the
1.1115 + // same length, then use the binary GCD algorithm to find the GCD.
1.1116 + MutableBigInteger a = this;
1.1117 + MutableBigInteger q = new MutableBigInteger();
1.1118 +
1.1119 + while (b.intLen != 0) {
1.1120 + if (Math.abs(a.intLen - b.intLen) < 2)
1.1121 + return a.binaryGCD(b);
1.1122 +
1.1123 + MutableBigInteger r = a.divide(b, q);
1.1124 + a = b;
1.1125 + b = r;
1.1126 + }
1.1127 + return a;
1.1128 + }
1.1129 +
1.1130 + /**
1.1131 + * Calculate GCD of this and v.
1.1132 + * Assumes that this and v are not zero.
1.1133 + */
1.1134 + private MutableBigInteger binaryGCD(MutableBigInteger v) {
1.1135 + // Algorithm B from Knuth section 4.5.2
1.1136 + MutableBigInteger u = this;
1.1137 + MutableBigInteger r = new MutableBigInteger();
1.1138 +
1.1139 + // step B1
1.1140 + int s1 = u.getLowestSetBit();
1.1141 + int s2 = v.getLowestSetBit();
1.1142 + int k = (s1 < s2) ? s1 : s2;
1.1143 + if (k != 0) {
1.1144 + u.rightShift(k);
1.1145 + v.rightShift(k);
1.1146 + }
1.1147 +
1.1148 + // step B2
1.1149 + boolean uOdd = (k==s1);
1.1150 + MutableBigInteger t = uOdd ? v: u;
1.1151 + int tsign = uOdd ? -1 : 1;
1.1152 +
1.1153 + int lb;
1.1154 + while ((lb = t.getLowestSetBit()) >= 0) {
1.1155 + // steps B3 and B4
1.1156 + t.rightShift(lb);
1.1157 + // step B5
1.1158 + if (tsign > 0)
1.1159 + u = t;
1.1160 + else
1.1161 + v = t;
1.1162 +
1.1163 + // Special case one word numbers
1.1164 + if (u.intLen < 2 && v.intLen < 2) {
1.1165 + int x = u.value[u.offset];
1.1166 + int y = v.value[v.offset];
1.1167 + x = binaryGcd(x, y);
1.1168 + r.value[0] = x;
1.1169 + r.intLen = 1;
1.1170 + r.offset = 0;
1.1171 + if (k > 0)
1.1172 + r.leftShift(k);
1.1173 + return r;
1.1174 + }
1.1175 +
1.1176 + // step B6
1.1177 + if ((tsign = u.difference(v)) == 0)
1.1178 + break;
1.1179 + t = (tsign >= 0) ? u : v;
1.1180 + }
1.1181 +
1.1182 + if (k > 0)
1.1183 + u.leftShift(k);
1.1184 + return u;
1.1185 + }
1.1186 +
1.1187 + /**
1.1188 + * Calculate GCD of a and b interpreted as unsigned integers.
1.1189 + */
1.1190 + static int binaryGcd(int a, int b) {
1.1191 + if (b==0)
1.1192 + return a;
1.1193 + if (a==0)
1.1194 + return b;
1.1195 +
1.1196 + // Right shift a & b till their last bits equal to 1.
1.1197 + int aZeros = Integer.numberOfTrailingZeros(a);
1.1198 + int bZeros = Integer.numberOfTrailingZeros(b);
1.1199 + a >>>= aZeros;
1.1200 + b >>>= bZeros;
1.1201 +
1.1202 + int t = (aZeros < bZeros ? aZeros : bZeros);
1.1203 +
1.1204 + while (a != b) {
1.1205 + if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
1.1206 + a -= b;
1.1207 + a >>>= Integer.numberOfTrailingZeros(a);
1.1208 + } else {
1.1209 + b -= a;
1.1210 + b >>>= Integer.numberOfTrailingZeros(b);
1.1211 + }
1.1212 + }
1.1213 + return a<<t;
1.1214 + }
1.1215 +
1.1216 + /**
1.1217 + * Returns the modInverse of this mod p. This and p are not affected by
1.1218 + * the operation.
1.1219 + */
1.1220 + MutableBigInteger mutableModInverse(MutableBigInteger p) {
1.1221 + // Modulus is odd, use Schroeppel's algorithm
1.1222 + if (p.isOdd())
1.1223 + return modInverse(p);
1.1224 +
1.1225 + // Base and modulus are even, throw exception
1.1226 + if (isEven())
1.1227 + throw new ArithmeticException("BigInteger not invertible.");
1.1228 +
1.1229 + // Get even part of modulus expressed as a power of 2
1.1230 + int powersOf2 = p.getLowestSetBit();
1.1231 +
1.1232 + // Construct odd part of modulus
1.1233 + MutableBigInteger oddMod = new MutableBigInteger(p);
1.1234 + oddMod.rightShift(powersOf2);
1.1235 +
1.1236 + if (oddMod.isOne())
1.1237 + return modInverseMP2(powersOf2);
1.1238 +
1.1239 + // Calculate 1/a mod oddMod
1.1240 + MutableBigInteger oddPart = modInverse(oddMod);
1.1241 +
1.1242 + // Calculate 1/a mod evenMod
1.1243 + MutableBigInteger evenPart = modInverseMP2(powersOf2);
1.1244 +
1.1245 + // Combine the results using Chinese Remainder Theorem
1.1246 + MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1.1247 + MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1.1248 +
1.1249 + MutableBigInteger temp1 = new MutableBigInteger();
1.1250 + MutableBigInteger temp2 = new MutableBigInteger();
1.1251 + MutableBigInteger result = new MutableBigInteger();
1.1252 +
1.1253 + oddPart.leftShift(powersOf2);
1.1254 + oddPart.multiply(y1, result);
1.1255 +
1.1256 + evenPart.multiply(oddMod, temp1);
1.1257 + temp1.multiply(y2, temp2);
1.1258 +
1.1259 + result.add(temp2);
1.1260 + return result.divide(p, temp1);
1.1261 + }
1.1262 +
1.1263 + /*
1.1264 + * Calculate the multiplicative inverse of this mod 2^k.
1.1265 + */
1.1266 + MutableBigInteger modInverseMP2(int k) {
1.1267 + if (isEven())
1.1268 + throw new ArithmeticException("Non-invertible. (GCD != 1)");
1.1269 +
1.1270 + if (k > 64)
1.1271 + return euclidModInverse(k);
1.1272 +
1.1273 + int t = inverseMod32(value[offset+intLen-1]);
1.1274 +
1.1275 + if (k < 33) {
1.1276 + t = (k == 32 ? t : t & ((1 << k) - 1));
1.1277 + return new MutableBigInteger(t);
1.1278 + }
1.1279 +
1.1280 + long pLong = (value[offset+intLen-1] & LONG_MASK);
1.1281 + if (intLen > 1)
1.1282 + pLong |= ((long)value[offset+intLen-2] << 32);
1.1283 + long tLong = t & LONG_MASK;
1.1284 + tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
1.1285 + tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
1.1286 +
1.1287 + MutableBigInteger result = new MutableBigInteger(new int[2]);
1.1288 + result.value[0] = (int)(tLong >>> 32);
1.1289 + result.value[1] = (int)tLong;
1.1290 + result.intLen = 2;
1.1291 + result.normalize();
1.1292 + return result;
1.1293 + }
1.1294 +
1.1295 + /*
1.1296 + * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
1.1297 + */
1.1298 + static int inverseMod32(int val) {
1.1299 + // Newton's iteration!
1.1300 + int t = val;
1.1301 + t *= 2 - val*t;
1.1302 + t *= 2 - val*t;
1.1303 + t *= 2 - val*t;
1.1304 + t *= 2 - val*t;
1.1305 + return t;
1.1306 + }
1.1307 +
1.1308 + /*
1.1309 + * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
1.1310 + */
1.1311 + static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
1.1312 + // Copy the mod to protect original
1.1313 + return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
1.1314 + }
1.1315 +
1.1316 + /**
1.1317 + * Calculate the multiplicative inverse of this mod mod, where mod is odd.
1.1318 + * This and mod are not changed by the calculation.
1.1319 + *
1.1320 + * This method implements an algorithm due to Richard Schroeppel, that uses
1.1321 + * the same intermediate representation as Montgomery Reduction
1.1322 + * ("Montgomery Form"). The algorithm is described in an unpublished
1.1323 + * manuscript entitled "Fast Modular Reciprocals."
1.1324 + */
1.1325 + private MutableBigInteger modInverse(MutableBigInteger mod) {
1.1326 + MutableBigInteger p = new MutableBigInteger(mod);
1.1327 + MutableBigInteger f = new MutableBigInteger(this);
1.1328 + MutableBigInteger g = new MutableBigInteger(p);
1.1329 + SignedMutableBigInteger c = new SignedMutableBigInteger(1);
1.1330 + SignedMutableBigInteger d = new SignedMutableBigInteger();
1.1331 + MutableBigInteger temp = null;
1.1332 + SignedMutableBigInteger sTemp = null;
1.1333 +
1.1334 + int k = 0;
1.1335 + // Right shift f k times until odd, left shift d k times
1.1336 + if (f.isEven()) {
1.1337 + int trailingZeros = f.getLowestSetBit();
1.1338 + f.rightShift(trailingZeros);
1.1339 + d.leftShift(trailingZeros);
1.1340 + k = trailingZeros;
1.1341 + }
1.1342 +
1.1343 + // The Almost Inverse Algorithm
1.1344 + while(!f.isOne()) {
1.1345 + // If gcd(f, g) != 1, number is not invertible modulo mod
1.1346 + if (f.isZero())
1.1347 + throw new ArithmeticException("BigInteger not invertible.");
1.1348 +
1.1349 + // If f < g exchange f, g and c, d
1.1350 + if (f.compare(g) < 0) {
1.1351 + temp = f; f = g; g = temp;
1.1352 + sTemp = d; d = c; c = sTemp;
1.1353 + }
1.1354 +
1.1355 + // If f == g (mod 4)
1.1356 + if (((f.value[f.offset + f.intLen - 1] ^
1.1357 + g.value[g.offset + g.intLen - 1]) & 3) == 0) {
1.1358 + f.subtract(g);
1.1359 + c.signedSubtract(d);
1.1360 + } else { // If f != g (mod 4)
1.1361 + f.add(g);
1.1362 + c.signedAdd(d);
1.1363 + }
1.1364 +
1.1365 + // Right shift f k times until odd, left shift d k times
1.1366 + int trailingZeros = f.getLowestSetBit();
1.1367 + f.rightShift(trailingZeros);
1.1368 + d.leftShift(trailingZeros);
1.1369 + k += trailingZeros;
1.1370 + }
1.1371 +
1.1372 + while (c.sign < 0)
1.1373 + c.signedAdd(p);
1.1374 +
1.1375 + return fixup(c, p, k);
1.1376 + }
1.1377 +
1.1378 + /*
1.1379 + * The Fixup Algorithm
1.1380 + * Calculates X such that X = C * 2^(-k) (mod P)
1.1381 + * Assumes C<P and P is odd.
1.1382 + */
1.1383 + static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
1.1384 + int k) {
1.1385 + MutableBigInteger temp = new MutableBigInteger();
1.1386 + // Set r to the multiplicative inverse of p mod 2^32
1.1387 + int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
1.1388 +
1.1389 + for(int i=0, numWords = k >> 5; i<numWords; i++) {
1.1390 + // V = R * c (mod 2^j)
1.1391 + int v = r * c.value[c.offset + c.intLen-1];
1.1392 + // c = c + (v * p)
1.1393 + p.mul(v, temp);
1.1394 + c.add(temp);
1.1395 + // c = c / 2^j
1.1396 + c.intLen--;
1.1397 + }
1.1398 + int numBits = k & 0x1f;
1.1399 + if (numBits != 0) {
1.1400 + // V = R * c (mod 2^j)
1.1401 + int v = r * c.value[c.offset + c.intLen-1];
1.1402 + v &= ((1<<numBits) - 1);
1.1403 + // c = c + (v * p)
1.1404 + p.mul(v, temp);
1.1405 + c.add(temp);
1.1406 + // c = c / 2^j
1.1407 + c.rightShift(numBits);
1.1408 + }
1.1409 +
1.1410 + // In theory, c may be greater than p at this point (Very rare!)
1.1411 + while (c.compare(p) >= 0)
1.1412 + c.subtract(p);
1.1413 +
1.1414 + return c;
1.1415 + }
1.1416 +
1.1417 + /**
1.1418 + * Uses the extended Euclidean algorithm to compute the modInverse of base
1.1419 + * mod a modulus that is a power of 2. The modulus is 2^k.
1.1420 + */
1.1421 + MutableBigInteger euclidModInverse(int k) {
1.1422 + MutableBigInteger b = new MutableBigInteger(1);
1.1423 + b.leftShift(k);
1.1424 + MutableBigInteger mod = new MutableBigInteger(b);
1.1425 +
1.1426 + MutableBigInteger a = new MutableBigInteger(this);
1.1427 + MutableBigInteger q = new MutableBigInteger();
1.1428 + MutableBigInteger r = b.divide(a, q);
1.1429 +
1.1430 + MutableBigInteger swapper = b;
1.1431 + // swap b & r
1.1432 + b = r;
1.1433 + r = swapper;
1.1434 +
1.1435 + MutableBigInteger t1 = new MutableBigInteger(q);
1.1436 + MutableBigInteger t0 = new MutableBigInteger(1);
1.1437 + MutableBigInteger temp = new MutableBigInteger();
1.1438 +
1.1439 + while (!b.isOne()) {
1.1440 + r = a.divide(b, q);
1.1441 +
1.1442 + if (r.intLen == 0)
1.1443 + throw new ArithmeticException("BigInteger not invertible.");
1.1444 +
1.1445 + swapper = r;
1.1446 + a = swapper;
1.1447 +
1.1448 + if (q.intLen == 1)
1.1449 + t1.mul(q.value[q.offset], temp);
1.1450 + else
1.1451 + q.multiply(t1, temp);
1.1452 + swapper = q;
1.1453 + q = temp;
1.1454 + temp = swapper;
1.1455 + t0.add(q);
1.1456 +
1.1457 + if (a.isOne())
1.1458 + return t0;
1.1459 +
1.1460 + r = b.divide(a, q);
1.1461 +
1.1462 + if (r.intLen == 0)
1.1463 + throw new ArithmeticException("BigInteger not invertible.");
1.1464 +
1.1465 + swapper = b;
1.1466 + b = r;
1.1467 +
1.1468 + if (q.intLen == 1)
1.1469 + t0.mul(q.value[q.offset], temp);
1.1470 + else
1.1471 + q.multiply(t0, temp);
1.1472 + swapper = q; q = temp; temp = swapper;
1.1473 +
1.1474 + t1.add(q);
1.1475 + }
1.1476 + mod.subtract(t1);
1.1477 + return mod;
1.1478 + }
1.1479 +
1.1480 +}