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 */
031public class Precision {
032    /**
033     * <p>
034     * Largest double-precision floating-point number such that
035     * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
036     * bound on the relative error due to rounding real numbers to double
037     * precision floating-point numbers.
038     * </p>
039     * <p>
040     * In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
041     * </p>
042     *
043     * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
044     */
045    public static final double EPSILON;
046
047    /**
048     * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow.
049     * <br/>
050     * In IEEE 754 arithmetic, this is also the smallest normalized
051     * number 2<sup>-1022</sup>.
052     */
053    public static final double SAFE_MIN;
054
055    /** Exponent offset in IEEE754 representation. */
056    private static final long EXPONENT_OFFSET = 1023l;
057
058    /** Offset to order signed double numbers lexicographically. */
059    private static final long SGN_MASK = 0x8000000000000000L;
060    /** Offset to order signed double numbers lexicographically. */
061    private static final int SGN_MASK_FLOAT = 0x80000000;
062    /** Positive zero. */
063    private static final double POSITIVE_ZERO = 0d;
064    /** Positive zero bits. */
065    private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0);
066    /** Negative zero bits. */
067    private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0);
068    /** Positive zero bits. */
069    private static final int POSITIVE_ZERO_FLOAT_BITS   = Float.floatToRawIntBits(+0.0f);
070    /** Negative zero bits. */
071    private static final int NEGATIVE_ZERO_FLOAT_BITS   = Float.floatToRawIntBits(-0.0f);
072
073    static {
074        /*
075         *  This was previously expressed as = 0x1.0p-53;
076         *  However, OpenJDK (Sparc Solaris) cannot handle such small
077         *  constants: MATH-721
078         */
079        EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53l) << 52);
080
081        /*
082         * This was previously expressed as = 0x1.0p-1022;
083         * However, OpenJDK (Sparc Solaris) cannot handle such small
084         * constants: MATH-721
085         */
086        SAFE_MIN = Double.longBitsToDouble((EXPONENT_OFFSET - 1022l) << 52);
087    }
088
089    /**
090     * Private constructor.
091     */
092    private Precision() {}
093
094    /**
095     * Compares two numbers given some amount of allowed error.
096     *
097     * @param x the first number
098     * @param y the second number
099     * @param eps the amount of error to allow when checking for equality
100     * @return <ul><li>0 if  {@link #equals(double, double, double) equals(x, y, eps)}</li>
101     *       <li>&lt; 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x &lt; y</li>
102     *       <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x > y</li></ul>
103     */
104    public static int compareTo(double x, double y, double eps) {
105        if (equals(x, y, eps)) {
106            return 0;
107        } else if (x < y) {
108            return -1;
109        }
110        return 1;
111    }
112
113    /**
114     * Compares two numbers given some amount of allowed error.
115     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
116     * (or fewer) floating point numbers between them, i.e. two adjacent floating
117     * point numbers are considered equal.
118     * Adapted from <a
119     * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
120     * Bruce Dawson</a>
121     *
122     * @param x first value
123     * @param y second value
124     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
125     * values between {@code x} and {@code y}.
126     * @return <ul><li>0 if  {@link #equals(double, double, int) equals(x, y, maxUlps)}</li>
127     *       <li>&lt; 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} &amp;&amp; x &lt; y</li>
128     *       <li>> 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} &amp;&amp; x > y</li></ul>
129     */
130    public static int compareTo(final double x, final double y, final int maxUlps) {
131        if (equals(x, y, maxUlps)) {
132            return 0;
133        } else if (x < y) {
134            return -1;
135        }
136        return 1;
137    }
138
139    /**
140     * Returns true iff they are equal as defined by
141     * {@link #equals(float,float,int) equals(x, y, 1)}.
142     *
143     * @param x first value
144     * @param y second value
145     * @return {@code true} if the values are equal.
146     */
147    public static boolean equals(float x, float y) {
148        return equals(x, y, 1);
149    }
150
151    /**
152     * Returns true if both arguments are NaN or neither is NaN and they are
153     * equal as defined by {@link #equals(float,float) equals(x, y, 1)}.
154     *
155     * @param x first value
156     * @param y second value
157     * @return {@code true} if the values are equal or both are NaN.
158     * @since 2.2
159     */
160    public static boolean equalsIncludingNaN(float x, float y) {
161        return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, 1);
162    }
163
164    /**
165     * Returns true if both arguments are equal or within the range of allowed
166     * error (inclusive).
167     *
168     * @param x first value
169     * @param y second value
170     * @param eps the amount of absolute error to allow.
171     * @return {@code true} if the values are equal or within range of each other.
172     * @since 2.2
173     */
174    public static boolean equals(float x, float y, float eps) {
175        return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
176    }
177
178    /**
179     * Returns true if both arguments are NaN or are equal or within the range
180     * of allowed error (inclusive).
181     *
182     * @param x first value
183     * @param y second value
184     * @param eps the amount of absolute error to allow.
185     * @return {@code true} if the values are equal or within range of each other,
186     * or both are NaN.
187     * @since 2.2
188     */
189    public static boolean equalsIncludingNaN(float x, float y, float eps) {
190        return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
191    }
192
193    /**
194     * Returns true if both arguments are equal or within the range of allowed
195     * error (inclusive).
196     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
197     * (or fewer) floating point numbers between them, i.e. two adjacent floating
198     * point numbers are considered equal.
199     * Adapted from <a
200     * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
201     * Bruce Dawson</a>
202     *
203     * @param x first value
204     * @param y second value
205     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
206     * values between {@code x} and {@code y}.
207     * @return {@code true} if there are fewer than {@code maxUlps} floating
208     * point values between {@code x} and {@code y}.
209     * @since 2.2
210     */
211    public static boolean equals(final float x, final float y, final int maxUlps) {
212
213        final int xInt = Float.floatToRawIntBits(x);
214        final int yInt = Float.floatToRawIntBits(y);
215
216        final boolean isEqual;
217        if (((xInt ^ yInt) & SGN_MASK_FLOAT) == 0) {
218            // number have same sign, there is no risk of overflow
219            isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
220        } else {
221            // number have opposite signs, take care of overflow
222            final int deltaPlus;
223            final int deltaMinus;
224            if (xInt < yInt) {
225                deltaPlus  = yInt - POSITIVE_ZERO_FLOAT_BITS;
226                deltaMinus = xInt - NEGATIVE_ZERO_FLOAT_BITS;
227            } else {
228                deltaPlus  = xInt - POSITIVE_ZERO_FLOAT_BITS;
229                deltaMinus = yInt - NEGATIVE_ZERO_FLOAT_BITS;
230            }
231
232            if (deltaPlus > maxUlps) {
233                isEqual = false;
234            } else {
235                isEqual = deltaMinus <= (maxUlps - deltaPlus);
236            }
237
238        }
239
240        return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
241
242    }
243
244    /**
245     * Returns true if both arguments are NaN or if they are equal as defined
246     * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
247     *
248     * @param x first value
249     * @param y second value
250     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
251     * values between {@code x} and {@code y}.
252     * @return {@code true} if both arguments are NaN or if there are less than
253     * {@code maxUlps} floating point values between {@code x} and {@code y}.
254     * @since 2.2
255     */
256    public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
257        return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, maxUlps);
258    }
259
260    /**
261     * Returns true iff they are equal as defined by
262     * {@link #equals(double,double,int) equals(x, y, 1)}.
263     *
264     * @param x first value
265     * @param y second value
266     * @return {@code true} if the values are equal.
267     */
268    public static boolean equals(double x, double y) {
269        return equals(x, y, 1);
270    }
271
272    /**
273     * Returns true if both arguments are NaN or neither is NaN and they are
274     * equal as defined by {@link #equals(double,double) equals(x, y, 1)}.
275     *
276     * @param x first value
277     * @param y second value
278     * @return {@code true} if the values are equal or both are NaN.
279     * @since 2.2
280     */
281    public static boolean equalsIncludingNaN(double x, double y) {
282        return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, 1);
283    }
284
285    /**
286     * Returns {@code true} if there is no double value strictly between the
287     * arguments or the difference between them is within the range of allowed
288     * error (inclusive).
289     *
290     * @param x First value.
291     * @param y Second value.
292     * @param eps Amount of allowed absolute error.
293     * @return {@code true} if the values are two adjacent floating point
294     * numbers or they are within range of each other.
295     */
296    public static boolean equals(double x, double y, double eps) {
297        return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
298    }
299
300    /**
301     * Returns {@code true} if there is no double value strictly between the
302     * arguments or the relative difference between them is smaller or equal
303     * to the given tolerance.
304     *
305     * @param x First value.
306     * @param y Second value.
307     * @param eps Amount of allowed relative error.
308     * @return {@code true} if the values are two adjacent floating point
309     * numbers or they are within range of each other.
310     * @since 3.1
311     */
312    public static boolean equalsWithRelativeTolerance(double x, double y, double eps) {
313        if (equals(x, y, 1)) {
314            return true;
315        }
316
317        final double absoluteMax = FastMath.max(FastMath.abs(x), FastMath.abs(y));
318        final double relativeDifference = FastMath.abs((x - y) / absoluteMax);
319
320        return relativeDifference <= eps;
321    }
322
323    /**
324     * Returns true if both arguments are NaN or are equal or within the range
325     * of allowed error (inclusive).
326     *
327     * @param x first value
328     * @param y second value
329     * @param eps the amount of absolute error to allow.
330     * @return {@code true} if the values are equal or within range of each other,
331     * or both are NaN.
332     * @since 2.2
333     */
334    public static boolean equalsIncludingNaN(double x, double y, double eps) {
335        return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
336    }
337
338    /**
339     * Returns true if both arguments are equal or within the range of allowed
340     * error (inclusive).
341     * <p>
342     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
343     * (or fewer) floating point numbers between them, i.e. two adjacent
344     * floating point numbers are considered equal.
345     * </p>
346     * <p>
347     * Adapted from <a
348     * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
349     * Bruce Dawson</a>
350     * </p>
351     *
352     * @param x first value
353     * @param y second value
354     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
355     * values between {@code x} and {@code y}.
356     * @return {@code true} if there are fewer than {@code maxUlps} floating
357     * point values between {@code x} and {@code y}.
358     */
359    public static boolean equals(final double x, final double y, final int maxUlps) {
360
361        final long xInt = Double.doubleToRawLongBits(x);
362        final long yInt = Double.doubleToRawLongBits(y);
363
364        final boolean isEqual;
365        if (((xInt ^ yInt) & SGN_MASK) == 0l) {
366            // number have same sign, there is no risk of overflow
367            isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
368        } else {
369            // number have opposite signs, take care of overflow
370            final long deltaPlus;
371            final long deltaMinus;
372            if (xInt < yInt) {
373                deltaPlus  = yInt - POSITIVE_ZERO_DOUBLE_BITS;
374                deltaMinus = xInt - NEGATIVE_ZERO_DOUBLE_BITS;
375            } else {
376                deltaPlus  = xInt - POSITIVE_ZERO_DOUBLE_BITS;
377                deltaMinus = yInt - NEGATIVE_ZERO_DOUBLE_BITS;
378            }
379
380            if (deltaPlus > maxUlps) {
381                isEqual = false;
382            } else {
383                isEqual = deltaMinus <= (maxUlps - deltaPlus);
384            }
385
386        }
387
388        return isEqual && !Double.isNaN(x) && !Double.isNaN(y);
389
390    }
391
392    /**
393     * Returns true if both arguments are NaN or if they are equal as defined
394     * by {@link #equals(double,double,int) equals(x, y, maxUlps)}.
395     *
396     * @param x first value
397     * @param y second value
398     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
399     * values between {@code x} and {@code y}.
400     * @return {@code true} if both arguments are NaN or if there are less than
401     * {@code maxUlps} floating point values between {@code x} and {@code y}.
402     * @since 2.2
403     */
404    public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
405        return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, maxUlps);
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 double round(double 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     * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is
426     * returned unchanged, regardless of the other parameters.
427     *
428     * @param x Value to round.
429     * @param scale Number of digits to the right of the decimal point.
430     * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
431     * @return the rounded value.
432     * @throws ArithmeticException if {@code roundingMethod == ROUND_UNNECESSARY}
433     * and the specified scaling operation would require rounding.
434     * @throws IllegalArgumentException if {@code roundingMethod} does not
435     * represent a valid rounding mode.
436     * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
437     */
438    public static double round(double x, int scale, int roundingMethod) {
439        try {
440            final double rounded = (new BigDecimal(Double.toString(x))
441                   .setScale(scale, roundingMethod))
442                   .doubleValue();
443            // MATH-1089: negative values rounded to zero should result in negative zero
444            return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded;
445        } catch (NumberFormatException ex) {
446            if (Double.isInfinite(x)) {
447                return x;
448            } else {
449                return Double.NaN;
450            }
451        }
452    }
453
454    /**
455     * Rounds the given value to the specified number of decimal places.
456     * The value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
457     *
458     * @param x Value to round.
459     * @param scale Number of digits to the right of the decimal point.
460     * @return the rounded value.
461     * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
462     */
463    public static float round(float x, int scale) {
464        return round(x, scale, BigDecimal.ROUND_HALF_UP);
465    }
466
467    /**
468     * Rounds the given value to the specified number of decimal places.
469     * The value is rounded using the given method which is any method defined
470     * in {@link BigDecimal}.
471     *
472     * @param x Value to round.
473     * @param scale Number of digits to the right of the decimal point.
474     * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
475     * @return the rounded value.
476     * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
477     * @throws MathArithmeticException if an exact operation is required but result is not exact
478     * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method.
479     */
480    public static float round(float x, int scale, int roundingMethod)
481        throws MathArithmeticException, MathIllegalArgumentException {
482        final float sign = FastMath.copySign(1f, x);
483        final float factor = (float) FastMath.pow(10.0f, scale) * sign;
484        return (float) roundUnscaled(x * factor, sign, roundingMethod) / factor;
485    }
486
487    /**
488     * Rounds the given non-negative value to the "nearest" integer. Nearest is
489     * determined by the rounding method specified. Rounding methods are defined
490     * in {@link BigDecimal}.
491     *
492     * @param unscaled Value to round.
493     * @param sign Sign of the original, scaled value.
494     * @param roundingMethod Rounding method, as defined in {@link BigDecimal}.
495     * @return the rounded value.
496     * @throws MathArithmeticException if an exact operation is required but result is not exact
497     * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method.
498     * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0)
499     */
500    private static double roundUnscaled(double unscaled,
501                                        double sign,
502                                        int roundingMethod)
503        throws MathArithmeticException, MathIllegalArgumentException {
504        switch (roundingMethod) {
505        case BigDecimal.ROUND_CEILING :
506            if (sign == -1) {
507                unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
508            } else {
509                unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
510            }
511            break;
512        case BigDecimal.ROUND_DOWN :
513            unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
514            break;
515        case BigDecimal.ROUND_FLOOR :
516            if (sign == -1) {
517                unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
518            } else {
519                unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
520            }
521            break;
522        case BigDecimal.ROUND_HALF_DOWN : {
523            unscaled = FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY);
524            double fraction = unscaled - FastMath.floor(unscaled);
525            if (fraction > 0.5) {
526                unscaled = FastMath.ceil(unscaled);
527            } else {
528                unscaled = FastMath.floor(unscaled);
529            }
530            break;
531        }
532        case BigDecimal.ROUND_HALF_EVEN : {
533            double fraction = unscaled - FastMath.floor(unscaled);
534            if (fraction > 0.5) {
535                unscaled = FastMath.ceil(unscaled);
536            } else if (fraction < 0.5) {
537                unscaled = FastMath.floor(unscaled);
538            } else {
539                // The following equality test is intentional and needed for rounding purposes
540                if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(FastMath.floor(unscaled) / 2.0)) { // even
541                    unscaled = FastMath.floor(unscaled);
542                } else { // odd
543                    unscaled = FastMath.ceil(unscaled);
544                }
545            }
546            break;
547        }
548        case BigDecimal.ROUND_HALF_UP : {
549            unscaled = FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY);
550            double fraction = unscaled - FastMath.floor(unscaled);
551            if (fraction >= 0.5) {
552                unscaled = FastMath.ceil(unscaled);
553            } else {
554                unscaled = FastMath.floor(unscaled);
555            }
556            break;
557        }
558        case BigDecimal.ROUND_UNNECESSARY :
559            if (unscaled != FastMath.floor(unscaled)) {
560                throw new MathArithmeticException();
561            }
562            break;
563        case BigDecimal.ROUND_UP :
564            // do not round if the discarded fraction is equal to zero
565            if (unscaled != FastMath.floor(unscaled)) {
566                unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
567            }
568            break;
569        default :
570            throw new MathIllegalArgumentException(LocalizedFormats.INVALID_ROUNDING_METHOD,
571                                                   roundingMethod,
572                                                   "ROUND_CEILING", BigDecimal.ROUND_CEILING,
573                                                   "ROUND_DOWN", BigDecimal.ROUND_DOWN,
574                                                   "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
575                                                   "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
576                                                   "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
577                                                   "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
578                                                   "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
579                                                   "ROUND_UP", BigDecimal.ROUND_UP);
580        }
581        return unscaled;
582    }
583
584
585    /**
586     * Computes a number {@code delta} close to {@code originalDelta} with
587     * the property that <pre><code>
588     *   x + delta - x
589     * </code></pre>
590     * is exactly machine-representable.
591     * This is useful when computing numerical derivatives, in order to reduce
592     * roundoff errors.
593     *
594     * @param x Value.
595     * @param originalDelta Offset value.
596     * @return a number {@code delta} so that {@code x + delta} and {@code x}
597     * differ by a representable floating number.
598     */
599    public static double representableDelta(double x,
600                                            double originalDelta) {
601        return x + originalDelta - x;
602    }
603}