001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.util;
019
020 import java.math.BigDecimal;
021
022 import org.apache.commons.math3.exception.MathArithmeticException;
023 import org.apache.commons.math3.exception.MathIllegalArgumentException;
024 import org.apache.commons.math3.exception.util.LocalizedFormats;
025
026 /**
027 * Utilities for comparing numbers.
028 *
029 * @since 3.0
030 * @version $Id: Precision.java 1422313 2012-12-15 18:53:41Z psteitz $
031 */
032 public class Precision {
033 /**
034 * <p>
035 * Largest double-precision floating-point number such that
036 * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
037 * bound on the relative error due to rounding real numbers to double
038 * precision floating-point numbers.
039 * </p>
040 * <p>
041 * In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
042 * </p>
043 *
044 * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
045 */
046 public static final double EPSILON;
047
048 /**
049 * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow.
050 * <br/>
051 * In IEEE 754 arithmetic, this is also the smallest normalized
052 * number 2<sup>-1022</sup>.
053 */
054 public static final double SAFE_MIN;
055
056 /** Exponent offset in IEEE754 representation. */
057 private static final long EXPONENT_OFFSET = 1023l;
058
059 /** Offset to order signed double numbers lexicographically. */
060 private static final long SGN_MASK = 0x8000000000000000L;
061 /** Offset to order signed double numbers lexicographically. */
062 private static final int SGN_MASK_FLOAT = 0x80000000;
063
064 static {
065 /*
066 * This was previously expressed as = 0x1.0p-53;
067 * However, OpenJDK (Sparc Solaris) cannot handle such small
068 * constants: MATH-721
069 */
070 EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53l) << 52);
071
072 /*
073 * This was previously expressed as = 0x1.0p-1022;
074 * However, OpenJDK (Sparc Solaris) cannot handle such small
075 * constants: MATH-721
076 */
077 SAFE_MIN = Double.longBitsToDouble((EXPONENT_OFFSET - 1022l) << 52);
078 }
079
080 /**
081 * Private constructor.
082 */
083 private Precision() {}
084
085 /**
086 * Compares two numbers given some amount of allowed error.
087 *
088 * @param x the first number
089 * @param y the second number
090 * @param eps the amount of error to allow when checking for equality
091 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li>
092 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li>
093 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul>
094 */
095 public static int compareTo(double x, double y, double eps) {
096 if (equals(x, y, eps)) {
097 return 0;
098 } else if (x < y) {
099 return -1;
100 }
101 return 1;
102 }
103
104 /**
105 * Compares two numbers given some amount of allowed error.
106 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
107 * (or fewer) floating point numbers between them, i.e. two adjacent floating
108 * point numbers are considered equal.
109 * Adapted from <a
110 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
111 * Bruce Dawson</a>
112 *
113 * @param x first value
114 * @param y second value
115 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
116 * values between {@code x} and {@code y}.
117 * @return <ul><li>0 if {@link #equals(double, double, int) equals(x, y, maxUlps)}</li>
118 * <li>< 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} && x < y</li>
119 * <li>> 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} && x > y</li></ul>
120 */
121 public static int compareTo(final double x, final double y, final int maxUlps) {
122 if (equals(x, y, maxUlps)) {
123 return 0;
124 } else if (x < y) {
125 return -1;
126 }
127 return 1;
128 }
129
130 /**
131 * Returns true iff they are equal as defined by
132 * {@link #equals(float,float,int) equals(x, y, 1)}.
133 *
134 * @param x first value
135 * @param y second value
136 * @return {@code true} if the values are equal.
137 */
138 public static boolean equals(float x, float y) {
139 return equals(x, y, 1);
140 }
141
142 /**
143 * Returns true if both arguments are NaN or neither is NaN and they are
144 * equal as defined by {@link #equals(float,float) equals(x, y, 1)}.
145 *
146 * @param x first value
147 * @param y second value
148 * @return {@code true} if the values are equal or both are NaN.
149 * @since 2.2
150 */
151 public static boolean equalsIncludingNaN(float x, float y) {
152 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, 1);
153 }
154
155 /**
156 * Returns true if both arguments are equal or within the range of allowed
157 * error (inclusive).
158 *
159 * @param x first value
160 * @param y second value
161 * @param eps the amount of absolute error to allow.
162 * @return {@code true} if the values are equal or within range of each other.
163 * @since 2.2
164 */
165 public static boolean equals(float x, float y, float eps) {
166 return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
167 }
168
169 /**
170 * Returns true if both arguments are NaN or are equal or within the range
171 * of allowed error (inclusive).
172 *
173 * @param x first value
174 * @param y second value
175 * @param eps the amount of absolute error to allow.
176 * @return {@code true} if the values are equal or within range of each other,
177 * or both are NaN.
178 * @since 2.2
179 */
180 public static boolean equalsIncludingNaN(float x, float y, float eps) {
181 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
182 }
183
184 /**
185 * Returns true if both arguments are equal or within the range of allowed
186 * error (inclusive).
187 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
188 * (or fewer) floating point numbers between them, i.e. two adjacent floating
189 * point numbers are considered equal.
190 * Adapted from <a
191 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
192 * Bruce Dawson</a>
193 *
194 * @param x first value
195 * @param y second value
196 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
197 * values between {@code x} and {@code y}.
198 * @return {@code true} if there are fewer than {@code maxUlps} floating
199 * point values between {@code x} and {@code y}.
200 * @since 2.2
201 */
202 public static boolean equals(float x, float y, int maxUlps) {
203 int xInt = Float.floatToIntBits(x);
204 int yInt = Float.floatToIntBits(y);
205
206 // Make lexicographically ordered as a two's-complement integer.
207 if (xInt < 0) {
208 xInt = SGN_MASK_FLOAT - xInt;
209 }
210 if (yInt < 0) {
211 yInt = SGN_MASK_FLOAT - yInt;
212 }
213
214 final boolean isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
215
216 return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
217 }
218
219 /**
220 * Returns true if both arguments are NaN or if they are equal as defined
221 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
222 *
223 * @param x first value
224 * @param y second value
225 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
226 * values between {@code x} and {@code y}.
227 * @return {@code true} if both arguments are NaN or if there are less than
228 * {@code maxUlps} floating point values between {@code x} and {@code y}.
229 * @since 2.2
230 */
231 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
232 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, maxUlps);
233 }
234
235 /**
236 * Returns true iff they are equal as defined by
237 * {@link #equals(double,double,int) equals(x, y, 1)}.
238 *
239 * @param x first value
240 * @param y second value
241 * @return {@code true} if the values are equal.
242 */
243 public static boolean equals(double x, double y) {
244 return equals(x, y, 1);
245 }
246
247 /**
248 * Returns true if both arguments are NaN or neither is NaN and they are
249 * equal as defined by {@link #equals(double,double) equals(x, y, 1)}.
250 *
251 * @param x first value
252 * @param y second value
253 * @return {@code true} if the values are equal or both are NaN.
254 * @since 2.2
255 */
256 public static boolean equalsIncludingNaN(double x, double y) {
257 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, 1);
258 }
259
260 /**
261 * Returns {@code true} if there is no double value strictly between the
262 * arguments or the difference between them is within the range of allowed
263 * error (inclusive).
264 *
265 * @param x First value.
266 * @param y Second value.
267 * @param eps Amount of allowed absolute error.
268 * @return {@code true} if the values are two adjacent floating point
269 * numbers or they are within range of each other.
270 */
271 public static boolean equals(double x, double y, double eps) {
272 return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
273 }
274
275 /**
276 * Returns {@code true} if there is no double value strictly between the
277 * arguments or the reltaive difference between them is smaller or equal
278 * to the given tolerance.
279 *
280 * @param x First value.
281 * @param y Second value.
282 * @param eps Amount of allowed relative error.
283 * @return {@code true} if the values are two adjacent floating point
284 * numbers or they are within range of each other.
285 * @since 3.1
286 */
287 public static boolean equalsWithRelativeTolerance(double x, double y, double eps) {
288 if (equals(x, y, 1)) {
289 return true;
290 }
291
292 final double absoluteMax = FastMath.max(FastMath.abs(x), FastMath.abs(y));
293 final double relativeDifference = FastMath.abs((x - y) / absoluteMax);
294
295 return relativeDifference <= eps;
296 }
297
298 /**
299 * Returns true if both arguments are NaN or are equal or within the range
300 * of allowed error (inclusive).
301 *
302 * @param x first value
303 * @param y second value
304 * @param eps the amount of absolute error to allow.
305 * @return {@code true} if the values are equal or within range of each other,
306 * or both are NaN.
307 * @since 2.2
308 */
309 public static boolean equalsIncludingNaN(double x, double y, double eps) {
310 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
311 }
312
313 /**
314 * Returns true if both arguments are equal or within the range of allowed
315 * error (inclusive).
316 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
317 * (or fewer) floating point numbers between them, i.e. two adjacent floating
318 * point numbers are considered equal.
319 * Adapted from <a
320 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
321 * Bruce Dawson</a>
322 *
323 * @param x first value
324 * @param y second value
325 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
326 * values between {@code x} and {@code y}.
327 * @return {@code true} if there are fewer than {@code maxUlps} floating
328 * point values between {@code x} and {@code y}.
329 */
330 public static boolean equals(double x, double y, int maxUlps) {
331 long xInt = Double.doubleToLongBits(x);
332 long yInt = Double.doubleToLongBits(y);
333
334 // Make lexicographically ordered as a two's-complement integer.
335 if (xInt < 0) {
336 xInt = SGN_MASK - xInt;
337 }
338 if (yInt < 0) {
339 yInt = SGN_MASK - yInt;
340 }
341
342 final boolean isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
343
344 return isEqual && !Double.isNaN(x) && !Double.isNaN(y);
345 }
346
347 /**
348 * Returns true if both arguments are NaN or if they are equal as defined
349 * by {@link #equals(double,double,int) equals(x, y, maxUlps)}.
350 *
351 * @param x first value
352 * @param y second value
353 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
354 * values between {@code x} and {@code y}.
355 * @return {@code true} if both arguments are NaN or if there are less than
356 * {@code maxUlps} floating point values between {@code x} and {@code y}.
357 * @since 2.2
358 */
359 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
360 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, maxUlps);
361 }
362
363 /**
364 * Rounds the given value to the specified number of decimal places.
365 * The value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
366 *
367 * @param x Value to round.
368 * @param scale Number of digits to the right of the decimal point.
369 * @return the rounded value.
370 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
371 */
372 public static double round(double x, int scale) {
373 return round(x, scale, BigDecimal.ROUND_HALF_UP);
374 }
375
376 /**
377 * Rounds the given value to the specified number of decimal places.
378 * The value is rounded using the given method which is any method defined
379 * in {@link BigDecimal}.
380 * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is
381 * returned unchanged, regardless of the other parameters.
382 *
383 * @param x Value to round.
384 * @param scale Number of digits to the right of the decimal point.
385 * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
386 * @return the rounded value.
387 * @throws ArithmeticException if {@code roundingMethod == ROUND_UNNECESSARY}
388 * and the specified scaling operation would require rounding.
389 * @throws IllegalArgumentException if {@code roundingMethod} does not
390 * represent a valid rounding mode.
391 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
392 */
393 public static double round(double x, int scale, int roundingMethod) {
394 try {
395 return (new BigDecimal
396 (Double.toString(x))
397 .setScale(scale, roundingMethod))
398 .doubleValue();
399 } catch (NumberFormatException ex) {
400 if (Double.isInfinite(x)) {
401 return x;
402 } else {
403 return Double.NaN;
404 }
405 }
406 }
407
408 /**
409 * Rounds the given value to the specified number of decimal places.
410 * The value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
411 *
412 * @param x Value to round.
413 * @param scale Number of digits to the right of the decimal point.
414 * @return the rounded value.
415 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
416 */
417 public static float round(float x, int scale) {
418 return round(x, scale, BigDecimal.ROUND_HALF_UP);
419 }
420
421 /**
422 * Rounds the given value to the specified number of decimal places.
423 * The value is rounded using the given method which is any method defined
424 * in {@link BigDecimal}.
425 *
426 * @param x Value to round.
427 * @param scale Number of digits to the right of the decimal point.
428 * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
429 * @return the rounded value.
430 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
431 * @throws MathArithmeticException if an exact operation is required but result is not exact
432 * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method.
433 */
434 public static float round(float x, int scale, int roundingMethod)
435 throws MathArithmeticException, MathIllegalArgumentException {
436 final float sign = FastMath.copySign(1f, x);
437 final float factor = (float) FastMath.pow(10.0f, scale) * sign;
438 return (float) roundUnscaled(x * factor, sign, roundingMethod) / factor;
439 }
440
441 /**
442 * Rounds the given non-negative value to the "nearest" integer. Nearest is
443 * determined by the rounding method specified. Rounding methods are defined
444 * in {@link BigDecimal}.
445 *
446 * @param unscaled Value to round.
447 * @param sign Sign of the original, scaled value.
448 * @param roundingMethod Rounding method, as defined in {@link BigDecimal}.
449 * @return the rounded value.
450 * @throws MathArithmeticException if an exact operation is required but result is not exact
451 * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method.
452 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
453 */
454 private static double roundUnscaled(double unscaled,
455 double sign,
456 int roundingMethod)
457 throws MathArithmeticException, MathIllegalArgumentException {
458 switch (roundingMethod) {
459 case BigDecimal.ROUND_CEILING :
460 if (sign == -1) {
461 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
462 } else {
463 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
464 }
465 break;
466 case BigDecimal.ROUND_DOWN :
467 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
468 break;
469 case BigDecimal.ROUND_FLOOR :
470 if (sign == -1) {
471 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
472 } else {
473 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
474 }
475 break;
476 case BigDecimal.ROUND_HALF_DOWN : {
477 unscaled = FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY);
478 double fraction = unscaled - FastMath.floor(unscaled);
479 if (fraction > 0.5) {
480 unscaled = FastMath.ceil(unscaled);
481 } else {
482 unscaled = FastMath.floor(unscaled);
483 }
484 break;
485 }
486 case BigDecimal.ROUND_HALF_EVEN : {
487 double fraction = unscaled - FastMath.floor(unscaled);
488 if (fraction > 0.5) {
489 unscaled = FastMath.ceil(unscaled);
490 } else if (fraction < 0.5) {
491 unscaled = FastMath.floor(unscaled);
492 } else {
493 // The following equality test is intentional and needed for rounding purposes
494 if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(Math
495 .floor(unscaled) / 2.0)) { // even
496 unscaled = FastMath.floor(unscaled);
497 } else { // odd
498 unscaled = FastMath.ceil(unscaled);
499 }
500 }
501 break;
502 }
503 case BigDecimal.ROUND_HALF_UP : {
504 unscaled = FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY);
505 double fraction = unscaled - FastMath.floor(unscaled);
506 if (fraction >= 0.5) {
507 unscaled = FastMath.ceil(unscaled);
508 } else {
509 unscaled = FastMath.floor(unscaled);
510 }
511 break;
512 }
513 case BigDecimal.ROUND_UNNECESSARY :
514 if (unscaled != FastMath.floor(unscaled)) {
515 throw new MathArithmeticException();
516 }
517 break;
518 case BigDecimal.ROUND_UP :
519 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
520 break;
521 default :
522 throw new MathIllegalArgumentException(LocalizedFormats.INVALID_ROUNDING_METHOD,
523 roundingMethod,
524 "ROUND_CEILING", BigDecimal.ROUND_CEILING,
525 "ROUND_DOWN", BigDecimal.ROUND_DOWN,
526 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
527 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
528 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
529 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
530 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
531 "ROUND_UP", BigDecimal.ROUND_UP);
532 }
533 return unscaled;
534 }
535
536
537 /**
538 * Computes a number {@code delta} close to {@code originalDelta} with
539 * the property that <pre><code>
540 * x + delta - x
541 * </code></pre>
542 * is exactly machine-representable.
543 * This is useful when computing numerical derivatives, in order to reduce
544 * roundoff errors.
545 *
546 * @param x Value.
547 * @param originalDelta Offset value.
548 * @return a number {@code delta} so that {@code x + delta} and {@code x}
549 * differ by a representable floating number.
550 */
551 public static double representableDelta(double x,
552 double originalDelta) {
553 return x + originalDelta - x;
554 }
555 }