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