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π:<br> 95 * {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li> 96 * <li>normalize an angle between -π and +π<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 π 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π interval for the result 106 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π 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π 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π:<br> 117 * {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li> 118 * <li>normalize an angle between -π and +π<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 π 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π interval for the result 129 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π 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 & 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 & 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 }