Source code

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.numbers.core;
019
020import java.math.BigDecimal;
021import java.math.RoundingMode;
022
023/**
024 * Utilities for comparing numbers.
025 */
026public final class Precision {
027    /**
028     * <p>
029     * Largest double-precision floating-point number such that
030     * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
031     * bound on the relative error due to rounding real numbers to double
032     * precision floating-point numbers.
033     * </p>
034     * <p>
035     * In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
036     * </p>
037     *
038     * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
039     */
040    public static final double EPSILON;
041
042    /**
043     * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow.
044     * In IEEE 754 arithmetic, this is also the smallest normalized
045     * number 2<sup>-1022</sup>.
046     */
047    public static final double SAFE_MIN;
048
049    /** Exponent offset in IEEE754 representation. */
050    private static final long EXPONENT_OFFSET = 1023L;
051
052    /** Offset to order signed double numbers lexicographically. */
053    private static final long SGN_MASK = 0x8000000000000000L;
054    /** Offset to order signed double numbers lexicographically. */
055    private static final int SGN_MASK_FLOAT = 0x80000000;
056    /** Positive zero. */
057    private static final double POSITIVE_ZERO = 0d;
058    /** Positive zero bits. */
059    private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0);
060    /** Negative zero bits. */
061    private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0);
062    /** Positive zero bits. */
063    private static final int POSITIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(+0.0f);
064    /** Negative zero bits. */
065    private static final int NEGATIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(-0.0f);
066
067    static {
068        /*
069         *  This was previously expressed as = 0x1.0p-53
070         *  However, OpenJDK (Sparc Solaris) cannot handle such small
071         *  constants: MATH-721
072         */
073        EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
074
075        /*
076         * This was previously expressed as = 0x1.0p-1022
077         * However, OpenJDK (Sparc Solaris) cannot handle such small
078         * constants: MATH-721
079         */
080        SAFE_MIN = Double.longBitsToDouble((EXPONENT_OFFSET - 1022L) << 52);
081    }
082
083    /**
084     * Private constructor.
085     */
086    private Precision() {}
087
088    /**
089     * Compares two numbers given some amount of allowed error.
090     * The returned value is
091     * <ul>
092     *  <li>
093     *   0 if  {@link #equals(double,double,double) equals(x, y, eps)},
094     *  </li>
095     *  <li>
096     *   negative if !{@link #equals(double,double,double) equals(x, y, eps)} and {@code x < y},
097     *  </li>
098     *  <li>
099     *   positive if !{@link #equals(double,double,double) equals(x, y, eps)} and {@code x > y} or
100     *   either argument is {@code NaN}.
101     *  </li>
102     * </ul>
103     *
104     * @param x First value.
105     * @param y Second value.
106     * @param eps Allowed error when checking for equality.
107     * @return 0 if the value are considered equal, -1 if the first is smaller than
108     * the second, 1 is the first is larger than the second.
109     */
110    public static int compareTo(double x, double y, double eps) {
111        if (equals(x, y, eps)) {
112            return 0;
113        } else if (x < y) {
114            return -1;
115        }
116        return 1;
117    }
118
119    /**
120     * Compares two numbers given some amount of allowed error.
121     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
122     * (or fewer) floating point numbers between them, i.e. two adjacent floating
123     * point numbers are considered equal.
124     * Adapted from
125     * <a href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
126     * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN.
127     * The returned value is
128     * <ul>
129     *  <li>
130     *   zero if {@link #equals(double,double,int) equals(x, y, maxUlps)},
131     *  </li>
132     *  <li>
133     *   negative if !{@link #equals(double,double,int) equals(x, y, maxUlps)} and {@code x < y},
134     *  </li>
135     *  <li>
136     *   positive if !{@link #equals(double,double,int) equals(x, y, maxUlps)} and {@code x > y}
137     *       or either argument is {@code NaN}.
138     *  </li>
139     * </ul>
140     *
141     * @param x First value.
142     * @param y Second value.
143     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
144     * values between {@code x} and {@code y}.
145     * @return 0 if the value are considered equal, -1 if the first is smaller than
146     * the second, 1 is the first is larger than the second.
147     */
148    public static int compareTo(final double x, final double y, final int maxUlps) {
149        if (equals(x, y, maxUlps)) {
150            return 0;
151        } else if (x < y) {
152            return -1;
153        }
154        return 1;
155    }
156
157    /**
158     * Returns true iff they are equal as defined by
159     * {@link #equals(float,float,int) equals(x, y, 1)}.
160     *
161     * @param x first value
162     * @param y second value
163     * @return {@code true} if the values are equal.
164     */
165    public static boolean equals(float x, float y) {
166        return equals(x, y, 1);
167    }
168
169    /**
170     * Returns true if both arguments are NaN or they are
171     * equal as defined by {@link #equals(float,float) equals(x, y, 1)}.
172     *
173     * @param x first value
174     * @param y second value
175     * @return {@code true} if the values are equal or both are NaN.
176     */
177    public static boolean equalsIncludingNaN(float x, float y) {
178        final boolean xIsNan = Float.isNaN(x);
179        final boolean yIsNan = Float.isNaN(y);
180        return xIsNan || yIsNan ?
181            !(xIsNan ^ yIsNan) :
182            equals(x, y, 1);
183    }
184
185    /**
186     * Returns true if the arguments are equal or within the range of allowed
187     * error (inclusive).  Returns {@code false} if either of the arguments
188     * is NaN.
189     *
190     * @param x first value
191     * @param y second value
192     * @param eps the amount of absolute error to allow.
193     * @return {@code true} if the values are equal or within range of each other.
194     */
195    public static boolean equals(float x, float y, float eps) {
196        return equals(x, y, 1) || Math.abs(y - x) <= eps;
197    }
198
199    /**
200     * Returns true if the arguments are both NaN, are equal, or are within the range
201     * of allowed error (inclusive).
202     *
203     * @param x first value
204     * @param y second value
205     * @param eps the amount of absolute error to allow.
206     * @return {@code true} if the values are equal or within range of each other,
207     * or both are NaN.
208     */
209    public static boolean equalsIncludingNaN(float x, float y, float eps) {
210        return equalsIncludingNaN(x, y) || (Math.abs(y - x) <= eps);
211    }
212
213    /**
214     * Returns true if the arguments are equal or within the range of allowed
215     * error (inclusive).
216     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
217     * (or fewer) floating point numbers between them, i.e. two adjacent floating
218     * point numbers are considered equal.
219     * Adapted from <a
220     * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
221     * Bruce Dawson</a>.  Returns {@code false} if either of the arguments is NaN.
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 there are fewer than {@code maxUlps} floating
228     * point values between {@code x} and {@code y}.
229     */
230    public static boolean equals(final float x, final float y, final int maxUlps) {
231
232        final int xInt = Float.floatToRawIntBits(x);
233        final int yInt = Float.floatToRawIntBits(y);
234
235        final boolean isEqual;
236        if (((xInt ^ yInt) & SGN_MASK_FLOAT) == 0) {
237            // number have same sign, there is no risk of overflow
238            isEqual = Math.abs(xInt - yInt) <= maxUlps;
239        } else {
240            // number have opposite signs, take care of overflow
241            final int deltaPlus;
242            final int deltaMinus;
243            if (xInt < yInt) {
244                deltaPlus  = yInt - POSITIVE_ZERO_FLOAT_BITS;
245                deltaMinus = xInt - NEGATIVE_ZERO_FLOAT_BITS;
246            } else {
247                deltaPlus  = xInt - POSITIVE_ZERO_FLOAT_BITS;
248                deltaMinus = yInt - NEGATIVE_ZERO_FLOAT_BITS;
249            }
250
251            if (deltaPlus > maxUlps) {
252                isEqual = false;
253            } else {
254                isEqual = deltaMinus <= (maxUlps - deltaPlus);
255            }
256
257        }
258
259        return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
260
261    }
262
263    /**
264     * Returns true if the arguments are both NaN or if they are equal as defined
265     * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
266     *
267     * @param x first value
268     * @param y second value
269     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
270     * values between {@code x} and {@code y}.
271     * @return {@code true} if both arguments are NaN or if there are less than
272     * {@code maxUlps} floating point values between {@code x} and {@code y}.
273     */
274    public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
275        final boolean xIsNan = Float.isNaN(x);
276        final boolean yIsNan = Float.isNaN(y);
277        return xIsNan || yIsNan ?
278            !(xIsNan ^ yIsNan) :
279            equals(x, y, maxUlps);
280    }
281
282    /**
283     * Returns true iff they are equal as defined by
284     * {@link #equals(double,double,int) equals(x, y, 1)}.
285     *
286     * @param x first value
287     * @param y second value
288     * @return {@code true} if the values are equal.
289     */
290    public static boolean equals(double x, double y) {
291        return equals(x, y, 1);
292    }
293
294    /**
295     * Returns true if the arguments are both NaN or they are
296     * equal as defined by {@link #equals(double,double) equals(x, y, 1)}.
297     *
298     * @param x first value
299     * @param y second value
300     * @return {@code true} if the values are equal or both are NaN.
301     */
302    public static boolean equalsIncludingNaN(double x, double y) {
303        final boolean xIsNan = Double.isNaN(x);
304        final boolean yIsNan = Double.isNaN(y);
305        return xIsNan || yIsNan ?
306            !(xIsNan ^ yIsNan) :
307            equals(x, y, 1);
308    }
309
310    /**
311     * Returns {@code true} if there is no double value strictly between the
312     * arguments or the difference between them is within the range of allowed
313     * error (inclusive). Returns {@code false} if either of the arguments
314     * is NaN.
315     *
316     * @param x First value.
317     * @param y Second value.
318     * @param eps Amount of allowed absolute error.
319     * @return {@code true} if the values are two adjacent floating point
320     * numbers or they are within range of each other.
321     */
322    public static boolean equals(double x, double y, double eps) {
323        return equals(x, y, 1) || Math.abs(y - x) <= eps;
324    }
325
326    /**
327     * Returns {@code true} if there is no double value strictly between the
328     * arguments or the relative difference between them is less than or equal
329     * to the given tolerance. Returns {@code false} if either of the arguments
330     * is NaN.
331     *
332     * @param x First value.
333     * @param y Second value.
334     * @param eps Amount of allowed relative error.
335     * @return {@code true} if the values are two adjacent floating point
336     * numbers or they are within range of each other.
337     */
338    public static boolean equalsWithRelativeTolerance(double x, double y, double eps) {
339        if (equals(x, y, 1)) {
340            return true;
341        }
342
343        final double absoluteMax = Math.max(Math.abs(x), Math.abs(y));
344        final double relativeDifference = Math.abs((x - y) / absoluteMax);
345
346        return relativeDifference <= eps;
347    }
348
349    /**
350     * Returns true if the arguments are both NaN, are equal or are within the range
351     * of allowed error (inclusive).
352     *
353     * @param x first value
354     * @param y second value
355     * @param eps the amount of absolute error to allow.
356     * @return {@code true} if the values are equal or within range of each other,
357     * or both are NaN.
358     */
359    public static boolean equalsIncludingNaN(double x, double y, double eps) {
360        return equalsIncludingNaN(x, y) || (Math.abs(y - x) <= eps);
361    }
362
363    /**
364     * Returns true if the arguments are equal or within the range of allowed
365     * error (inclusive).
366     * <p>
367     * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
368     * (or fewer) floating point numbers between them, i.e. two adjacent
369     * floating point numbers are considered equal.
370     * </p>
371     * <p>
372     * Adapted from <a
373     * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
374     * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN.
375     * </p>
376     *
377     * @param x first value
378     * @param y second value
379     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
380     * values between {@code x} and {@code y}.
381     * @return {@code true} if there are fewer than {@code maxUlps} floating
382     * point values between {@code x} and {@code y}.
383     */
384    public static boolean equals(final double x, final double y, final int maxUlps) {
385
386        final long xInt = Double.doubleToRawLongBits(x);
387        final long yInt = Double.doubleToRawLongBits(y);
388
389        final boolean isEqual;
390        if (((xInt ^ yInt) & SGN_MASK) == 0L) {
391            // number have same sign, there is no risk of overflow
392            isEqual = Math.abs(xInt - yInt) <= maxUlps;
393        } else {
394            // number have opposite signs, take care of overflow
395            final long deltaPlus;
396            final long deltaMinus;
397            if (xInt < yInt) {
398                deltaPlus  = yInt - POSITIVE_ZERO_DOUBLE_BITS;
399                deltaMinus = xInt - NEGATIVE_ZERO_DOUBLE_BITS;
400            } else {
401                deltaPlus  = xInt - POSITIVE_ZERO_DOUBLE_BITS;
402                deltaMinus = yInt - NEGATIVE_ZERO_DOUBLE_BITS;
403            }
404
405            if (deltaPlus > maxUlps) {
406                isEqual = false;
407            } else {
408                isEqual = deltaMinus <= (maxUlps - deltaPlus);
409            }
410
411        }
412
413        return isEqual && !Double.isNaN(x) && !Double.isNaN(y);
414
415    }
416
417    /**
418     * Returns true if both arguments are NaN or if they are equal as defined
419     * by {@link #equals(double,double,int) equals(x, y, maxUlps)}.
420     *
421     * @param x first value
422     * @param y second value
423     * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
424     * values between {@code x} and {@code y}.
425     * @return {@code true} if both arguments are NaN or if there are less than
426     * {@code maxUlps} floating point values between {@code x} and {@code y}.
427     */
428    public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
429        final boolean xIsNan = Double.isNaN(x);
430        final boolean yIsNan = Double.isNaN(y);
431        return xIsNan || yIsNan ?
432            !(xIsNan ^ yIsNan) :
433            equals(x, y, maxUlps);
434    }
435
436    /**
437     * Rounds the given value to the specified number of decimal places.
438     * The value is rounded using the {@link RoundingMode#HALF_UP} method.
439     *
440     * @param x Value to round.
441     * @param scale Number of digits to the right of the decimal point.
442     * @return the rounded value.
443     */
444    public static double round(double x, int scale) {
445        return round(x, scale, RoundingMode.HALF_UP);
446    }
447
448    /**
449     * Rounds the given value to the specified number of decimal places.
450     * The value is rounded using the given method which is any method defined
451     * in {@link BigDecimal}.
452     * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is
453     * returned unchanged, regardless of the other parameters.
454     *
455     * @param x Value to round.
456     * @param scale Number of digits to the right of the decimal point.
457     * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
458     * @return the rounded value.
459     * @throws ArithmeticException if {@code roundingMethod} is
460     * {@link RoundingMode#UNNECESSARY} and the specified scaling operation
461     * would require rounding.
462     */
463    public static double round(double x,
464                               int scale,
465                               RoundingMode roundingMethod) {
466        try {
467            final double rounded = (new BigDecimal(Double.toString(x))
468                   .setScale(scale, roundingMethod))
469                   .doubleValue();
470            // MATH-1089: negative values rounded to zero should result in negative zero
471            return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded;
472        } catch (NumberFormatException ex) {
473            if (Double.isInfinite(x)) {
474                return x;
475            }
476            return Double.NaN;
477        }
478    }
479
480    /**
481     * Computes a number {@code delta} close to {@code originalDelta} with
482     * the property that <pre><code>
483     *   x + delta - x
484     * </code></pre>
485     * is exactly machine-representable.
486     * This is useful when computing numerical derivatives, in order to reduce
487     * roundoff errors.
488     *
489     * @param x Value.
490     * @param originalDelta Offset value.
491     * @return a number {@code delta} so that {@code x + delta} and {@code x}
492     * differ by a representable floating number.
493     */
494    public static double representableDelta(double x,
495                                            double originalDelta) {
496        return x + originalDelta - x;
497    }
498}