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
018package org.apache.commons.math3.util;
019
020import java.math.BigDecimal;
021
022import org.apache.commons.math3.exception.MathArithmeticException;
023import org.apache.commons.math3.exception.MathIllegalArgumentException;
024import 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 */
032public 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}