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
017 package com.google.common.math;
018
019 import static com.google.common.base.Preconditions.checkArgument;
020 import static com.google.common.math.DoubleUtils.IMPLICIT_BIT;
021 import static com.google.common.math.DoubleUtils.SIGNIFICAND_BITS;
022 import static com.google.common.math.DoubleUtils.getSignificand;
023 import static com.google.common.math.DoubleUtils.isFinite;
024 import static com.google.common.math.DoubleUtils.isNormal;
025 import static com.google.common.math.DoubleUtils.scaleNormalize;
026 import static com.google.common.math.MathPreconditions.checkInRange;
027 import static com.google.common.math.MathPreconditions.checkNonNegative;
028 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
029 import static java.lang.Math.abs;
030 import static java.lang.Math.copySign;
031 import static java.lang.Math.getExponent;
032 import static java.lang.Math.log;
033 import static java.lang.Math.rint;
034
035 import com.google.common.annotations.Beta;
036 import com.google.common.annotations.VisibleForTesting;
037 import com.google.common.primitives.Booleans;
038
039 import java.math.BigInteger;
040 import java.math.RoundingMode;
041
042 /**
043 * A class for arithmetic on doubles that is not covered by {@link java.lang.Math}.
044 *
045 * @author Louis Wasserman
046 * @since 11.0
047 */
048 @Beta
049 public final class DoubleMath {
050 /*
051 * This method returns a value y such that rounding y DOWN (towards zero) gives the same result
052 * as rounding x according to the specified mode.
053 */
054 static double roundIntermediate(double x, RoundingMode mode) {
055 if (!isFinite(x)) {
056 throw new ArithmeticException("input is infinite or NaN");
057 }
058 switch (mode) {
059 case UNNECESSARY:
060 checkRoundingUnnecessary(isMathematicalInteger(x));
061 return x;
062
063 case FLOOR:
064 if (x >= 0.0 || isMathematicalInteger(x)) {
065 return x;
066 } else {
067 return x - 1.0;
068 }
069
070 case CEILING:
071 if (x <= 0.0 || isMathematicalInteger(x)) {
072 return x;
073 } else {
074 return x + 1.0;
075 }
076
077 case DOWN:
078 return x;
079
080 case UP:
081 if (isMathematicalInteger(x)) {
082 return x;
083 } else {
084 return x + Math.copySign(1.0, x);
085 }
086
087 case HALF_EVEN:
088 return rint(x);
089
090 case HALF_UP: {
091 double z = rint(x);
092 if (abs(x - z) == 0.5) {
093 return x + copySign(0.5, x);
094 } else {
095 return z;
096 }
097 }
098
099 case HALF_DOWN: {
100 double z = rint(x);
101 if (abs(x - z) == 0.5) {
102 return x;
103 } else {
104 return z;
105 }
106 }
107
108 default:
109 throw new AssertionError();
110 }
111 }
112
113 /**
114 * Returns the {@code int} value that is equal to {@code x} rounded with the specified rounding
115 * mode, if possible.
116 *
117 * @throws ArithmeticException if
118 * <ul>
119 * <li>{@code x} is infinite or NaN
120 * <li>{@code x}, after being rounded to a mathematical integer using the specified
121 * rounding mode, is either less than {@code Integer.MIN_VALUE} or greater than {@code
122 * Integer.MAX_VALUE}
123 * <li>{@code x} is not a mathematical integer and {@code mode} is
124 * {@link RoundingMode#UNNECESSARY}
125 * </ul>
126 */
127 public static int roundToInt(double x, RoundingMode mode) {
128 double z = roundIntermediate(x, mode);
129 checkInRange(z > MIN_INT_AS_DOUBLE - 1.0 & z < MAX_INT_AS_DOUBLE + 1.0);
130 return (int) z;
131 }
132
133 private static final double MIN_INT_AS_DOUBLE = -0x1p31;
134 private static final double MAX_INT_AS_DOUBLE = 0x1p31 - 1.0;
135
136 /**
137 * Returns the {@code long} value that is equal to {@code x} rounded with the specified rounding
138 * mode, if possible.
139 *
140 * @throws ArithmeticException if
141 * <ul>
142 * <li>{@code x} is infinite or NaN
143 * <li>{@code x}, after being rounded to a mathematical integer using the specified
144 * rounding mode, is either less than {@code Long.MIN_VALUE} or greater than {@code
145 * Long.MAX_VALUE}
146 * <li>{@code x} is not a mathematical integer and {@code mode} is
147 * {@link RoundingMode#UNNECESSARY}
148 * </ul>
149 */
150 public static long roundToLong(double x, RoundingMode mode) {
151 double z = roundIntermediate(x, mode);
152 checkInRange(MIN_LONG_AS_DOUBLE - z < 1.0 & z < MAX_LONG_AS_DOUBLE_PLUS_ONE);
153 return (long) z;
154 }
155
156 private static final double MIN_LONG_AS_DOUBLE = -0x1p63;
157 /*
158 * We cannot store Long.MAX_VALUE as a double without losing precision. Instead, we store
159 * Long.MAX_VALUE + 1 == -Long.MIN_VALUE, and then offset all comparisons by 1.
160 */
161 private static final double MAX_LONG_AS_DOUBLE_PLUS_ONE = 0x1p63;
162
163 /**
164 * Returns the {@code BigInteger} value that is equal to {@code x} rounded with the specified
165 * rounding mode, if possible.
166 *
167 * @throws ArithmeticException if
168 * <ul>
169 * <li>{@code x} is infinite or NaN
170 * <li>{@code x} is not a mathematical integer and {@code mode} is
171 * {@link RoundingMode#UNNECESSARY}
172 * </ul>
173 */
174 public static BigInteger roundToBigInteger(double x, RoundingMode mode) {
175 x = roundIntermediate(x, mode);
176 if (MIN_LONG_AS_DOUBLE - x < 1.0 & x < MAX_LONG_AS_DOUBLE_PLUS_ONE) {
177 return BigInteger.valueOf((long) x);
178 }
179 int exponent = getExponent(x);
180 long significand = getSignificand(x);
181 BigInteger result = BigInteger.valueOf(significand).shiftLeft(exponent - SIGNIFICAND_BITS);
182 return (x < 0) ? result.negate() : result;
183 }
184
185 /**
186 * Returns {@code true} if {@code x} is exactly equal to {@code 2^k} for some finite integer
187 * {@code k}.
188 */
189 public static boolean isPowerOfTwo(double x) {
190 return x > 0.0 && isFinite(x) && LongMath.isPowerOfTwo(getSignificand(x));
191 }
192
193 /**
194 * Returns the base 2 logarithm of a double value.
195 *
196 * <p>Special cases:
197 * <ul>
198 * <li>If {@code x} is NaN or less than zero, the result is NaN.
199 * <li>If {@code x} is positive infinity, the result is positive infinity.
200 * <li>If {@code x} is positive or negative zero, the result is negative infinity.
201 * </ul>
202 *
203 * <p>The computed result is within 1 ulp of the exact result.
204 *
205 * <p>If the result of this method will be immediately rounded to an {@code int},
206 * {@link #log2(double, RoundingMode)} is faster.
207 */
208 public static double log2(double x) {
209 return log(x) / LN_2; // surprisingly within 1 ulp according to tests
210 }
211
212 private static final double LN_2 = log(2);
213
214 /**
215 * Returns the base 2 logarithm of a double value, rounded with the specified rounding mode to an
216 * {@code int}.
217 *
218 * <p>Regardless of the rounding mode, this is faster than {@code (int) log2(x)}.
219 *
220 * @throws IllegalArgumentException if {@code x <= 0.0}, {@code x} is NaN, or {@code x} is
221 * infinite
222 */
223 @SuppressWarnings("fallthrough")
224 public static int log2(double x, RoundingMode mode) {
225 checkArgument(x > 0.0 && isFinite(x), "x must be positive and finite");
226 int exponent = getExponent(x);
227 if (!isNormal(x)) {
228 return log2(x * IMPLICIT_BIT, mode) - SIGNIFICAND_BITS;
229 // Do the calculation on a normal value.
230 }
231 // x is positive, finite, and normal
232 boolean increment;
233 switch (mode) {
234 case UNNECESSARY:
235 checkRoundingUnnecessary(isPowerOfTwo(x));
236 // fall through
237 case FLOOR:
238 increment = false;
239 break;
240 case CEILING:
241 increment = !isPowerOfTwo(x);
242 break;
243 case DOWN:
244 increment = exponent < 0 & !isPowerOfTwo(x);
245 break;
246 case UP:
247 increment = exponent >= 0 & !isPowerOfTwo(x);
248 break;
249 case HALF_DOWN:
250 case HALF_EVEN:
251 case HALF_UP:
252 double xScaled = scaleNormalize(x);
253 // sqrt(2) is irrational, and the spec is relative to the "exact numerical result,"
254 // so log2(x) is never exactly exponent + 0.5.
255 increment = (xScaled * xScaled) > 2.0;
256 break;
257 default:
258 throw new AssertionError();
259 }
260 return increment ? exponent + 1 : exponent;
261 }
262
263 /**
264 * Returns {@code true} if {@code x} represents a mathematical integer.
265 *
266 * <p>This is equivalent to, but not necessarily implemented as, the expression {@code
267 * !Double.isNaN(x) && !Double.isInfinite(x) && x == Math.rint(x)}.
268 */
269 public static boolean isMathematicalInteger(double x) {
270 return isFinite(x)
271 && (x == 0.0 ||
272 SIGNIFICAND_BITS - Long.numberOfTrailingZeros(getSignificand(x)) <= getExponent(x));
273 }
274
275 /**
276 * Returns {@code n!}, that is, the product of the first {@code n} positive
277 * integers, {@code 1} if {@code n == 0}, or e n!}, or
278 * {@link Double#POSITIVE_INFINITY} if {@code n! > Double.MAX_VALUE}.
279 *
280 * <p>The result is within 1 ulp of the true value.
281 *
282 * @throws IllegalArgumentException if {@code n < 0}
283 */
284 public static double factorial(int n) {
285 checkNonNegative("n", n);
286 if (n > MAX_FACTORIAL) {
287 return Double.POSITIVE_INFINITY;
288 } else {
289 // Multiplying the last (n & 0xf) values into their own accumulator gives a more accurate
290 // result than multiplying by EVERY_SIXTEENTH_FACTORIAL[n >> 4] directly.
291 double accum = 1.0;
292 for (int i = 1 + (n & ~0xf); i <= n; i++) {
293 accum *= i;
294 }
295 return accum * EVERY_SIXTEENTH_FACTORIAL[n >> 4];
296 }
297 }
298
299 @VisibleForTesting
300 static final int MAX_FACTORIAL = 170;
301
302 @VisibleForTesting
303 static final double[] EVERY_SIXTEENTH_FACTORIAL = {
304 0x1.0p0,
305 0x1.30777758p44,
306 0x1.956ad0aae33a4p117,
307 0x1.ee69a78d72cb6p202,
308 0x1.fe478ee34844ap295,
309 0x1.c619094edabffp394,
310 0x1.3638dd7bd6347p498,
311 0x1.7cac197cfe503p605,
312 0x1.1e5dfc140e1e5p716,
313 0x1.8ce85fadb707ep829,
314 0x1.95d5f3d928edep945};
315
316 /**
317 * Returns {@code true} if {@code a} and {@code b} are within {@code tolerance} of each other.
318 *
319 * <p>Technically speaking, this is equivalent to
320 * {@code Math.abs(a - b) <= tolerance || Double.valueOf(a).equals(Double.valueOf(b))}.
321 *
322 * <p>Notable special cases include:
323 * <ul>
324 * <li>All NaNs are fuzzily equal.
325 * <li>If {@code a == b}, then {@code a} and {@code b} are always fuzzily equal.
326 * <li>Positive and negative zero are always fuzzily equal.
327 * <li>If {@code tolerance} is zero, and neither {@code a} nor {@code b} is NaN, then
328 * {@code a} and {@code b} are fuzzily equal if and only if {@code a == b}.
329 * <li>With {@link Double#POSITIVE_INFINITY} tolerance, all non-NaN values are fuzzily equal.
330 * <li>With finite tolerance, {@code Double.POSITIVE_INFINITY} and {@code
331 * Double.NEGATIVE_INFINITY} are fuzzily equal only to themselves.
332 * </li>
333 *
334 * <p>This is reflexive and symmetric, but <em>not</em> transitive, so it is <em>not</em> an
335 * equivalence relation and <em>not</em> suitable for use in {@link Object#equals}
336 * implementations.
337 *
338 * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
339 * @since 13.0
340 */
341 @Beta
342 public static boolean fuzzyEquals(double a, double b, double tolerance) {
343 MathPreconditions.checkNonNegative("tolerance", tolerance);
344 return
345 Math.copySign(a - b, 1.0) <= tolerance
346 // copySign(x, 1.0) is a branch-free version of abs(x), but with different NaN semantics
347 || (a == b) // needed to ensure that infinities equal themselves
348 || ((a != a) && (b != b)); // x != x is equivalent to Double.isNaN(x), but faster
349 }
350
351 /**
352 * Compares {@code a} and {@code b} "fuzzily," with a tolerance for nearly-equal values.
353 *
354 * <p>This method is equivalent to
355 * {@code fuzzyEquals(a, b, tolerance) ? 0 : Double.compare(a, b)}. In particular, like
356 * {@link Double#compare(double, double)}, it treats all NaN values as equal and greater than all
357 * other values (including {@link Double#POSITIVE_INFINITY}).
358 *
359 * <p>This is <em>not</em> a total ordering and is <em>not</em> suitable for use in
360 * {@link Comparable#compareTo} implementations. In particular, it is not transitive.
361 *
362 * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
363 * @since 13.0
364 */
365 @Beta
366 public static int fuzzyCompare(double a, double b, double tolerance) {
367 if (fuzzyEquals(a, b, tolerance)) {
368 return 0;
369 } else if (a < b) {
370 return -1;
371 } else if (a > b) {
372 return 1;
373 } else {
374 return Booleans.compare(Double.isNaN(a), Double.isNaN(b));
375 }
376 }
377
378 private DoubleMath() {}
379 }