https://github.com/google/guava/guava: Guava is a suite of core and expanded libraries that include utility classes, google's collections, io classes, and much much more. 
/*
 * Copyright (C) 2011 The Guava Authors
 *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
 * in compliance with the License. You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software distributed under the License
 * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
 * or implied. See the License for the specific language governing permissions and limitations under
 * the License.
 */

package com.google.common.math;

import static com.google.common.base.Preconditions.checkArgument;
import static com.google.common.base.Preconditions.checkNotNull;
import static com.google.common.math.MathPreconditions.checkNonNegative;
import static com.google.common.math.MathPreconditions.checkPositive;
import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
import static java.math.RoundingMode.CEILING;
import static java.math.RoundingMode.FLOOR;
import static java.math.RoundingMode.HALF_EVEN;

import com.google.common.annotations.Beta;
import com.google.common.annotations.GwtCompatible;
import com.google.common.annotations.GwtIncompatible;
import com.google.common.annotations.VisibleForTesting;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;


A class for arithmetic on values of type BigInteger. 
The implementations of many methods in this class are based on material from Henry S. Warren,
Jr.'s Hacker's Delight, (Addison Wesley, 2002).

Similar functionality for int and for long can be found in IntMath and LongMath respectively. 
Author:Louis Wasserman
Since:11.0
/**
 * A class for arithmetic on values of type {@code BigInteger}.
 *
 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
 *
 * <p>Similar functionality for {@code int} and for {@code long} can be found in {@link IntMath} and
 * {@link LongMath} respectively.
 *
 * @author Louis Wasserman
 * @since 11.0
 */
@GwtCompatible(emulated = true)
public final class BigIntegerMath {
  
Returns the smallest power of two greater than or equal to x. This is equivalent to BigInteger.valueOf(2).pow(log2(x, CEILING)). 
Throws:
Since:20.0
/**
   * Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to
   * {@code BigInteger.valueOf(2).pow(log2(x, CEILING))}.
   *
   * @throws IllegalArgumentException if {@code x <= 0}
   * @since 20.0
   */
  @Beta
  public static BigInteger ceilingPowerOfTwo(BigInteger x) {
    return BigInteger.ZERO.setBit(log2(x, RoundingMode.CEILING));
  }

  
Returns the largest power of two less than or equal to x. This is equivalent to 
BigInteger.valueOf(2).pow(log2(x, FLOOR)). 
Throws:
Since:20.0
/**
   * Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code
   * BigInteger.valueOf(2).pow(log2(x, FLOOR))}.
   *
   * @throws IllegalArgumentException if {@code x <= 0}
   * @since 20.0
   */
  @Beta
  public static BigInteger floorPowerOfTwo(BigInteger x) {
    return BigInteger.ZERO.setBit(log2(x, RoundingMode.FLOOR));
  }

  
Returns true if x represents a power of two. 
/** Returns {@code true} if {@code x} represents a power of two. */
  public static boolean isPowerOfTwo(BigInteger x) {
    checkNotNull(x);
    return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1;
  }

  
Returns the base-2 logarithm of x, rounded according to the specified rounding mode. 
Throws:
/**
   * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
   *
   * @throws IllegalArgumentException if {@code x <= 0}
   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
   *     is not a power of two
   */
  @SuppressWarnings("fallthrough")
  // TODO(kevinb): remove after this warning is disabled globally
  public static int log2(BigInteger x, RoundingMode mode) {
    checkPositive("x", checkNotNull(x));
    int logFloor = x.bitLength() - 1;
    switch (mode) {
      case UNNECESSARY:
        checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through
      case DOWN:
      case FLOOR:
        return logFloor;

      case UP:
      case CEILING:
        return isPowerOfTwo(x) ? logFloor : logFloor + 1;

      case HALF_DOWN:
      case HALF_UP:
      case HALF_EVEN:
        if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) {
          BigInteger halfPower =
              SQRT2_PRECOMPUTED_BITS.shiftRight(SQRT2_PRECOMPUTE_THRESHOLD - logFloor);
          if (x.compareTo(halfPower) <= 0) {
            return logFloor;
          } else {
            return logFloor + 1;
          }
        }
        // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
        //
        // To determine which side of logFloor.5 the logarithm is,
        // we compare x^2 to 2^(2 * logFloor + 1).
        BigInteger x2 = x.pow(2);
        int logX2Floor = x2.bitLength() - 1;
        return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1;

      default:
        throw new AssertionError();
    }
  }

  /*
   * The maximum number of bits in a square root for which we'll precompute an explicit half power
   * of two. This can be any value, but higher values incur more class load time and linearly
   * increasing memory consumption.
   */
  @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256;

  @VisibleForTesting
  static final BigInteger SQRT2_PRECOMPUTED_BITS =
      new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16);

  
Returns the base-10 logarithm of x, rounded according to the specified rounding mode. 
Throws:
/**
   * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
   *
   * @throws IllegalArgumentException if {@code x <= 0}
   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
   *     is not a power of ten
   */
  @GwtIncompatible // TODO
  @SuppressWarnings("fallthrough")
  public static int log10(BigInteger x, RoundingMode mode) {
    checkPositive("x", x);
    if (fitsInLong(x)) {
      return LongMath.log10(x.longValue(), mode);
    }

    int approxLog10 = (int) (log2(x, FLOOR) * LN_2 / LN_10);
    BigInteger approxPow = BigInteger.TEN.pow(approxLog10);
    int approxCmp = approxPow.compareTo(x);

    /*
     * We adjust approxLog10 and approxPow until they're equal to floor(log10(x)) and
     * 10^floor(log10(x)).
     */

    if (approxCmp > 0) {
      /*
       * The code is written so that even completely incorrect approximations will still yield the
       * correct answer eventually, but in practice this branch should almost never be entered, and
       * even then the loop should not run more than once.
       */
      do {
        approxLog10--;
        approxPow = approxPow.divide(BigInteger.TEN);
        approxCmp = approxPow.compareTo(x);
      } while (approxCmp > 0);
    } else {
      BigInteger nextPow = BigInteger.TEN.multiply(approxPow);
      int nextCmp = nextPow.compareTo(x);
      while (nextCmp <= 0) {
        approxLog10++;
        approxPow = nextPow;
        approxCmp = nextCmp;
        nextPow = BigInteger.TEN.multiply(approxPow);
        nextCmp = nextPow.compareTo(x);
      }
    }

    int floorLog = approxLog10;
    BigInteger floorPow = approxPow;
    int floorCmp = approxCmp;

    switch (mode) {
      case UNNECESSARY:
        checkRoundingUnnecessary(floorCmp == 0);
        // fall through
      case FLOOR:
      case DOWN:
        return floorLog;

      case CEILING:
      case UP:
        return floorPow.equals(x) ? floorLog : floorLog + 1;

      case HALF_DOWN:
      case HALF_UP:
      case HALF_EVEN:
        // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5
        BigInteger x2 = x.pow(2);
        BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN);
        return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1;
      default:
        throw new AssertionError();
    }
  }

  private static final double LN_10 = Math.log(10);
  private static final double LN_2 = Math.log(2);

  
Returns the square root of x, rounded with the specified rounding mode. 
Throws:
/**
   * Returns the square root of {@code x}, rounded with the specified rounding mode.
   *
   * @throws IllegalArgumentException if {@code x < 0}
   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
   *     sqrt(x)} is not an integer
   */
  @GwtIncompatible // TODO
  @SuppressWarnings("fallthrough")
  public static BigInteger sqrt(BigInteger x, RoundingMode mode) {
    checkNonNegative("x", x);
    if (fitsInLong(x)) {
      return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode));
    }
    BigInteger sqrtFloor = sqrtFloor(x);
    switch (mode) {
      case UNNECESSARY:
        checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through
      case FLOOR:
      case DOWN:
        return sqrtFloor;
      case CEILING:
      case UP:
        int sqrtFloorInt = sqrtFloor.intValue();
        boolean sqrtFloorIsExact =
            (sqrtFloorInt * sqrtFloorInt == x.intValue()) // fast check mod 2^32
                && sqrtFloor.pow(2).equals(x); // slow exact check
        return sqrtFloorIsExact ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
      case HALF_DOWN:
      case HALF_UP:
      case HALF_EVEN:
        BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor);
        /*
         * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
         * and halfSquare are integers, this is equivalent to testing whether or not x <=
         * halfSquare.
         */
        return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
      default:
        throw new AssertionError();
    }
  }

  @GwtIncompatible // TODO
  private static BigInteger sqrtFloor(BigInteger x) {
    /*
     * Adapted from Hacker's Delight, Figure 11-1.
     *
     * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and
     * then we can get a double approximation of the square root. Then, we iteratively improve this
     * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2.
     * This iteration has the following two properties:
     *
     * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is
     * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean,
     * and the arithmetic mean is always higher than the geometric mean.
     *
     * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles
     * with each iteration, so this algorithm takes O(log(digits)) iterations.
     *
     * We start out with a double-precision approximation, which may be higher or lower than the
     * true value. Therefore, we perform at least one Newton iteration to get a guess that's
     * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point.
     */
    BigInteger sqrt0;
    int log2 = log2(x, FLOOR);
    if (log2 < Double.MAX_EXPONENT) {
      sqrt0 = sqrtApproxWithDoubles(x);
    } else {
      int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even!
      /*
       * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be
       * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift).
       */
      sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1);
    }
    BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
    if (sqrt0.equals(sqrt1)) {
      return sqrt0;
    }
    do {
      sqrt0 = sqrt1;
      sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
    } while (sqrt1.compareTo(sqrt0) < 0);
    return sqrt0;
  }

  @GwtIncompatible // TODO
  private static BigInteger sqrtApproxWithDoubles(BigInteger x) {
    return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN);
  }

  
Returns the result of dividing p by q, rounding using the specified 
RoundingMode. 
Throws:
  • ArithmeticException – if q == 0, or if mode == UNNECESSARY and a is not an integer multiple of b
/**
   * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
   * RoundingMode}.
   *
   * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
   *     is not an integer multiple of {@code b}
   */
  @GwtIncompatible // TODO
  public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode) {
    BigDecimal pDec = new BigDecimal(p);
    BigDecimal qDec = new BigDecimal(q);
    return pDec.divide(qDec, 0, mode).toBigIntegerExact();
  }

  
Returns n!, that is, the product of the first n positive integers, or 1 if n == 0. 
Warning: the result takes O(n log n) space, so use cautiously.

This uses an efficient binary recursive algorithm to compute the factorial with balanced
multiplies. It also removes all the 2s from the intermediate products (shifting them back in at
the end).

Throws:
/**
   * Returns {@code n!}, that is, the product of the first {@code n} positive integers, or {@code 1}
   * if {@code n == 0}.
   *
   * <p><b>Warning:</b> the result takes <i>O(n log n)</i> space, so use cautiously.
   *
   * <p>This uses an efficient binary recursive algorithm to compute the factorial with balanced
   * multiplies. It also removes all the 2s from the intermediate products (shifting them back in at
   * the end).
   *
   * @throws IllegalArgumentException if {@code n < 0}
   */
  public static BigInteger factorial(int n) {
    checkNonNegative("n", n);

    // If the factorial is small enough, just use LongMath to do it.
    if (n < LongMath.factorials.length) {
      return BigInteger.valueOf(LongMath.factorials[n]);
    }

    // Pre-allocate space for our list of intermediate BigIntegers.
    int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING);
    ArrayList<BigInteger> bignums = new ArrayList<>(approxSize);

    // Start from the pre-computed maximum long factorial.
    int startingNumber = LongMath.factorials.length;
    long product = LongMath.factorials[startingNumber - 1];
    // Strip off 2s from this value.
    int shift = Long.numberOfTrailingZeros(product);
    product >>= shift;

    // Use floor(log2(num)) + 1 to prevent overflow of multiplication.
    int productBits = LongMath.log2(product, FLOOR) + 1;
    int bits = LongMath.log2(startingNumber, FLOOR) + 1;
    // Check for the next power of two boundary, to save us a CLZ operation.
    int nextPowerOfTwo = 1 << (bits - 1);

    // Iteratively multiply the longs as big as they can go.
    for (long num = startingNumber; num <= n; num++) {
      // Check to see if the floor(log2(num)) + 1 has changed.
      if ((num & nextPowerOfTwo) != 0) {
        nextPowerOfTwo <<= 1;
        bits++;
      }
      // Get rid of the 2s in num.
      int tz = Long.numberOfTrailingZeros(num);
      long normalizedNum = num >> tz;
      shift += tz;
      // Adjust floor(log2(num)) + 1.
      int normalizedBits = bits - tz;
      // If it won't fit in a long, then we store off the intermediate product.
      if (normalizedBits + productBits >= Long.SIZE) {
        bignums.add(BigInteger.valueOf(product));
        product = 1;
        productBits = 0;
      }
      product *= normalizedNum;
      productBits = LongMath.log2(product, FLOOR) + 1;
    }
    // Check for leftovers.
    if (product > 1) {
      bignums.add(BigInteger.valueOf(product));
    }
    // Efficiently multiply all the intermediate products together.
    return listProduct(bignums).shiftLeft(shift);
  }

  static BigInteger listProduct(List<BigInteger> nums) {
    return listProduct(nums, 0, nums.size());
  }

  static BigInteger listProduct(List<BigInteger> nums, int start, int end) {
    switch (end - start) {
      case 0:
        return BigInteger.ONE;
      case 1:
        return nums.get(start);
      case 2:
        return nums.get(start).multiply(nums.get(start + 1));
      case 3:
        return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2));
      default:
        // Otherwise, split the list in half and recursively do this.
        int m = (end + start) >>> 1;
        return listProduct(nums, start, m).multiply(listProduct(nums, m, end));
    }
  }

  
Returns n choose k, also known as the binomial coefficient of n and k, that is, n! / (k! (n - k)!). 
Warning: the result can take as much as O(k log n) space.

Throws:
/**
   * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
   * {@code k}, that is, {@code n! / (k! (n - k)!)}.
   *
   * <p><b>Warning:</b> the result can take as much as <i>O(k log n)</i> space.
   *
   * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
   */
  public static BigInteger binomial(int n, int k) {
    checkNonNegative("n", n);
    checkNonNegative("k", k);
    checkArgument(k <= n, "k (%s) > n (%s)", k, n);
    if (k > (n >> 1)) {
      k = n - k;
    }
    if (k < LongMath.biggestBinomials.length && n <= LongMath.biggestBinomials[k]) {
      return BigInteger.valueOf(LongMath.binomial(n, k));
    }

    BigInteger accum = BigInteger.ONE;

    long numeratorAccum = n;
    long denominatorAccum = 1;

    int bits = LongMath.log2(n, RoundingMode.CEILING);

    int numeratorBits = bits;

    for (int i = 1; i < k; i++) {
      int p = n - i;
      int q = i + 1;

      // log2(p) >= bits - 1, because p >= n/2

      if (numeratorBits + bits >= Long.SIZE - 1) {
        // The numerator is as big as it can get without risking overflow.
        // Multiply numeratorAccum / denominatorAccum into accum.
        accum =
            accum
                .multiply(BigInteger.valueOf(numeratorAccum))
                .divide(BigInteger.valueOf(denominatorAccum));
        numeratorAccum = p;
        denominatorAccum = q;
        numeratorBits = bits;
      } else {
        // We can definitely multiply into the long accumulators without overflowing them.
        numeratorAccum *= p;
        denominatorAccum *= q;
        numeratorBits += bits;
      }
    }
    return accum
        .multiply(BigInteger.valueOf(numeratorAccum))
        .divide(BigInteger.valueOf(denominatorAccum));
  }

  // Returns true if BigInteger.valueOf(x.longValue()).equals(x).
  @GwtIncompatible // TODO
  static boolean fitsInLong(BigInteger x) {
    return x.bitLength() <= Long.SIZE - 1;
  }

  private BigIntegerMath() {}
}