001/*
002 * Copyright (C) 2011 The Guava Authors
003 *
004 * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
005 * in compliance with the License. You may obtain a copy of the License at
006 *
007 * http://www.apache.org/licenses/LICENSE-2.0
008 *
009 * Unless required by applicable law or agreed to in writing, software distributed under the License
010 * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
011 * or implied. See the License for the specific language governing permissions and limitations under
012 * the License.
013 */
014
015package com.google.common.math;
016
017import static com.google.common.base.Preconditions.checkArgument;
018import static com.google.common.base.Preconditions.checkNotNull;
019import static com.google.common.math.MathPreconditions.checkNonNegative;
020import static com.google.common.math.MathPreconditions.checkPositive;
021import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
022import static java.math.RoundingMode.CEILING;
023import static java.math.RoundingMode.FLOOR;
024import static java.math.RoundingMode.HALF_EVEN;
025
026import com.google.common.annotations.GwtCompatible;
027import com.google.common.annotations.GwtIncompatible;
028import com.google.common.annotations.VisibleForTesting;
029import java.math.BigDecimal;
030import java.math.BigInteger;
031import java.math.RoundingMode;
032import java.util.ArrayList;
033import java.util.List;
034
035/**
036 * A class for arithmetic on values of type {@code BigInteger}.
037 *
038 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
039 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
040 *
041 * <p>Similar functionality for {@code int} and for {@code long} can be found in {@link IntMath} and
042 * {@link LongMath} respectively.
043 *
044 * @author Louis Wasserman
045 * @since 11.0
046 */
047@GwtCompatible(emulated = true)
048@ElementTypesAreNonnullByDefault
049public final class BigIntegerMath {
050  /**
051   * Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to
052   * {@code BigInteger.valueOf(2).pow(log2(x, CEILING))}.
053   *
054   * @throws IllegalArgumentException if {@code x <= 0}
055   * @since 20.0
056   */
057  public static BigInteger ceilingPowerOfTwo(BigInteger x) {
058    return BigInteger.ZERO.setBit(log2(x, CEILING));
059  }
060
061  /**
062   * Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code
063   * BigInteger.valueOf(2).pow(log2(x, FLOOR))}.
064   *
065   * @throws IllegalArgumentException if {@code x <= 0}
066   * @since 20.0
067   */
068  public static BigInteger floorPowerOfTwo(BigInteger x) {
069    return BigInteger.ZERO.setBit(log2(x, FLOOR));
070  }
071
072  /** Returns {@code true} if {@code x} represents a power of two. */
073  public static boolean isPowerOfTwo(BigInteger x) {
074    checkNotNull(x);
075    return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1;
076  }
077
078  /**
079   * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
080   *
081   * @throws IllegalArgumentException if {@code x <= 0}
082   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
083   *     is not a power of two
084   */
085  @SuppressWarnings("fallthrough")
086  // TODO(kevinb): remove after this warning is disabled globally
087  public static int log2(BigInteger x, RoundingMode mode) {
088    checkPositive("x", checkNotNull(x));
089    int logFloor = x.bitLength() - 1;
090    switch (mode) {
091      case UNNECESSARY:
092        checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through
093      case DOWN:
094      case FLOOR:
095        return logFloor;
096
097      case UP:
098      case CEILING:
099        return isPowerOfTwo(x) ? logFloor : logFloor + 1;
100
101      case HALF_DOWN:
102      case HALF_UP:
103      case HALF_EVEN:
104        if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) {
105          BigInteger halfPower =
106              SQRT2_PRECOMPUTED_BITS.shiftRight(SQRT2_PRECOMPUTE_THRESHOLD - logFloor);
107          if (x.compareTo(halfPower) <= 0) {
108            return logFloor;
109          } else {
110            return logFloor + 1;
111          }
112        }
113        // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
114        //
115        // To determine which side of logFloor.5 the logarithm is,
116        // we compare x^2 to 2^(2 * logFloor + 1).
117        BigInteger x2 = x.pow(2);
118        int logX2Floor = x2.bitLength() - 1;
119        return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1;
120
121      default:
122        throw new AssertionError();
123    }
124  }
125
126  /*
127   * The maximum number of bits in a square root for which we'll precompute an explicit half power
128   * of two. This can be any value, but higher values incur more class load time and linearly
129   * increasing memory consumption.
130   */
131  @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256;
132
133  @VisibleForTesting
134  static final BigInteger SQRT2_PRECOMPUTED_BITS =
135      new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16);
136
137  /**
138   * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
139   *
140   * @throws IllegalArgumentException if {@code x <= 0}
141   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
142   *     is not a power of ten
143   */
144  @GwtIncompatible // TODO
145  @SuppressWarnings("fallthrough")
146  public static int log10(BigInteger x, RoundingMode mode) {
147    checkPositive("x", x);
148    if (fitsInLong(x)) {
149      return LongMath.log10(x.longValue(), mode);
150    }
151
152    int approxLog10 = (int) (log2(x, FLOOR) * LN_2 / LN_10);
153    BigInteger approxPow = BigInteger.TEN.pow(approxLog10);
154    int approxCmp = approxPow.compareTo(x);
155
156    /*
157     * We adjust approxLog10 and approxPow until they're equal to floor(log10(x)) and
158     * 10^floor(log10(x)).
159     */
160
161    if (approxCmp > 0) {
162      /*
163       * The code is written so that even completely incorrect approximations will still yield the
164       * correct answer eventually, but in practice this branch should almost never be entered, and
165       * even then the loop should not run more than once.
166       */
167      do {
168        approxLog10--;
169        approxPow = approxPow.divide(BigInteger.TEN);
170        approxCmp = approxPow.compareTo(x);
171      } while (approxCmp > 0);
172    } else {
173      BigInteger nextPow = BigInteger.TEN.multiply(approxPow);
174      int nextCmp = nextPow.compareTo(x);
175      while (nextCmp <= 0) {
176        approxLog10++;
177        approxPow = nextPow;
178        approxCmp = nextCmp;
179        nextPow = BigInteger.TEN.multiply(approxPow);
180        nextCmp = nextPow.compareTo(x);
181      }
182    }
183
184    int floorLog = approxLog10;
185    BigInteger floorPow = approxPow;
186    int floorCmp = approxCmp;
187
188    switch (mode) {
189      case UNNECESSARY:
190        checkRoundingUnnecessary(floorCmp == 0);
191        // fall through
192      case FLOOR:
193      case DOWN:
194        return floorLog;
195
196      case CEILING:
197      case UP:
198        return floorPow.equals(x) ? floorLog : floorLog + 1;
199
200      case HALF_DOWN:
201      case HALF_UP:
202      case HALF_EVEN:
203        // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5
204        BigInteger x2 = x.pow(2);
205        BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN);
206        return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1;
207      default:
208        throw new AssertionError();
209    }
210  }
211
212  private static final double LN_10 = Math.log(10);
213  private static final double LN_2 = Math.log(2);
214
215  /**
216   * Returns the square root of {@code x}, rounded with the specified rounding mode.
217   *
218   * @throws IllegalArgumentException if {@code x < 0}
219   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
220   *     sqrt(x)} is not an integer
221   */
222  @GwtIncompatible // TODO
223  @SuppressWarnings("fallthrough")
224  public static BigInteger sqrt(BigInteger x, RoundingMode mode) {
225    checkNonNegative("x", x);
226    if (fitsInLong(x)) {
227      return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode));
228    }
229    BigInteger sqrtFloor = sqrtFloor(x);
230    switch (mode) {
231      case UNNECESSARY:
232        checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through
233      case FLOOR:
234      case DOWN:
235        return sqrtFloor;
236      case CEILING:
237      case UP:
238        int sqrtFloorInt = sqrtFloor.intValue();
239        boolean sqrtFloorIsExact =
240            (sqrtFloorInt * sqrtFloorInt == x.intValue()) // fast check mod 2^32
241                && sqrtFloor.pow(2).equals(x); // slow exact check
242        return sqrtFloorIsExact ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
243      case HALF_DOWN:
244      case HALF_UP:
245      case HALF_EVEN:
246        BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor);
247        /*
248         * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
249         * and halfSquare are integers, this is equivalent to testing whether or not x <=
250         * halfSquare.
251         */
252        return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
253      default:
254        throw new AssertionError();
255    }
256  }
257
258  @GwtIncompatible // TODO
259  private static BigInteger sqrtFloor(BigInteger x) {
260    /*
261     * Adapted from Hacker's Delight, Figure 11-1.
262     *
263     * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and
264     * then we can get a double approximation of the square root. Then, we iteratively improve this
265     * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2.
266     * This iteration has the following two properties:
267     *
268     * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is
269     * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean,
270     * and the arithmetic mean is always higher than the geometric mean.
271     *
272     * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles
273     * with each iteration, so this algorithm takes O(log(digits)) iterations.
274     *
275     * We start out with a double-precision approximation, which may be higher or lower than the
276     * true value. Therefore, we perform at least one Newton iteration to get a guess that's
277     * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point.
278     */
279    BigInteger sqrt0;
280    int log2 = log2(x, FLOOR);
281    if (log2 < Double.MAX_EXPONENT) {
282      sqrt0 = sqrtApproxWithDoubles(x);
283    } else {
284      int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even!
285      /*
286       * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be
287       * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift).
288       */
289      sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1);
290    }
291    BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
292    if (sqrt0.equals(sqrt1)) {
293      return sqrt0;
294    }
295    do {
296      sqrt0 = sqrt1;
297      sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
298    } while (sqrt1.compareTo(sqrt0) < 0);
299    return sqrt0;
300  }
301
302  @GwtIncompatible // TODO
303  private static BigInteger sqrtApproxWithDoubles(BigInteger x) {
304    return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN);
305  }
306
307  /**
308   * Returns {@code x}, rounded to a {@code double} with the specified rounding mode. If {@code x}
309   * is precisely representable as a {@code double}, its {@code double} value will be returned;
310   * otherwise, the rounding will choose between the two nearest representable values with {@code
311   * mode}.
312   *
313   * <p>For the case of {@link RoundingMode#HALF_DOWN}, {@code HALF_UP}, and {@code HALF_EVEN},
314   * infinite {@code double} values are considered infinitely far away. For example, 2^2000 is not
315   * representable as a double, but {@code roundToDouble(BigInteger.valueOf(2).pow(2000), HALF_UP)}
316   * will return {@code Double.MAX_VALUE}, not {@code Double.POSITIVE_INFINITY}.
317   *
318   * <p>For the case of {@link RoundingMode#HALF_EVEN}, this implementation uses the IEEE 754
319   * default rounding mode: if the two nearest representable values are equally near, the one with
320   * the least significant bit zero is chosen. (In such cases, both of the nearest representable
321   * values are even integers; this method returns the one that is a multiple of a greater power of
322   * two.)
323   *
324   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
325   *     is not precisely representable as a {@code double}
326   * @since 30.0
327   */
328  @GwtIncompatible
329  public static double roundToDouble(BigInteger x, RoundingMode mode) {
330    return BigIntegerToDoubleRounder.INSTANCE.roundToDouble(x, mode);
331  }
332
333  @GwtIncompatible
334  private static class BigIntegerToDoubleRounder extends ToDoubleRounder<BigInteger> {
335    static final BigIntegerToDoubleRounder INSTANCE = new BigIntegerToDoubleRounder();
336
337    private BigIntegerToDoubleRounder() {}
338
339    @Override
340    double roundToDoubleArbitrarily(BigInteger bigInteger) {
341      return DoubleUtils.bigToDouble(bigInteger);
342    }
343
344    @Override
345    int sign(BigInteger bigInteger) {
346      return bigInteger.signum();
347    }
348
349    @Override
350    BigInteger toX(double d, RoundingMode mode) {
351      return DoubleMath.roundToBigInteger(d, mode);
352    }
353
354    @Override
355    BigInteger minus(BigInteger a, BigInteger b) {
356      return a.subtract(b);
357    }
358  }
359
360  /**
361   * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
362   * RoundingMode}.
363   *
364   * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
365   *     is not an integer multiple of {@code b}
366   */
367  @GwtIncompatible // TODO
368  public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode) {
369    BigDecimal pDec = new BigDecimal(p);
370    BigDecimal qDec = new BigDecimal(q);
371    return pDec.divide(qDec, 0, mode).toBigIntegerExact();
372  }
373
374  /**
375   * Returns {@code n!}, that is, the product of the first {@code n} positive integers, or {@code 1}
376   * if {@code n == 0}.
377   *
378   * <p><b>Warning:</b> the result takes <i>O(n log n)</i> space, so use cautiously.
379   *
380   * <p>This uses an efficient binary recursive algorithm to compute the factorial with balanced
381   * multiplies. It also removes all the 2s from the intermediate products (shifting them back in at
382   * the end).
383   *
384   * @throws IllegalArgumentException if {@code n < 0}
385   */
386  public static BigInteger factorial(int n) {
387    checkNonNegative("n", n);
388
389    // If the factorial is small enough, just use LongMath to do it.
390    if (n < LongMath.factorials.length) {
391      return BigInteger.valueOf(LongMath.factorials[n]);
392    }
393
394    // Pre-allocate space for our list of intermediate BigIntegers.
395    int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING);
396    ArrayList<BigInteger> bignums = new ArrayList<>(approxSize);
397
398    // Start from the pre-computed maximum long factorial.
399    int startingNumber = LongMath.factorials.length;
400    long product = LongMath.factorials[startingNumber - 1];
401    // Strip off 2s from this value.
402    int shift = Long.numberOfTrailingZeros(product);
403    product >>= shift;
404
405    // Use floor(log2(num)) + 1 to prevent overflow of multiplication.
406    int productBits = LongMath.log2(product, FLOOR) + 1;
407    int bits = LongMath.log2(startingNumber, FLOOR) + 1;
408    // Check for the next power of two boundary, to save us a CLZ operation.
409    int nextPowerOfTwo = 1 << (bits - 1);
410
411    // Iteratively multiply the longs as big as they can go.
412    for (long num = startingNumber; num <= n; num++) {
413      // Check to see if the floor(log2(num)) + 1 has changed.
414      if ((num & nextPowerOfTwo) != 0) {
415        nextPowerOfTwo <<= 1;
416        bits++;
417      }
418      // Get rid of the 2s in num.
419      int tz = Long.numberOfTrailingZeros(num);
420      long normalizedNum = num >> tz;
421      shift += tz;
422      // Adjust floor(log2(num)) + 1.
423      int normalizedBits = bits - tz;
424      // If it won't fit in a long, then we store off the intermediate product.
425      if (normalizedBits + productBits >= Long.SIZE) {
426        bignums.add(BigInteger.valueOf(product));
427        product = 1;
428        productBits = 0;
429      }
430      product *= normalizedNum;
431      productBits = LongMath.log2(product, FLOOR) + 1;
432    }
433    // Check for leftovers.
434    if (product > 1) {
435      bignums.add(BigInteger.valueOf(product));
436    }
437    // Efficiently multiply all the intermediate products together.
438    return listProduct(bignums).shiftLeft(shift);
439  }
440
441  static BigInteger listProduct(List<BigInteger> nums) {
442    return listProduct(nums, 0, nums.size());
443  }
444
445  static BigInteger listProduct(List<BigInteger> nums, int start, int end) {
446    switch (end - start) {
447      case 0:
448        return BigInteger.ONE;
449      case 1:
450        return nums.get(start);
451      case 2:
452        return nums.get(start).multiply(nums.get(start + 1));
453      case 3:
454        return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2));
455      default:
456        // Otherwise, split the list in half and recursively do this.
457        int m = (end + start) >>> 1;
458        return listProduct(nums, start, m).multiply(listProduct(nums, m, end));
459    }
460  }
461
462  /**
463   * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
464   * {@code k}, that is, {@code n! / (k! (n - k)!)}.
465   *
466   * <p><b>Warning:</b> the result can take as much as <i>O(k log n)</i> space.
467   *
468   * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
469   */
470  public static BigInteger binomial(int n, int k) {
471    checkNonNegative("n", n);
472    checkNonNegative("k", k);
473    checkArgument(k <= n, "k (%s) > n (%s)", k, n);
474    if (k > (n >> 1)) {
475      k = n - k;
476    }
477    if (k < LongMath.biggestBinomials.length && n <= LongMath.biggestBinomials[k]) {
478      return BigInteger.valueOf(LongMath.binomial(n, k));
479    }
480
481    BigInteger accum = BigInteger.ONE;
482
483    long numeratorAccum = n;
484    long denominatorAccum = 1;
485
486    int bits = LongMath.log2(n, CEILING);
487
488    int numeratorBits = bits;
489
490    for (int i = 1; i < k; i++) {
491      int p = n - i;
492      int q = i + 1;
493
494      // log2(p) >= bits - 1, because p >= n/2
495
496      if (numeratorBits + bits >= Long.SIZE - 1) {
497        // The numerator is as big as it can get without risking overflow.
498        // Multiply numeratorAccum / denominatorAccum into accum.
499        accum =
500            accum
501                .multiply(BigInteger.valueOf(numeratorAccum))
502                .divide(BigInteger.valueOf(denominatorAccum));
503        numeratorAccum = p;
504        denominatorAccum = q;
505        numeratorBits = bits;
506      } else {
507        // We can definitely multiply into the long accumulators without overflowing them.
508        numeratorAccum *= p;
509        denominatorAccum *= q;
510        numeratorBits += bits;
511      }
512    }
513    return accum
514        .multiply(BigInteger.valueOf(numeratorAccum))
515        .divide(BigInteger.valueOf(denominatorAccum));
516  }
517
518  // Returns true if BigInteger.valueOf(x.longValue()).equals(x).
519  @GwtIncompatible // TODO
520  static boolean fitsInLong(BigInteger x) {
521    return x.bitLength() <= Long.SIZE - 1;
522  }
523
524  private BigIntegerMath() {}
525}