View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.util;
24  
25  import java.util.Arrays;
26  
27  import org.hipparchus.CalculusFieldElement;
28  import org.hipparchus.FieldElement;
29  import org.hipparchus.exception.Localizable;
30  import org.hipparchus.exception.LocalizedCoreFormats;
31  import org.hipparchus.exception.MathIllegalArgumentException;
32  import org.hipparchus.exception.MathRuntimeException;
33  import org.hipparchus.exception.NullArgumentException;
34  
35  /**
36   * Miscellaneous utility functions.
37   *
38   * @see ArithmeticUtils
39   * @see Precision
40   * @see MathArrays
41   */
42  public final class MathUtils {
43      /** \(2\pi\) */
44      public static final double TWO_PI = 2 * FastMath.PI;
45  
46      /** \(\pi^2\) */
47      public static final double PI_SQUARED = FastMath.PI * FastMath.PI;
48  
49      /** \(\pi/2\). */
50      public static final double SEMI_PI = 0.5 * FastMath.PI;
51  
52      /**
53       * Class contains only static methods.
54       */
55      private MathUtils() {}
56  
57  
58      /**
59       * Returns an integer hash code representing the given double value.
60       *
61       * @param value the value to be hashed
62       * @return the hash code
63       */
64      public static int hash(double value) {
65          return Double.hashCode(value);
66      }
67  
68      /**
69       * Returns {@code true} if the values are equal according to semantics of
70       * {@link Double#equals(Object)}.
71       *
72       * @param x Value
73       * @param y Value
74       * @return {@code Double.valueOf(x).equals(Double.valueOf(y))}
75       */
76      public static boolean equals(double x, double y) {
77          return Double.valueOf(x).equals(Double.valueOf(y));
78      }
79  
80      /**
81       * Returns an integer hash code representing the given double array.
82       *
83       * @param value the value to be hashed (may be null)
84       * @return the hash code
85       */
86      public static int hash(double[] value) {
87          return Arrays.hashCode(value);
88      }
89  
90      /**
91       * Normalize an angle in a 2π wide interval around a center value.
92       * <p>This method has three main uses:</p>
93       * <ul>
94       *   <li>normalize an angle between 0 and 2&pi;:<br>
95       *       {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
96       *   <li>normalize an angle between -&pi; and +&pi;<br>
97       *       {@code a = MathUtils.normalizeAngle(a, 0.0);}</li>
98       *   <li>compute the angle between two defining angular positions:<br>
99       *       {@code angle = MathUtils.normalizeAngle(end, start) - start;}</li>
100      * </ul>
101      * <p>Note that due to numerical accuracy and since &pi; cannot be represented
102      * exactly, the result interval is <em>closed</em>, it cannot be half-closed
103      * as would be more satisfactory in a purely mathematical view.</p>
104      * @param a angle to normalize
105      * @param center center of the desired 2&pi; interval for the result
106      * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
107      */
108      public static double normalizeAngle(double a, double center) {
109          return a - TWO_PI * FastMath.floor((a + FastMath.PI - center) / TWO_PI);
110      }
111 
112      /**
113       * Normalize an angle in a 2&pi; wide interval around a center value.
114       * <p>This method has three main uses:</p>
115       * <ul>
116       *   <li>normalize an angle between 0 and 2&pi;:<br>
117       *       {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
118       *   <li>normalize an angle between -&pi; and +&pi;<br>
119       *       {@code a = MathUtils.normalizeAngle(a, zero);}</li>
120       *   <li>compute the angle between two defining angular positions:<br>
121       *       {@code angle = MathUtils.normalizeAngle(end, start).subtract(start);}</li>
122       * </ul>
123       * <p>Note that due to numerical accuracy and since &pi; cannot be represented
124       * exactly, the result interval is <em>closed</em>, it cannot be half-closed
125       * as would be more satisfactory in a purely mathematical view.</p>
126       * @param <T> the type of the field elements
127       * @param a angle to normalize
128       * @param center center of the desired 2&pi; interval for the result
129       * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
130       */
131       public static <T extends CalculusFieldElement<T>> T normalizeAngle(T a, T center) {
132           return a.subtract(FastMath.floor(a.add(FastMath.PI).subtract(center).divide(TWO_PI)).multiply(TWO_PI));
133       }
134 
135      /** Find the maximum of two field elements.
136       * @param <T> the type of the field elements
137       * @param e1 first element
138       * @param e2 second element
139       * @return max(a1, e2)
140       */
141      public static <T extends CalculusFieldElement<T>> T max(final T e1, final T e2) {
142          return e1.subtract(e2).getReal() >= 0 ? e1 : e2;
143      }
144 
145      /** Find the minimum of two field elements.
146       * @param <T> the type of the field elements
147       * @param e1 first element
148       * @param e2 second element
149       * @return min(a1, e2)
150       */
151      public static <T extends CalculusFieldElement<T>> T min(final T e1, final T e2) {
152          return e1.subtract(e2).getReal() >= 0 ? e2 : e1;
153      }
154 
155     /**
156      * <p>Reduce {@code |a - offset|} to the primary interval
157      * {@code [0, |period|)}.</p>
158      *
159      * <p>Specifically, the value returned is <br>
160      * {@code a - |period| * floor((a - offset) / |period|) - offset}.</p>
161      *
162      * <p>If any of the parameters are {@code NaN} or infinite, the result is
163      * {@code NaN}.</p>
164      *
165      * @param a Value to reduce.
166      * @param period Period.
167      * @param offset Value that will be mapped to {@code 0}.
168      * @return the value, within the interval {@code [0 |period|)},
169      * that corresponds to {@code a}.
170      */
171     public static double reduce(double a,
172                                 double period,
173                                 double offset) {
174         final double p = FastMath.abs(period);
175         return a - p * FastMath.floor((a - offset) / p) - offset;
176     }
177 
178     /**
179      * Returns the first argument with the sign of the second argument.
180      *
181      * @param magnitude Magnitude of the returned value.
182      * @param sign Sign of the returned value.
183      * @return a value with magnitude equal to {@code magnitude} and with the
184      * same sign as the {@code sign} argument.
185      * @throws MathRuntimeException if {@code magnitude == Byte.MIN_VALUE}
186      * and {@code sign >= 0}.
187      */
188     public static byte copySign(byte magnitude, byte sign)
189         throws MathRuntimeException {
190         if ((magnitude >= 0 && sign >= 0) ||
191             (magnitude < 0 && sign < 0)) { // Sign is OK.
192             return magnitude;
193         } else if (sign >= 0 &&
194                    magnitude == Byte.MIN_VALUE) {
195             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
196         } else {
197             return (byte) -magnitude; // Flip sign.
198         }
199     }
200 
201     /**
202      * Returns the first argument with the sign of the second argument.
203      *
204      * @param magnitude Magnitude of the returned value.
205      * @param sign Sign of the returned value.
206      * @return a value with magnitude equal to {@code magnitude} and with the
207      * same sign as the {@code sign} argument.
208      * @throws MathRuntimeException if {@code magnitude == Short.MIN_VALUE}
209      * and {@code sign >= 0}.
210      */
211     public static short copySign(short magnitude, short sign)
212             throws MathRuntimeException {
213         if ((magnitude >= 0 && sign >= 0) ||
214             (magnitude < 0 && sign < 0)) { // Sign is OK.
215             return magnitude;
216         } else if (sign >= 0 &&
217                    magnitude == Short.MIN_VALUE) {
218             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
219         } else {
220             return (short) -magnitude; // Flip sign.
221         }
222     }
223 
224     /**
225      * Returns the first argument with the sign of the second argument.
226      *
227      * @param magnitude Magnitude of the returned value.
228      * @param sign Sign of the returned value.
229      * @return a value with magnitude equal to {@code magnitude} and with the
230      * same sign as the {@code sign} argument.
231      * @throws MathRuntimeException if {@code magnitude == Integer.MIN_VALUE}
232      * and {@code sign >= 0}.
233      */
234     public static int copySign(int magnitude, int sign)
235             throws MathRuntimeException {
236         if ((magnitude >= 0 && sign >= 0) ||
237             (magnitude < 0 && sign < 0)) { // Sign is OK.
238             return magnitude;
239         } else if (sign >= 0 &&
240                    magnitude == Integer.MIN_VALUE) {
241             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
242         } else {
243             return -magnitude; // Flip sign.
244         }
245     }
246 
247     /**
248      * Returns the first argument with the sign of the second argument.
249      *
250      * @param magnitude Magnitude of the returned value.
251      * @param sign Sign of the returned value.
252      * @return a value with magnitude equal to {@code magnitude} and with the
253      * same sign as the {@code sign} argument.
254      * @throws MathRuntimeException if {@code magnitude == Long.MIN_VALUE}
255      * and {@code sign >= 0}.
256      */
257     public static long copySign(long magnitude, long sign)
258         throws MathRuntimeException {
259         if ((magnitude >= 0 && sign >= 0) ||
260             (magnitude < 0 && sign < 0)) { // Sign is OK.
261             return magnitude;
262         } else if (sign >= 0 &&
263                    magnitude == Long.MIN_VALUE) {
264             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
265         } else {
266             return -magnitude; // Flip sign.
267         }
268     }
269     /**
270      * Check that the argument is a real number.
271      *
272      * @param x Argument.
273      * @throws MathIllegalArgumentException if {@code x} is not a
274      * finite real number.
275      */
276     public static void checkFinite(final double x)
277         throws MathIllegalArgumentException {
278         if (Double.isInfinite(x) || Double.isNaN(x)) {
279             throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_FINITE_NUMBER, x);
280         }
281     }
282 
283     /**
284      * Check that all the elements are real numbers.
285      *
286      * @param val Arguments.
287      * @throws MathIllegalArgumentException if any values of the array is not a
288      * finite real number.
289      */
290     public static void checkFinite(final double[] val)
291         throws MathIllegalArgumentException {
292         for (int i = 0; i < val.length; i++) {
293             final double x = val[i];
294             if (Double.isInfinite(x) || Double.isNaN(x)) {
295                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_FINITE_NUMBER, x);
296             }
297         }
298     }
299 
300     /**
301      * Checks that an object is not null.
302      *
303      * @param o Object to be checked.
304      * @param pattern Message pattern.
305      * @param args Arguments to replace the placeholders in {@code pattern}.
306      * @throws NullArgumentException if {@code o} is {@code null}.
307      */
308     public static void checkNotNull(Object o,
309                                     Localizable pattern,
310                                     Object ... args)
311         throws NullArgumentException {
312         if (o == null) {
313             throw new NullArgumentException(pattern, args);
314         }
315     }
316 
317     /**
318      * Checks that an object is not null.
319      *
320      * @param o Object to be checked.
321      * @throws NullArgumentException if {@code o} is {@code null}.
322      */
323     public static void checkNotNull(Object o)
324         throws NullArgumentException {
325         if (o == null) {
326             throw new NullArgumentException(LocalizedCoreFormats.NULL_NOT_ALLOWED);
327         }
328     }
329 
330     /**
331      * Checks that the given value is strictly within the range [lo, hi].
332      *
333      * @param value value to be checked.
334      * @param lo the lower bound (inclusive).
335      * @param hi the upper bound (inclusive).
336      * @throws MathIllegalArgumentException if {@code value} is strictly outside [lo, hi].
337      */
338     public static void checkRangeInclusive(long value, long lo, long hi) {
339         if (value < lo || value > hi) {
340             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
341                                                    value, lo, hi);
342         }
343     }
344 
345     /**
346      * Checks that the given value is strictly within the range [lo, hi].
347      *
348      * @param value value to be checked.
349      * @param lo the lower bound (inclusive).
350      * @param hi the upper bound (inclusive).
351      * @throws MathIllegalArgumentException if {@code value} is strictly outside [lo, hi].
352      */
353     public static void checkRangeInclusive(double value, double lo, double hi) {
354         if (value < lo || value > hi) {
355             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
356                                                    value, lo, hi);
357         }
358     }
359 
360     /**
361      * Checks that the given dimensions match.
362      *
363      * @param dimension the first dimension.
364      * @param otherDimension the second dimension.
365      * @throws MathIllegalArgumentException if length != otherLength.
366      */
367     public static void checkDimension(int dimension, int otherDimension) {
368         if (dimension != otherDimension) {
369             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
370                                                    dimension, otherDimension);
371         }
372     }
373 
374     /**
375      * Sums {@code a} and {@code b} using Møller's 2Sum algorithm.
376      * <p>
377      * References:
378      * <ul>
379      * <li>Møller, Ole. "Quasi double-precision in floating point addition." BIT
380      * 5, 37–50 (1965).</li>
381      * <li>Shewchuk, Richard J. "Adaptive Precision Floating-Point Arithmetic
382      * and Fast Robust Geometric Predicates." Discrete &amp; Computational Geometry
383      * 18, 305–363 (1997).</li>
384      * <li><a href=
385      * "https://en.wikipedia.org/wiki/2Sum">https://en.wikipedia.org/wiki/2Sum</a></li>
386      * </ul>
387      * @param a first summand
388      * @param b second summand
389      * @return sum and residual error in the sum
390      */
391     public static SumAndResidual twoSum(final double a, final double b) {
392         final double s = a + b;
393         final double aPrime = s - b;
394         final double bPrime = s - aPrime;
395         final double deltaA = a - aPrime;
396         final double deltaB = b - bPrime;
397         final double t = deltaA + deltaB;
398         return new SumAndResidual(s, t);
399     }
400 
401     /**
402      * Sums {@code a} and {@code b} using Møller's 2Sum algorithm.
403      * <p>
404      * References:
405      * <ul>
406      * <li>Møller, Ole. "Quasi double-precision in floating point addition." BIT
407      * 5, 37–50 (1965).</li>
408      * <li>Shewchuk, Richard J. "Adaptive Precision Floating-Point Arithmetic
409      * and Fast Robust Geometric Predicates." Discrete &amp; Computational Geometry
410      * 18, 305–363 (1997).</li>
411      * <li><a href=
412      * "https://en.wikipedia.org/wiki/2Sum">https://en.wikipedia.org/wiki/2Sum</a></li>
413      * </ul>
414      * @param <T> field element type
415      * @param a first summand
416      * @param b second summand
417      * @return sum and residual error in the sum
418      */
419     public static <T extends FieldElement<T>> FieldSumAndResidual<T> twoSum(final T a, final T b) {
420         final T s = a.add(b);
421         final T aPrime = s.subtract(b);
422         final T bPrime = s.subtract(aPrime);
423         final T deltaA = a.subtract(aPrime);
424         final T deltaB = b.subtract(bPrime);
425         final T t = deltaA.add(deltaB);
426         return new FieldSumAndResidual<>(s, t);
427     }
428 
429     /**
430      * Result class for {@link MathUtils#twoSum(double, double)} containing the
431      * sum and the residual error in the sum.
432      */
433     public static final class SumAndResidual {
434 
435         /** Sum. */
436         private final double sum;
437         /** Residual error in the sum. */
438         private final double residual;
439 
440         /**
441          * Constructs a {@link SumAndResidual} instance.
442          * @param sum sum
443          * @param residual residual error in the sum
444          */
445         private SumAndResidual(final double sum, final double residual) {
446             this.sum = sum;
447             this.residual = residual;
448         }
449 
450         /**
451          * Returns the sum.
452          * @return sum
453          */
454         public double getSum() {
455             return sum;
456         }
457 
458         /**
459          * Returns the residual error in the sum.
460          * @return residual error in the sum
461          */
462         public double getResidual() {
463             return residual;
464         }
465 
466     }
467 
468     /**
469      * Result class for
470      * {@link MathUtils#twoSum(FieldElement, FieldElement)} containing
471      * the sum and the residual error in the sum.
472      * @param <T> field element type
473      */
474     public static final class FieldSumAndResidual<T extends FieldElement<T>> {
475 
476         /** Sum. */
477         private final T sum;
478         /** Residual error in the sum. */
479         private final T residual;
480 
481         /**
482          * Constructs a {@link FieldSumAndResidual} instance.
483          * @param sum sum
484          * @param residual residual error in the sum
485          */
486         private FieldSumAndResidual(final T sum, final T residual) {
487             this.sum = sum;
488             this.residual = residual;
489         }
490 
491         /**
492          * Returns the sum.
493          * @return sum
494          */
495         public T getSum() {
496             return sum;
497         }
498 
499         /**
500          * Returns the residual error in the sum.
501          * @return residual error in the sum
502          */
503         public T getResidual() {
504             return residual;
505         }
506 
507     }
508 
509 }