• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (C) 2011 The Guava Authors
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 package com.google.common.math;
18 
19 import static com.google.common.base.Preconditions.checkArgument;
20 import static com.google.common.base.Preconditions.checkNotNull;
21 import static com.google.common.math.MathPreconditions.checkNonNegative;
22 import static com.google.common.math.MathPreconditions.checkPositive;
23 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
24 import static java.math.RoundingMode.CEILING;
25 import static java.math.RoundingMode.FLOOR;
26 import static java.math.RoundingMode.HALF_EVEN;
27 
28 import com.google.common.annotations.Beta;
29 import com.google.common.annotations.VisibleForTesting;
30 
31 import java.math.BigDecimal;
32 import java.math.BigInteger;
33 import java.math.RoundingMode;
34 import java.util.ArrayList;
35 import java.util.List;
36 
37 /**
38  * A class for arithmetic on values of type {@code BigInteger}.
39  *
40  * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
41  * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
42  *
43  * <p>Similar functionality for {@code int} and for {@code long} can be found in
44  * {@link IntMath} and {@link LongMath} respectively.
45  *
46  * @author Louis Wasserman
47  * @since 11.0
48  */
49 @Beta
50 public final class BigIntegerMath {
51   /**
52    * Returns {@code true} if {@code x} represents a power of two.
53    */
isPowerOfTwo(BigInteger x)54   public static boolean isPowerOfTwo(BigInteger x) {
55     checkNotNull(x);
56     return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1;
57   }
58 
59   /**
60    * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
61    *
62    * @throws IllegalArgumentException if {@code x <= 0}
63    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
64    *         is not a power of two
65    */
66   @SuppressWarnings("fallthrough")
log2(BigInteger x, RoundingMode mode)67   public static int log2(BigInteger x, RoundingMode mode) {
68     checkPositive("x", checkNotNull(x));
69     int logFloor = x.bitLength() - 1;
70     switch (mode) {
71       case UNNECESSARY:
72         checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through
73       case DOWN:
74       case FLOOR:
75         return logFloor;
76 
77       case UP:
78       case CEILING:
79         return isPowerOfTwo(x) ? logFloor : logFloor + 1;
80 
81       case HALF_DOWN:
82       case HALF_UP:
83       case HALF_EVEN:
84         if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) {
85           BigInteger halfPower = SQRT2_PRECOMPUTED_BITS.shiftRight(
86               SQRT2_PRECOMPUTE_THRESHOLD - logFloor);
87           if (x.compareTo(halfPower) <= 0) {
88             return logFloor;
89           } else {
90             return logFloor + 1;
91           }
92         }
93         /*
94          * Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
95          *
96          * To determine which side of logFloor.5 the logarithm is, we compare x^2 to 2^(2 *
97          * logFloor + 1).
98          */
99         BigInteger x2 = x.pow(2);
100         int logX2Floor = x2.bitLength() - 1;
101         return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1;
102 
103       default:
104         throw new AssertionError();
105     }
106   }
107 
108   /*
109    * The maximum number of bits in a square root for which we'll precompute an explicit half power
110    * of two. This can be any value, but higher values incur more class load time and linearly
111    * increasing memory consumption.
112    */
113   @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256;
114 
115   @VisibleForTesting static final BigInteger SQRT2_PRECOMPUTED_BITS =
116       new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16);
117 
118   /**
119    * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
120    *
121    * @throws IllegalArgumentException if {@code x <= 0}
122    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
123    *         is not a power of ten
124    */
125   @SuppressWarnings("fallthrough")
log10(BigInteger x, RoundingMode mode)126   public static int log10(BigInteger x, RoundingMode mode) {
127     checkPositive("x", x);
128     if (fitsInLong(x)) {
129       return LongMath.log10(x.longValue(), mode);
130     }
131 
132     // capacity of 10 suffices for all x <= 10^(2^10).
133     List<BigInteger> powersOf10 = new ArrayList<BigInteger>(10);
134     BigInteger powerOf10 = BigInteger.TEN;
135     while (x.compareTo(powerOf10) >= 0) {
136       powersOf10.add(powerOf10);
137       powerOf10 = powerOf10.pow(2);
138     }
139     BigInteger floorPow = BigInteger.ONE;
140     int floorLog = 0;
141     for (int i = powersOf10.size() - 1; i >= 0; i--) {
142       BigInteger powOf10 = powersOf10.get(i);
143       floorLog *= 2;
144       BigInteger tenPow = powOf10.multiply(floorPow);
145       if (x.compareTo(tenPow) >= 0) {
146         floorPow = tenPow;
147         floorLog++;
148       }
149     }
150     switch (mode) {
151       case UNNECESSARY:
152         checkRoundingUnnecessary(floorPow.equals(x));
153         // fall through
154       case FLOOR:
155       case DOWN:
156         return floorLog;
157 
158       case CEILING:
159       case UP:
160         return floorPow.equals(x) ? floorLog : floorLog + 1;
161 
162       case HALF_DOWN:
163       case HALF_UP:
164       case HALF_EVEN:
165         // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5
166         BigInteger x2 = x.pow(2);
167         BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN);
168         return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1;
169       default:
170         throw new AssertionError();
171     }
172   }
173 
174   /**
175    * Returns the square root of {@code x}, rounded with the specified rounding mode.
176    *
177    * @throws IllegalArgumentException if {@code x < 0}
178    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and
179    *         {@code sqrt(x)} is not an integer
180    */
181   @SuppressWarnings("fallthrough")
sqrt(BigInteger x, RoundingMode mode)182   public static BigInteger sqrt(BigInteger x, RoundingMode mode) {
183     checkNonNegative("x", x);
184     if (fitsInLong(x)) {
185       return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode));
186     }
187     BigInteger sqrtFloor = sqrtFloor(x);
188     switch (mode) {
189       case UNNECESSARY:
190         checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through
191       case FLOOR:
192       case DOWN:
193         return sqrtFloor;
194       case CEILING:
195       case UP:
196         return sqrtFloor.pow(2).equals(x) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
197       case HALF_DOWN:
198       case HALF_UP:
199       case HALF_EVEN:
200         BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor);
201         /*
202          * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both
203          * x and halfSquare are integers, this is equivalent to testing whether or not x <=
204          * halfSquare.
205          */
206         return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
207       default:
208         throw new AssertionError();
209     }
210   }
211 
sqrtFloor(BigInteger x)212   private static BigInteger sqrtFloor(BigInteger x) {
213     /*
214      * Adapted from Hacker's Delight, Figure 11-1.
215      *
216      * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and
217      * then we can get a double approximation of the square root. Then, we iteratively improve this
218      * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2.
219      * This iteration has the following two properties:
220      *
221      * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is
222      * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean,
223      * and the arithmetic mean is always higher than the geometric mean.
224      *
225      * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles
226      * with each iteration, so this algorithm takes O(log(digits)) iterations.
227      *
228      * We start out with a double-precision approximation, which may be higher or lower than the
229      * true value. Therefore, we perform at least one Newton iteration to get a guess that's
230      * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point.
231      */
232     BigInteger sqrt0;
233     int log2 = log2(x, FLOOR);
234     if(log2 < DoubleUtils.MAX_DOUBLE_EXPONENT) {
235       sqrt0 = sqrtApproxWithDoubles(x);
236     } else {
237       int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even!
238       /*
239        * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be
240        * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift).
241        */
242       sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1);
243     }
244     BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
245     if (sqrt0.equals(sqrt1)) {
246       return sqrt0;
247     }
248     do {
249       sqrt0 = sqrt1;
250       sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
251     } while (sqrt1.compareTo(sqrt0) < 0);
252     return sqrt0;
253   }
254 
sqrtApproxWithDoubles(BigInteger x)255   private static BigInteger sqrtApproxWithDoubles(BigInteger x) {
256     return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN);
257   }
258 
259   /**
260    * Returns the result of dividing {@code p} by {@code q}, rounding using the specified
261    * {@code RoundingMode}.
262    *
263    * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
264    *         is not an integer multiple of {@code b}
265    */
divide(BigInteger p, BigInteger q, RoundingMode mode)266   public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode){
267     BigDecimal pDec = new BigDecimal(p);
268     BigDecimal qDec = new BigDecimal(q);
269     return pDec.divide(qDec, 0, mode).toBigIntegerExact();
270   }
271 
272   /**
273    * Returns {@code n!}, that is, the product of the first {@code n} positive
274    * integers, or {@code 1} if {@code n == 0}.
275    *
276    * <p><b>Warning</b>: the result takes <i>O(n log n)</i> space, so use cautiously.
277    *
278    * <p>This uses an efficient binary recursive algorithm to compute the factorial
279    * with balanced multiplies.  It also removes all the 2s from the intermediate
280    * products (shifting them back in at the end).
281    *
282    * @throws IllegalArgumentException if {@code n < 0}
283    */
factorial(int n)284   public static BigInteger factorial(int n) {
285     checkNonNegative("n", n);
286 
287     // If the factorial is small enough, just use LongMath to do it.
288     if (n < LongMath.FACTORIALS.length) {
289       return BigInteger.valueOf(LongMath.FACTORIALS[n]);
290     }
291 
292     // Pre-allocate space for our list of intermediate BigIntegers.
293     int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING);
294     ArrayList<BigInteger> bignums = new ArrayList<BigInteger>(approxSize);
295 
296     // Start from the pre-computed maximum long factorial.
297     int startingNumber = LongMath.FACTORIALS.length;
298     long product = LongMath.FACTORIALS[startingNumber - 1];
299     // Strip off 2s from this value.
300     int shift = Long.numberOfTrailingZeros(product);
301     product >>= shift;
302 
303     // Use floor(log2(num)) + 1 to prevent overflow of multiplication.
304     int productBits = LongMath.log2(product, FLOOR) + 1;
305     int bits = LongMath.log2(startingNumber, FLOOR) + 1;
306     // Check for the next power of two boundary, to save us a CLZ operation.
307     int nextPowerOfTwo = 1 << (bits - 1);
308 
309     // Iteratively multiply the longs as big as they can go.
310     for (long num = startingNumber; num <= n; num++) {
311       // Check to see if the floor(log2(num)) + 1 has changed.
312       if ((num & nextPowerOfTwo) != 0) {
313         nextPowerOfTwo <<= 1;
314         bits++;
315       }
316       // Get rid of the 2s in num.
317       int tz = Long.numberOfTrailingZeros(num);
318       long normalizedNum = num >> tz;
319       shift += tz;
320       // Adjust floor(log2(num)) + 1.
321       int normalizedBits = bits - tz;
322       // If it won't fit in a long, then we store off the intermediate product.
323       if (normalizedBits + productBits >= Long.SIZE) {
324         bignums.add(BigInteger.valueOf(product));
325         product = 1;
326         productBits = 0;
327       }
328       product *= normalizedNum;
329       productBits = LongMath.log2(product, FLOOR) + 1;
330     }
331     // Check for leftovers.
332     if (product > 1) {
333       bignums.add(BigInteger.valueOf(product));
334     }
335     // Efficiently multiply all the intermediate products together.
336     return listProduct(bignums).shiftLeft(shift);
337   }
338 
listProduct(List<BigInteger> nums)339   static BigInteger listProduct(List<BigInteger> nums) {
340     return listProduct(nums, 0, nums.size());
341   }
342 
listProduct(List<BigInteger> nums, int start, int end)343   static BigInteger listProduct(List<BigInteger> nums, int start, int end) {
344     switch (end - start) {
345       case 0:
346         return BigInteger.ONE;
347       case 1:
348         return nums.get(start);
349       case 2:
350         return nums.get(start).multiply(nums.get(start + 1));
351       case 3:
352         return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2));
353       default:
354         // Otherwise, split the list in half and recursively do this.
355         int m = (end + start) >>> 1;
356         return listProduct(nums, start, m).multiply(listProduct(nums, m, end));
357     }
358   }
359 
360  /**
361    * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
362    * {@code k}, that is, {@code n! / (k! (n - k)!)}.
363    *
364    * <p><b>Warning</b>: the result can take as much as <i>O(k log n)</i> space.
365    *
366    * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
367    */
binomial(int n, int k)368   public static BigInteger binomial(int n, int k) {
369     checkNonNegative("n", n);
370     checkNonNegative("k", k);
371     checkArgument(k <= n, "k (%s) > n (%s)", k, n);
372     if (k > (n >> 1)) {
373       k = n - k;
374     }
375     if (k < LongMath.BIGGEST_BINOMIALS.length && n <= LongMath.BIGGEST_BINOMIALS[k]) {
376       return BigInteger.valueOf(LongMath.binomial(n, k));
377     }
378     BigInteger result = BigInteger.ONE;
379     for (int i = 0; i < k; i++) {
380       result = result.multiply(BigInteger.valueOf(n - i));
381       result = result.divide(BigInteger.valueOf(i + 1));
382     }
383     return result;
384   }
385 
386   // Returns true if BigInteger.valueOf(x.longValue()).equals(x).
fitsInLong(BigInteger x)387   static boolean fitsInLong(BigInteger x) {
388     return x.bitLength() <= Long.SIZE - 1;
389   }
390 
BigIntegerMath()391   private BigIntegerMath() {}
392 }
393