View Javadoc
1   /*
2    * Copyright (C) 2011 The Guava Authors
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
5    * in compliance with the License. You may obtain a copy of the License at
6    *
7    * http://www.apache.org/licenses/LICENSE-2.0
8    *
9    * Unless required by applicable law or agreed to in writing, software distributed under the License
10   * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
11   * or implied. See the License for the specific language governing permissions and limitations under
12   * the License.
13   */
14  
15  package com.google.common.math;
16  
17  import static com.google.common.base.Preconditions.checkArgument;
18  import static com.google.common.math.DoubleUtils.IMPLICIT_BIT;
19  import static com.google.common.math.DoubleUtils.SIGNIFICAND_BITS;
20  import static com.google.common.math.DoubleUtils.getSignificand;
21  import static com.google.common.math.DoubleUtils.isFinite;
22  import static com.google.common.math.DoubleUtils.isNormal;
23  import static com.google.common.math.DoubleUtils.scaleNormalize;
24  import static com.google.common.math.MathPreconditions.checkInRange;
25  import static com.google.common.math.MathPreconditions.checkNonNegative;
26  import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
27  import static java.lang.Math.abs;
28  import static java.lang.Math.copySign;
29  import static java.lang.Math.getExponent;
30  import static java.lang.Math.log;
31  import static java.lang.Math.rint;
32  
33  import com.google.common.annotations.GwtCompatible;
34  import com.google.common.annotations.GwtIncompatible;
35  import com.google.common.annotations.VisibleForTesting;
36  import com.google.common.primitives.Booleans;
37  import com.google.errorprone.annotations.CanIgnoreReturnValue;
38  import java.math.BigInteger;
39  import java.math.RoundingMode;
40  import java.util.Iterator;
41  
42  /**
43   * A class for arithmetic on doubles that is not covered by {@link java.lang.Math}.
44   *
45   * @author Louis Wasserman
46   * @since 11.0
47   */
48  @GwtCompatible(emulated = true)
49  public final class DoubleMath {
50    /*
51     * This method returns a value y such that rounding y DOWN (towards zero) gives the same result as
52     * rounding x according to the specified mode.
53     */
54    @GwtIncompatible // #isMathematicalInteger, com.google.common.math.DoubleUtils
55    static double roundIntermediate(double x, RoundingMode mode) {
56      if (!isFinite(x)) {
57        throw new ArithmeticException("input is infinite or NaN");
58      }
59      switch (mode) {
60        case UNNECESSARY:
61          checkRoundingUnnecessary(isMathematicalInteger(x));
62          return x;
63  
64        case FLOOR:
65          if (x >= 0.0 || isMathematicalInteger(x)) {
66            return x;
67          } else {
68            return (long) x - 1;
69          }
70  
71        case CEILING:
72          if (x <= 0.0 || isMathematicalInteger(x)) {
73            return x;
74          } else {
75            return (long) x + 1;
76          }
77  
78        case DOWN:
79          return x;
80  
81        case UP:
82          if (isMathematicalInteger(x)) {
83            return x;
84          } else {
85            return (long) x + (x > 0 ? 1 : -1);
86          }
87  
88        case HALF_EVEN:
89          return rint(x);
90  
91        case HALF_UP:
92          {
93            double z = rint(x);
94            if (abs(x - z) == 0.5) {
95              return x + copySign(0.5, x);
96            } else {
97              return z;
98            }
99          }
100 
101       case HALF_DOWN:
102         {
103           double z = rint(x);
104           if (abs(x - z) == 0.5) {
105             return x;
106           } else {
107             return z;
108           }
109         }
110 
111       default:
112         throw new AssertionError();
113     }
114   }
115 
116   /**
117    * Returns the {@code int} value that is equal to {@code x} rounded with the specified rounding
118    * mode, if possible.
119    *
120    * @throws ArithmeticException if
121    *     <ul>
122    *     <li>{@code x} is infinite or NaN
123    *     <li>{@code x}, after being rounded to a mathematical integer using the specified rounding
124    *         mode, is either less than {@code Integer.MIN_VALUE} or greater than {@code
125    *         Integer.MAX_VALUE}
126    *     <li>{@code x} is not a mathematical integer and {@code mode} is
127    *         {@link RoundingMode#UNNECESSARY}
128    *     </ul>
129    */
130   @GwtIncompatible // #roundIntermediate
131   public static int roundToInt(double x, RoundingMode mode) {
132     double z = roundIntermediate(x, mode);
133     checkInRange(z > MIN_INT_AS_DOUBLE - 1.0 & z < MAX_INT_AS_DOUBLE + 1.0);
134     return (int) z;
135   }
136 
137   private static final double MIN_INT_AS_DOUBLE = -0x1p31;
138   private static final double MAX_INT_AS_DOUBLE = 0x1p31 - 1.0;
139 
140   /**
141    * Returns the {@code long} value that is equal to {@code x} rounded with the specified rounding
142    * mode, if possible.
143    *
144    * @throws ArithmeticException if
145    *     <ul>
146    *     <li>{@code x} is infinite or NaN
147    *     <li>{@code x}, after being rounded to a mathematical integer using the specified rounding
148    *         mode, is either less than {@code Long.MIN_VALUE} or greater than {@code
149    *         Long.MAX_VALUE}
150    *     <li>{@code x} is not a mathematical integer and {@code mode} is
151    *         {@link RoundingMode#UNNECESSARY}
152    *     </ul>
153    */
154   @GwtIncompatible // #roundIntermediate
155   public static long roundToLong(double x, RoundingMode mode) {
156     double z = roundIntermediate(x, mode);
157     checkInRange(MIN_LONG_AS_DOUBLE - z < 1.0 & z < MAX_LONG_AS_DOUBLE_PLUS_ONE);
158     return (long) z;
159   }
160 
161   private static final double MIN_LONG_AS_DOUBLE = -0x1p63;
162   /*
163    * We cannot store Long.MAX_VALUE as a double without losing precision. Instead, we store
164    * Long.MAX_VALUE + 1 == -Long.MIN_VALUE, and then offset all comparisons by 1.
165    */
166   private static final double MAX_LONG_AS_DOUBLE_PLUS_ONE = 0x1p63;
167 
168   /**
169    * Returns the {@code BigInteger} value that is equal to {@code x} rounded with the specified
170    * rounding mode, if possible.
171    *
172    * @throws ArithmeticException if
173    *     <ul>
174    *     <li>{@code x} is infinite or NaN
175    *     <li>{@code x} is not a mathematical integer and {@code mode} is
176    *         {@link RoundingMode#UNNECESSARY}
177    *     </ul>
178    */
179   // #roundIntermediate, java.lang.Math.getExponent, com.google.common.math.DoubleUtils
180   @GwtIncompatible
181   public static BigInteger roundToBigInteger(double x, RoundingMode mode) {
182     x = roundIntermediate(x, mode);
183     if (MIN_LONG_AS_DOUBLE - x < 1.0 & x < MAX_LONG_AS_DOUBLE_PLUS_ONE) {
184       return BigInteger.valueOf((long) x);
185     }
186     int exponent = getExponent(x);
187     long significand = getSignificand(x);
188     BigInteger result = BigInteger.valueOf(significand).shiftLeft(exponent - SIGNIFICAND_BITS);
189     return (x < 0) ? result.negate() : result;
190   }
191 
192   /**
193    * Returns {@code true} if {@code x} is exactly equal to {@code 2^k} for some finite integer
194    * {@code k}.
195    */
196   @GwtIncompatible // com.google.common.math.DoubleUtils
197   public static boolean isPowerOfTwo(double x) {
198     if (x > 0.0 && isFinite(x)) {
199       long significand = getSignificand(x);
200       return (significand & (significand - 1)) == 0;
201     }
202     return false;
203   }
204 
205   /**
206    * Returns the base 2 logarithm of a double value.
207    *
208    * <p>Special cases:
209    * <ul>
210    * <li>If {@code x} is NaN or less than zero, the result is NaN.
211    * <li>If {@code x} is positive infinity, the result is positive infinity.
212    * <li>If {@code x} is positive or negative zero, the result is negative infinity.
213    * </ul>
214    *
215    * <p>The computed result is within 1 ulp of the exact result.
216    *
217    * <p>If the result of this method will be immediately rounded to an {@code int},
218    * {@link #log2(double, RoundingMode)} is faster.
219    */
220   public static double log2(double x) {
221     return log(x) / LN_2; // surprisingly within 1 ulp according to tests
222   }
223 
224   private static final double LN_2 = log(2);
225 
226   /**
227    * Returns the base 2 logarithm of a double value, rounded with the specified rounding mode to an
228    * {@code int}.
229    *
230    * <p>Regardless of the rounding mode, this is faster than {@code (int) log2(x)}.
231    *
232    * @throws IllegalArgumentException if {@code x <= 0.0}, {@code x} is NaN, or {@code x} is
233    *     infinite
234    */
235   @GwtIncompatible // java.lang.Math.getExponent, com.google.common.math.DoubleUtils
236   @SuppressWarnings("fallthrough")
237   public static int log2(double x, RoundingMode mode) {
238     checkArgument(x > 0.0 && isFinite(x), "x must be positive and finite");
239     int exponent = getExponent(x);
240     if (!isNormal(x)) {
241       return log2(x * IMPLICIT_BIT, mode) - SIGNIFICAND_BITS;
242       // Do the calculation on a normal value.
243     }
244     // x is positive, finite, and normal
245     boolean increment;
246     switch (mode) {
247       case UNNECESSARY:
248         checkRoundingUnnecessary(isPowerOfTwo(x));
249         // fall through
250       case FLOOR:
251         increment = false;
252         break;
253       case CEILING:
254         increment = !isPowerOfTwo(x);
255         break;
256       case DOWN:
257         increment = exponent < 0 & !isPowerOfTwo(x);
258         break;
259       case UP:
260         increment = exponent >= 0 & !isPowerOfTwo(x);
261         break;
262       case HALF_DOWN:
263       case HALF_EVEN:
264       case HALF_UP:
265         double xScaled = scaleNormalize(x);
266         // sqrt(2) is irrational, and the spec is relative to the "exact numerical result,"
267         // so log2(x) is never exactly exponent + 0.5.
268         increment = (xScaled * xScaled) > 2.0;
269         break;
270       default:
271         throw new AssertionError();
272     }
273     return increment ? exponent + 1 : exponent;
274   }
275 
276   /**
277    * Returns {@code true} if {@code x} represents a mathematical integer.
278    *
279    * <p>This is equivalent to, but not necessarily implemented as, the expression {@code
280    * !Double.isNaN(x) && !Double.isInfinite(x) && x == Math.rint(x)}.
281    */
282   @GwtIncompatible // java.lang.Math.getExponent, com.google.common.math.DoubleUtils
283   public static boolean isMathematicalInteger(double x) {
284     return isFinite(x)
285         && (x == 0.0
286             || SIGNIFICAND_BITS - Long.numberOfTrailingZeros(getSignificand(x)) <= getExponent(x));
287   }
288 
289   /**
290    * Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
291    * {@code n == 0}, or {@code n!}, or {@link Double#POSITIVE_INFINITY} if
292    * {@code n! > Double.MAX_VALUE}.
293    *
294    * <p>The result is within 1 ulp of the true value.
295    *
296    * @throws IllegalArgumentException if {@code n < 0}
297    */
298   public static double factorial(int n) {
299     checkNonNegative("n", n);
300     if (n > MAX_FACTORIAL) {
301       return Double.POSITIVE_INFINITY;
302     } else {
303       // Multiplying the last (n & 0xf) values into their own accumulator gives a more accurate
304       // result than multiplying by everySixteenthFactorial[n >> 4] directly.
305       double accum = 1.0;
306       for (int i = 1 + (n & ~0xf); i <= n; i++) {
307         accum *= i;
308       }
309       return accum * everySixteenthFactorial[n >> 4];
310     }
311   }
312 
313   @VisibleForTesting static final int MAX_FACTORIAL = 170;
314 
315   @VisibleForTesting
316   static final double[] everySixteenthFactorial = {
317     0x1.0p0,
318     0x1.30777758p44,
319     0x1.956ad0aae33a4p117,
320     0x1.ee69a78d72cb6p202,
321     0x1.fe478ee34844ap295,
322     0x1.c619094edabffp394,
323     0x1.3638dd7bd6347p498,
324     0x1.7cac197cfe503p605,
325     0x1.1e5dfc140e1e5p716,
326     0x1.8ce85fadb707ep829,
327     0x1.95d5f3d928edep945
328   };
329 
330   /**
331    * Returns {@code true} if {@code a} and {@code b} are within {@code tolerance} of each other.
332    *
333    * <p>Technically speaking, this is equivalent to {@code Math.abs(a - b) <= tolerance ||
334    * Double.valueOf(a).equals(Double.valueOf(b))}.
335    *
336    * <p>Notable special cases include:
337    *
338    * <ul>
339    *   <li>All NaNs are fuzzily equal.
340    *   <li>If {@code a == b}, then {@code a} and {@code b} are always fuzzily equal.
341    *   <li>Positive and negative zero are always fuzzily equal.
342    *   <li>If {@code tolerance} is zero, and neither {@code a} nor {@code b} is NaN, then {@code a}
343    *       and {@code b} are fuzzily equal if and only if {@code a == b}.
344    *   <li>With {@link Double#POSITIVE_INFINITY} tolerance, all non-NaN values are fuzzily equal.
345    *   <li>With finite tolerance, {@code Double.POSITIVE_INFINITY} and {@code
346    *       Double.NEGATIVE_INFINITY} are fuzzily equal only to themselves.
347    * </ul>
348    *
349    * <p>This is reflexive and symmetric, but <em>not</em> transitive, so it is <em>not</em> an
350    * equivalence relation and <em>not</em> suitable for use in {@link Object#equals}
351    * implementations.
352    *
353    * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
354    * @since 13.0
355    */
356   public static boolean fuzzyEquals(double a, double b, double tolerance) {
357     MathPreconditions.checkNonNegative("tolerance", tolerance);
358     return Math.copySign(a - b, 1.0) <= tolerance
359         // copySign(x, 1.0) is a branch-free version of abs(x), but with different NaN semantics
360         || (a == b) // needed to ensure that infinities equal themselves
361         || (Double.isNaN(a) && Double.isNaN(b));
362   }
363 
364   /**
365    * Compares {@code a} and {@code b} "fuzzily," with a tolerance for nearly-equal values.
366    *
367    * <p>This method is equivalent to
368    * {@code fuzzyEquals(a, b, tolerance) ? 0 : Double.compare(a, b)}. In particular, like
369    * {@link Double#compare(double, double)}, it treats all NaN values as equal and greater than all
370    * other values (including {@link Double#POSITIVE_INFINITY}).
371    *
372    * <p>This is <em>not</em> a total ordering and is <em>not</em> suitable for use in
373    * {@link Comparable#compareTo} implementations. In particular, it is not transitive.
374    *
375    * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
376    * @since 13.0
377    */
378   public static int fuzzyCompare(double a, double b, double tolerance) {
379     if (fuzzyEquals(a, b, tolerance)) {
380       return 0;
381     } else if (a < b) {
382       return -1;
383     } else if (a > b) {
384       return 1;
385     } else {
386       return Booleans.compare(Double.isNaN(a), Double.isNaN(b));
387     }
388   }
389 
390   /**
391    * Returns the <a href="http://en.wikipedia.org/wiki/Arithmetic_mean">arithmetic mean</a> of
392    * {@code values}.
393    *
394    * <p>If these values are a sample drawn from a population, this is also an unbiased estimator of
395    * the arithmetic mean of the population.
396    *
397    * @param values a nonempty series of values
398    * @throws IllegalArgumentException if {@code values} is empty or contains any non-finite value
399    * @deprecated Use {@link Stats#meanOf} instead, noting the less strict handling of non-finite
400    *     values.
401    */
402   @Deprecated
403   // com.google.common.math.DoubleUtils
404   @GwtIncompatible
405   public static double mean(double... values) {
406     checkArgument(values.length > 0, "Cannot take mean of 0 values");
407     long count = 1;
408     double mean = checkFinite(values[0]);
409     for (int index = 1; index < values.length; ++index) {
410       checkFinite(values[index]);
411       count++;
412       // Art of Computer Programming vol. 2, Knuth, 4.2.2, (15)
413       mean += (values[index] - mean) / count;
414     }
415     return mean;
416   }
417 
418   /**
419    * Returns the <a href="http://en.wikipedia.org/wiki/Arithmetic_mean">arithmetic mean</a> of
420    * {@code values}.
421    *
422    * <p>If these values are a sample drawn from a population, this is also an unbiased estimator of
423    * the arithmetic mean of the population.
424    *
425    * @param values a nonempty series of values
426    * @throws IllegalArgumentException if {@code values} is empty
427    * @deprecated Use {@link Stats#meanOf} instead, noting the less strict handling of non-finite
428    *     values.
429    */
430   @Deprecated
431   public static double mean(int... values) {
432     checkArgument(values.length > 0, "Cannot take mean of 0 values");
433     // The upper bound on the the length of an array and the bounds on the int values mean that, in
434     // this case only, we can compute the sum as a long without risking overflow or loss of
435     // precision. So we do that, as it's slightly quicker than the Knuth algorithm.
436     long sum = 0;
437     for (int index = 0; index < values.length; ++index) {
438       sum += values[index];
439     }
440     return (double) sum / values.length;
441   }
442 
443   /**
444    * Returns the <a href="http://en.wikipedia.org/wiki/Arithmetic_mean">arithmetic mean</a> of
445    * {@code values}.
446    *
447    * <p>If these values are a sample drawn from a population, this is also an unbiased estimator of
448    * the arithmetic mean of the population.
449    *
450    * @param values a nonempty series of values, which will be converted to {@code double} values
451    *     (this may cause loss of precision for longs of magnitude over 2^53 (slightly over 9e15))
452    * @throws IllegalArgumentException if {@code values} is empty
453    * @deprecated Use {@link Stats#meanOf} instead, noting the less strict handling of non-finite
454    *     values.
455    */
456   @Deprecated
457   public static double mean(long... values) {
458     checkArgument(values.length > 0, "Cannot take mean of 0 values");
459     long count = 1;
460     double mean = values[0];
461     for (int index = 1; index < values.length; ++index) {
462       count++;
463       // Art of Computer Programming vol. 2, Knuth, 4.2.2, (15)
464       mean += (values[index] - mean) / count;
465     }
466     return mean;
467   }
468 
469   /**
470    * Returns the <a href="http://en.wikipedia.org/wiki/Arithmetic_mean">arithmetic mean</a> of
471    * {@code values}.
472    *
473    * <p>If these values are a sample drawn from a population, this is also an unbiased estimator of
474    * the arithmetic mean of the population.
475    *
476    * @param values a nonempty series of values, which will be converted to {@code double} values
477    *     (this may cause loss of precision)
478    * @throws IllegalArgumentException if {@code values} is empty or contains any non-finite value
479    * @deprecated Use {@link Stats#meanOf} instead, noting the less strict handling of non-finite
480    *     values.
481    */
482   @Deprecated
483   // com.google.common.math.DoubleUtils
484   @GwtIncompatible
485   public static double mean(Iterable<? extends Number> values) {
486     return mean(values.iterator());
487   }
488 
489   /**
490    * Returns the <a href="http://en.wikipedia.org/wiki/Arithmetic_mean">arithmetic mean</a> of
491    * {@code values}.
492    *
493    * <p>If these values are a sample drawn from a population, this is also an unbiased estimator of
494    * the arithmetic mean of the population.
495    *
496    * @param values a nonempty series of values, which will be converted to {@code double} values
497    *     (this may cause loss of precision)
498    * @throws IllegalArgumentException if {@code values} is empty or contains any non-finite value
499    * @deprecated Use {@link Stats#meanOf} instead, noting the less strict handling of non-finite
500    *     values.
501    */
502   @Deprecated
503   // com.google.common.math.DoubleUtils
504   @GwtIncompatible
505   public static double mean(Iterator<? extends Number> values) {
506     checkArgument(values.hasNext(), "Cannot take mean of 0 values");
507     long count = 1;
508     double mean = checkFinite(values.next().doubleValue());
509     while (values.hasNext()) {
510       double value = checkFinite(values.next().doubleValue());
511       count++;
512       // Art of Computer Programming vol. 2, Knuth, 4.2.2, (15)
513       mean += (value - mean) / count;
514     }
515     return mean;
516   }
517 
518   @GwtIncompatible // com.google.common.math.DoubleUtils
519   @CanIgnoreReturnValue
520   private static double checkFinite(double argument) {
521     checkArgument(isFinite(argument));
522     return argument;
523   }
524 
525   private DoubleMath() {}
526 }