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>&lt; 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x &lt; y</li>
093         *       <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; 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>&lt; 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} &amp;&amp; x &lt; y</li>
119         *       <li>> 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} &amp;&amp; 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    }