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