    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    }