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
018 package org.apache.commons.math.util;
019
020 import java.math.BigDecimal;
021 import java.math.BigInteger;
022 import java.util.Arrays;
023
024 import org.apache.commons.math.exception.MathArithmeticException;
025 import org.apache.commons.math.exception.MathIllegalArgumentException;
026 import org.apache.commons.math.exception.NotFiniteNumberException;
027 import org.apache.commons.math.exception.NotPositiveException;
028 import org.apache.commons.math.exception.NullArgumentException;
029 import org.apache.commons.math.exception.util.Localizable;
030 import org.apache.commons.math.exception.util.LocalizedFormats;
031
032 /**
033 * Some useful additions to the built-in functions in {@link Math}.
034 * @version $Id: MathUtils.java 1185841 2011-10-18 20:30:42Z erans $
035 */
036 public final class MathUtils {
037
038 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
039 public static final double EPSILON = 0x1.0p-53;
040
041 /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
042 * <p>In IEEE 754 arithmetic, this is also the smallest normalized
043 * number 2<sup>-1022</sup>.</p>
044 */
045 public static final double SAFE_MIN = 0x1.0p-1022;
046
047 /**
048 * 2 π.
049 * @since 2.1
050 */
051 public static final double TWO_PI = 2 * FastMath.PI;
052
053 /** -1.0 cast as a byte. */
054 private static final byte NB = (byte)-1;
055
056 /** -1.0 cast as a short. */
057 private static final short NS = (short)-1;
058
059 /** 1.0 cast as a byte. */
060 private static final byte PB = (byte)1;
061
062 /** 1.0 cast as a short. */
063 private static final short PS = (short)1;
064
065 /** 0.0 cast as a byte. */
066 private static final byte ZB = (byte)0;
067
068 /** 0.0 cast as a short. */
069 private static final short ZS = (short)0;
070
071 /**
072 * Private Constructor
073 */
074 private MathUtils() {
075 super();
076 }
077
078 /**
079 * Returns an integer hash code representing the given double value.
080 *
081 * @param value the value to be hashed
082 * @return the hash code
083 */
084 public static int hash(double value) {
085 return new Double(value).hashCode();
086 }
087
088 /**
089 * Returns an integer hash code representing the given double array.
090 *
091 * @param value the value to be hashed (may be null)
092 * @return the hash code
093 * @since 1.2
094 */
095 public static int hash(double[] value) {
096 return Arrays.hashCode(value);
097 }
098
099 /**
100 * For a byte value x, this method returns (byte)(+1) if x >= 0 and
101 * (byte)(-1) if x < 0.
102 *
103 * @param x the value, a byte
104 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
105 */
106 public static byte indicator(final byte x) {
107 return (x >= ZB) ? PB : NB;
108 }
109
110 /**
111 * For a double precision value x, this method returns +1.0 if x >= 0 and
112 * -1.0 if x < 0. Returns {@code NaN} if {@code x} is
113 * {@code NaN}.
114 *
115 * @param x the value, a double
116 * @return +1.0 or -1.0, depending on the sign of x
117 */
118 public static double indicator(final double x) {
119 if (Double.isNaN(x)) {
120 return Double.NaN;
121 }
122 return (x >= 0.0) ? 1.0 : -1.0;
123 }
124
125 /**
126 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
127 * 0. Returns {@code NaN} if {@code x} is {@code NaN}.
128 *
129 * @param x the value, a float
130 * @return +1.0F or -1.0F, depending on the sign of x
131 */
132 public static float indicator(final float x) {
133 if (Float.isNaN(x)) {
134 return Float.NaN;
135 }
136 return (x >= 0.0F) ? 1.0F : -1.0F;
137 }
138
139 /**
140 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
141 *
142 * @param x the value, an int
143 * @return +1 or -1, depending on the sign of x
144 */
145 public static int indicator(final int x) {
146 return (x >= 0) ? 1 : -1;
147 }
148
149 /**
150 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
151 *
152 * @param x the value, a long
153 * @return +1L or -1L, depending on the sign of x
154 */
155 public static long indicator(final long x) {
156 return (x >= 0L) ? 1L : -1L;
157 }
158
159 /**
160 * For a short value x, this method returns (short)(+1) if x >= 0 and
161 * (short)(-1) if x < 0.
162 *
163 * @param x the value, a short
164 * @return (short)(+1) or (short)(-1), depending on the sign of x
165 */
166 public static short indicator(final short x) {
167 return (x >= ZS) ? PS : NS;
168 }
169
170 /**
171 * <p>Returns the
172 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
173 * for base {@code b} of {@code x}.
174 * </p>
175 * <p>Returns {@code NaN} if either argument is negative. If
176 * {@code base} is 0 and {@code x} is positive, 0 is returned.
177 * If {@code base} is positive and {@code x} is 0,
178 * {@code Double.NEGATIVE_INFINITY} is returned. If both arguments
179 * are 0, the result is {@code NaN}.</p>
180 *
181 * @param base the base of the logarithm, must be greater than 0
182 * @param x argument, must be greater than 0
183 * @return the value of the logarithm - the number y such that base^y = x.
184 * @since 1.2
185 */
186 public static double log(double base, double x) {
187 return FastMath.log(x)/FastMath.log(base);
188 }
189
190 /**
191 * Normalize an angle in a 2&pi wide interval around a center value.
192 * <p>This method has three main uses:</p>
193 * <ul>
194 * <li>normalize an angle between 0 and 2π:<br/>
195 * {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
196 * <li>normalize an angle between -π and +π<br/>
197 * {@code a = MathUtils.normalizeAngle(a, 0.0);}</li>
198 * <li>compute the angle between two defining angular positions:<br>
199 * {@code angle = MathUtils.normalizeAngle(end, start) - start;}</li>
200 * </ul>
201 * <p>Note that due to numerical accuracy and since π cannot be represented
202 * exactly, the result interval is <em>closed</em>, it cannot be half-closed
203 * as would be more satisfactory in a purely mathematical view.</p>
204 * @param a angle to normalize
205 * @param center center of the desired 2π interval for the result
206 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π
207 * @since 1.2
208 */
209 public static double normalizeAngle(double a, double center) {
210 return a - TWO_PI * FastMath.floor((a + FastMath.PI - center) / TWO_PI);
211 }
212
213 /**
214 * <p>Reduce {@code |a - offset|} to the primary interval
215 * {@code [0, |period|)}.</p>
216 *
217 * <p>Specifically, the value returned is <br/>
218 * {@code a - |period| * floor((a - offset) / |period|) - offset}.</p>
219 *
220 * <p>If any of the parameters are {@code NaN} or infinite, the result is
221 * {@code NaN}.</p>
222 *
223 * @param a Value to reduce.
224 * @param period Period.
225 * @param offset Value that will be mapped to {@code 0}.
226 * @return the value, within the interval {@code [0 |period|)},
227 * that corresponds to {@code a}.
228 */
229 public static double reduce(double a,
230 double period,
231 double offset) {
232 final double p = FastMath.abs(period);
233 return a - p * FastMath.floor((a - offset) / p) - offset;
234 }
235
236 /**
237 * Round the given value to the specified number of decimal places. The
238 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
239 *
240 * @param x the value to round.
241 * @param scale the number of digits to the right of the decimal point.
242 * @return the rounded value.
243 * @since 1.1
244 */
245 public static double round(double x, int scale) {
246 return round(x, scale, BigDecimal.ROUND_HALF_UP);
247 }
248
249 /**
250 * <p>Round the given value to the specified number of decimal places. The
251 * value is rounded using the given method which is any method defined in
252 * {@link BigDecimal}.</p>
253 *
254 * <p>If {@code x} is infinite or NaN, then the value of {@code x} is
255 * returned unchanged, regardless of the other parameters.</p>
256 *
257 * @param x the value to round.
258 * @param scale the number of digits to the right of the decimal point.
259 * @param roundingMethod the rounding method as defined in
260 * {@link BigDecimal}.
261 * @return the rounded value.
262 * @throws ArithmeticException if roundingMethod==ROUND_UNNECESSARY and the
263 * specified scaling operation would require rounding.
264 * @throws IllegalArgumentException if roundingMethod does not represent a
265 * valid rounding mode.
266 * @since 1.1
267 */
268 public static double round(double x, int scale, int roundingMethod) {
269 try {
270 return (new BigDecimal
271 (Double.toString(x))
272 .setScale(scale, roundingMethod))
273 .doubleValue();
274 } catch (NumberFormatException ex) {
275 if (Double.isInfinite(x)) {
276 return x;
277 } else {
278 return Double.NaN;
279 }
280 }
281 }
282
283 /**
284 * Round the given value to the specified number of decimal places. The
285 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
286 *
287 * @param x the value to round.
288 * @param scale the number of digits to the right of the decimal point.
289 * @return the rounded value.
290 * @since 1.1
291 */
292 public static float round(float x, int scale) {
293 return round(x, scale, BigDecimal.ROUND_HALF_UP);
294 }
295
296 /**
297 * Round the given value to the specified number of decimal places. The
298 * value is rounded using the given method which is any method defined in
299 * {@link BigDecimal}.
300 *
301 * @param x the value to round.
302 * @param scale the number of digits to the right of the decimal point.
303 * @param roundingMethod the rounding method as defined in
304 * {@link BigDecimal}.
305 * @return the rounded value.
306 * @since 1.1
307 */
308 public static float round(float x, int scale, int roundingMethod) {
309 float sign = indicator(x);
310 float factor = (float)FastMath.pow(10.0f, scale) * sign;
311 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
312 }
313
314 /**
315 * Round the given non-negative value to the "nearest" integer. Nearest is
316 * determined by the rounding method specified. Rounding methods are defined
317 * in {@link BigDecimal}.
318 *
319 * @param unscaled Value to round.
320 * @param sign Sign of the original, scaled value.
321 * @param roundingMethod Rounding method, as defined in {@link BigDecimal}.
322 * @return the rounded value.
323 * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method.
324 * @since 1.1
325 */
326 private static double roundUnscaled(double unscaled,
327 double sign,
328 int roundingMethod) {
329 switch (roundingMethod) {
330 case BigDecimal.ROUND_CEILING :
331 if (sign == -1) {
332 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
333 } else {
334 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
335 }
336 break;
337 case BigDecimal.ROUND_DOWN :
338 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
339 break;
340 case BigDecimal.ROUND_FLOOR :
341 if (sign == -1) {
342 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
343 } else {
344 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY));
345 }
346 break;
347 case BigDecimal.ROUND_HALF_DOWN : {
348 unscaled = FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY);
349 double fraction = unscaled - FastMath.floor(unscaled);
350 if (fraction > 0.5) {
351 unscaled = FastMath.ceil(unscaled);
352 } else {
353 unscaled = FastMath.floor(unscaled);
354 }
355 break;
356 }
357 case BigDecimal.ROUND_HALF_EVEN : {
358 double fraction = unscaled - FastMath.floor(unscaled);
359 if (fraction > 0.5) {
360 unscaled = FastMath.ceil(unscaled);
361 } else if (fraction < 0.5) {
362 unscaled = FastMath.floor(unscaled);
363 } else {
364 // The following equality test is intentional and needed for rounding purposes
365 if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(Math
366 .floor(unscaled) / 2.0)) { // even
367 unscaled = FastMath.floor(unscaled);
368 } else { // odd
369 unscaled = FastMath.ceil(unscaled);
370 }
371 }
372 break;
373 }
374 case BigDecimal.ROUND_HALF_UP : {
375 unscaled = FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY);
376 double fraction = unscaled - FastMath.floor(unscaled);
377 if (fraction >= 0.5) {
378 unscaled = FastMath.ceil(unscaled);
379 } else {
380 unscaled = FastMath.floor(unscaled);
381 }
382 break;
383 }
384 case BigDecimal.ROUND_UNNECESSARY :
385 if (unscaled != FastMath.floor(unscaled)) {
386 throw new MathArithmeticException();
387 }
388 break;
389 case BigDecimal.ROUND_UP :
390 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY));
391 break;
392 default :
393 throw new MathIllegalArgumentException(LocalizedFormats.INVALID_ROUNDING_METHOD,
394 roundingMethod,
395 "ROUND_CEILING", BigDecimal.ROUND_CEILING,
396 "ROUND_DOWN", BigDecimal.ROUND_DOWN,
397 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
398 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
399 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
400 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
401 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
402 "ROUND_UP", BigDecimal.ROUND_UP);
403 }
404 return unscaled;
405 }
406
407 /**
408 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
409 * for byte value {@code x}.
410 * <p>
411 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
412 * x = 0, and (byte)(-1) if x < 0.</p>
413 *
414 * @param x the value, a byte
415 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
416 */
417 public static byte sign(final byte x) {
418 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
419 }
420
421 /**
422 * Returns the first argument with the sign of the second argument.
423 *
424 * @param magnitude Magnitude of the returned value.
425 * @param sign Sign of the returned value.
426 * @return a value with magnitude equal to {@code magnitude} and with the
427 * same sign as the {@code sign} argument.
428 * @throws MathArithmeticException if {@code magnitude == Byte.MIN_VALUE}
429 * and {@code sign >= 0}.
430 */
431 public static byte copySign(byte magnitude, byte sign) {
432 if ((magnitude >= 0 && sign >= 0) ||
433 (magnitude < 0 && sign < 0)) { // Sign is OK.
434 return magnitude;
435 } else if (sign >= 0 &&
436 magnitude == Byte.MIN_VALUE) {
437 throw new MathArithmeticException(LocalizedFormats.OVERFLOW);
438 } else {
439 return (byte) -magnitude; // Flip sign.
440 }
441 }
442
443 /**
444 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
445 * for int value {@code x}.
446 * <p>
447 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
448 * if x < 0.</p>
449 *
450 * @param x the value, an int
451 * @return +1, 0, or -1, depending on the sign of x
452 */
453 public static int sign(final int x) {
454 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
455 }
456
457 /**
458 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
459 * for long value {@code x}.
460 * <p>
461 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
462 * -1L if x < 0.</p>
463 *
464 * @param x the value, a long
465 * @return +1L, 0L, or -1L, depending on the sign of x
466 */
467 public static long sign(final long x) {
468 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
469 }
470
471 /**
472 * Compute the <a href="http://mathworld.wolfram.com/Sign.html">sign</a>
473 * of the argument.
474 *
475 * @param x the value, a short
476 * @return 1 if {@code x > 0}, 0 if {@code x == 0}, and -1 if {@code x < 0}.
477 */
478 public static short sign(final short x) {
479 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
480 }
481
482 /**
483 * Raise an int to an int power.
484 *
485 * @param k Number to raise.
486 * @param e Exponent (must be positive or zero).
487 * @return k<sup>e</sup>
488 * @throws NotPositiveException if {@code e < 0}.
489 */
490 public static int pow(final int k, int e) {
491 if (e < 0) {
492 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
493 }
494
495 int result = 1;
496 int k2p = k;
497 while (e != 0) {
498 if ((e & 0x1) != 0) {
499 result *= k2p;
500 }
501 k2p *= k2p;
502 e = e >> 1;
503 }
504
505 return result;
506 }
507
508 /**
509 * Raise an int to a long power.
510 *
511 * @param k Number to raise.
512 * @param e Exponent (must be positive or zero).
513 * @return k<sup>e</sup>
514 * @throws NotPositiveException if {@code e < 0}.
515 */
516 public static int pow(final int k, long e) {
517 if (e < 0) {
518 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
519 }
520
521 int result = 1;
522 int k2p = k;
523 while (e != 0) {
524 if ((e & 0x1) != 0) {
525 result *= k2p;
526 }
527 k2p *= k2p;
528 e = e >> 1;
529 }
530
531 return result;
532 }
533
534 /**
535 * Raise a long to an int power.
536 *
537 * @param k Number to raise.
538 * @param e Exponent (must be positive or zero).
539 * @return k<sup>e</sup>
540 * @throws NotPositiveException if {@code e < 0}.
541 */
542 public static long pow(final long k, int e) {
543 if (e < 0) {
544 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
545 }
546
547 long result = 1l;
548 long k2p = k;
549 while (e != 0) {
550 if ((e & 0x1) != 0) {
551 result *= k2p;
552 }
553 k2p *= k2p;
554 e = e >> 1;
555 }
556
557 return result;
558 }
559
560 /**
561 * Raise a long to a long power.
562 *
563 * @param k Number to raise.
564 * @param e Exponent (must be positive or zero).
565 * @return k<sup>e</sup>
566 * @throws NotPositiveException if {@code e < 0}.
567 */
568 public static long pow(final long k, long e) {
569 if (e < 0) {
570 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
571 }
572
573 long result = 1l;
574 long k2p = k;
575 while (e != 0) {
576 if ((e & 0x1) != 0) {
577 result *= k2p;
578 }
579 k2p *= k2p;
580 e = e >> 1;
581 }
582
583 return result;
584 }
585
586 /**
587 * Raise a BigInteger to an int power.
588 *
589 * @param k Number to raise.
590 * @param e Exponent (must be positive or zero).
591 * @return k<sup>e</sup>
592 * @throws NotPositiveException if {@code e < 0}.
593 */
594 public static BigInteger pow(final BigInteger k, int e) {
595 if (e < 0) {
596 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
597 }
598
599 return k.pow(e);
600 }
601
602 /**
603 * Raise a BigInteger to a long power.
604 *
605 * @param k Number to raise.
606 * @param e Exponent (must be positive or zero).
607 * @return k<sup>e</sup>
608 * @throws NotPositiveException if {@code e < 0}.
609 */
610 public static BigInteger pow(final BigInteger k, long e) {
611 if (e < 0) {
612 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
613 }
614
615 BigInteger result = BigInteger.ONE;
616 BigInteger k2p = k;
617 while (e != 0) {
618 if ((e & 0x1) != 0) {
619 result = result.multiply(k2p);
620 }
621 k2p = k2p.multiply(k2p);
622 e = e >> 1;
623 }
624
625 return result;
626
627 }
628
629 /**
630 * Raise a BigInteger to a BigInteger power.
631 *
632 * @param k Number to raise.
633 * @param e Exponent (must be positive or zero).
634 * @return k<sup>e</sup>
635 * @throws NotPositiveException if {@code e < 0}.
636 */
637 public static BigInteger pow(final BigInteger k, BigInteger e) {
638 if (e.compareTo(BigInteger.ZERO) < 0) {
639 throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
640 }
641
642 BigInteger result = BigInteger.ONE;
643 BigInteger k2p = k;
644 while (!BigInteger.ZERO.equals(e)) {
645 if (e.testBit(0)) {
646 result = result.multiply(k2p);
647 }
648 k2p = k2p.multiply(k2p);
649 e = e.shiftRight(1);
650 }
651
652 return result;
653 }
654
655 /**
656 * Check that the argument is a real number.
657 *
658 * @param x Argument.
659 * @throws NotFiniteNumberException if {@code x} is not a
660 * finite real number.
661 */
662 public static void checkFinite(final double x) {
663 if (Double.isInfinite(x) || Double.isNaN(x)) {
664 throw new NotFiniteNumberException(x);
665 }
666 }
667
668 /**
669 * Check that all the elements are real number.
670 *
671 * @param val Arguments.
672 * @throws NotFiniteNumberException if any values of the array is not a
673 * finite real number.
674 */
675 public static void checkFinite(final double[] val) {
676 for (int i = 0; i < val.length; i++) {
677 final double x = val[i];
678 if (Double.isInfinite(x) || Double.isNaN(x)) {
679 throw new NotFiniteNumberException(LocalizedFormats.ARRAY_ELEMENT, x, i);
680 }
681 }
682 }
683
684 /**
685 * Checks that an object is not null.
686 *
687 * @param o Object to be checked.
688 * @param pattern Message pattern.
689 * @param args Arguments to replace the placeholders in {@code pattern}.
690 * @throws NullArgumentException if {@code o} is {@code null}.
691 */
692 public static void checkNotNull(Object o,
693 Localizable pattern,
694 Object ... args) {
695 if (o == null) {
696 throw new NullArgumentException(pattern, args);
697 }
698 }
699
700 /**
701 * Checks that an object is not null.
702 *
703 * @param o Object to be checked.
704 * @throws NullArgumentException if {@code o} is {@code null}.
705 */
706 public static void checkNotNull(Object o)
707 throws NullArgumentException {
708 if (o == null) {
709 throw new NullArgumentException();
710 }
711 }
712 }