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