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.MathRuntimeException;
025
026 /**
027 * Some useful additions to the built-in functions in {@link Math}.
028 * @version $Revision: 927249 $ $Date: 2010-03-24 21:06:51 -0400 (Wed, 24 Mar 2010) $
029 */
030 public final class MathUtils {
031
032 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
033 public static final double EPSILON = 0x1.0p-53;
034
035 /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
036 * <p>In IEEE 754 arithmetic, this is also the smallest normalized
037 * number 2<sup>-1022</sup>.</p>
038 */
039 public static final double SAFE_MIN = 0x1.0p-1022;
040
041 /**
042 * 2 π.
043 * @since 2.1
044 */
045 public static final double TWO_PI = 2 * Math.PI;
046
047 /** -1.0 cast as a byte. */
048 private static final byte NB = (byte)-1;
049
050 /** -1.0 cast as a short. */
051 private static final short NS = (short)-1;
052
053 /** 1.0 cast as a byte. */
054 private static final byte PB = (byte)1;
055
056 /** 1.0 cast as a short. */
057 private static final short PS = (short)1;
058
059 /** 0.0 cast as a byte. */
060 private static final byte ZB = (byte)0;
061
062 /** 0.0 cast as a short. */
063 private static final short ZS = (short)0;
064
065 /** Gap between NaN and regular numbers. */
066 private static final int NAN_GAP = 4 * 1024 * 1024;
067
068 /** Offset to order signed double numbers lexicographically. */
069 private static final long SGN_MASK = 0x8000000000000000L;
070
071 /** All long-representable factorials */
072 private static final long[] FACTORIALS = new long[] {
073 1l, 1l, 2l,
074 6l, 24l, 120l,
075 720l, 5040l, 40320l,
076 362880l, 3628800l, 39916800l,
077 479001600l, 6227020800l, 87178291200l,
078 1307674368000l, 20922789888000l, 355687428096000l,
079 6402373705728000l, 121645100408832000l, 2432902008176640000l };
080
081 /**
082 * Private Constructor
083 */
084 private MathUtils() {
085 super();
086 }
087
088 /**
089 * Add two integers, checking for overflow.
090 *
091 * @param x an addend
092 * @param y an addend
093 * @return the sum <code>x+y</code>
094 * @throws ArithmeticException if the result can not be represented as an
095 * int
096 * @since 1.1
097 */
098 public static int addAndCheck(int x, int y) {
099 long s = (long)x + (long)y;
100 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
101 throw new ArithmeticException("overflow: add");
102 }
103 return (int)s;
104 }
105
106 /**
107 * Add two long integers, checking for overflow.
108 *
109 * @param a an addend
110 * @param b an addend
111 * @return the sum <code>a+b</code>
112 * @throws ArithmeticException if the result can not be represented as an
113 * long
114 * @since 1.2
115 */
116 public static long addAndCheck(long a, long b) {
117 return addAndCheck(a, b, "overflow: add");
118 }
119
120 /**
121 * Add two long integers, checking for overflow.
122 *
123 * @param a an addend
124 * @param b an addend
125 * @param msg the message to use for any thrown exception.
126 * @return the sum <code>a+b</code>
127 * @throws ArithmeticException if the result can not be represented as an
128 * long
129 * @since 1.2
130 */
131 private static long addAndCheck(long a, long b, String msg) {
132 long ret;
133 if (a > b) {
134 // use symmetry to reduce boundary cases
135 ret = addAndCheck(b, a, msg);
136 } else {
137 // assert a <= b
138
139 if (a < 0) {
140 if (b < 0) {
141 // check for negative overflow
142 if (Long.MIN_VALUE - b <= a) {
143 ret = a + b;
144 } else {
145 throw new ArithmeticException(msg);
146 }
147 } else {
148 // opposite sign addition is always safe
149 ret = a + b;
150 }
151 } else {
152 // assert a >= 0
153 // assert b >= 0
154
155 // check for positive overflow
156 if (a <= Long.MAX_VALUE - b) {
157 ret = a + b;
158 } else {
159 throw new ArithmeticException(msg);
160 }
161 }
162 }
163 return ret;
164 }
165
166 /**
167 * Returns an exact representation of the <a
168 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
169 * Coefficient</a>, "<code>n choose k</code>", the number of
170 * <code>k</code>-element subsets that can be selected from an
171 * <code>n</code>-element set.
172 * <p>
173 * <Strong>Preconditions</strong>:
174 * <ul>
175 * <li> <code>0 <= k <= n </code> (otherwise
176 * <code>IllegalArgumentException</code> is thrown)</li>
177 * <li> The result is small enough to fit into a <code>long</code>. The
178 * largest value of <code>n</code> for which all coefficients are
179 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
180 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is
181 * thrown.</li>
182 * </ul></p>
183 *
184 * @param n the size of the set
185 * @param k the size of the subsets to be counted
186 * @return <code>n choose k</code>
187 * @throws IllegalArgumentException if preconditions are not met.
188 * @throws ArithmeticException if the result is too large to be represented
189 * by a long integer.
190 */
191 public static long binomialCoefficient(final int n, final int k) {
192 checkBinomial(n, k);
193 if ((n == k) || (k == 0)) {
194 return 1;
195 }
196 if ((k == 1) || (k == n - 1)) {
197 return n;
198 }
199 // Use symmetry for large k
200 if (k > n / 2)
201 return binomialCoefficient(n, n - k);
202
203 // We use the formula
204 // (n choose k) = n! / (n-k)! / k!
205 // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
206 // which could be written
207 // (n choose k) == (n-1 choose k-1) * n / k
208 long result = 1;
209 if (n <= 61) {
210 // For n <= 61, the naive implementation cannot overflow.
211 int i = n - k + 1;
212 for (int j = 1; j <= k; j++) {
213 result = result * i / j;
214 i++;
215 }
216 } else if (n <= 66) {
217 // For n > 61 but n <= 66, the result cannot overflow,
218 // but we must take care not to overflow intermediate values.
219 int i = n - k + 1;
220 for (int j = 1; j <= k; j++) {
221 // We know that (result * i) is divisible by j,
222 // but (result * i) may overflow, so we split j:
223 // Filter out the gcd, d, so j/d and i/d are integer.
224 // result is divisible by (j/d) because (j/d)
225 // is relative prime to (i/d) and is a divisor of
226 // result * (i/d).
227 final long d = gcd(i, j);
228 result = (result / (j / d)) * (i / d);
229 i++;
230 }
231 } else {
232 // For n > 66, a result overflow might occur, so we check
233 // the multiplication, taking care to not overflow
234 // unnecessary.
235 int i = n - k + 1;
236 for (int j = 1; j <= k; j++) {
237 final long d = gcd(i, j);
238 result = mulAndCheck(result / (j / d), i / d);
239 i++;
240 }
241 }
242 return result;
243 }
244
245 /**
246 * Returns a <code>double</code> representation of the <a
247 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
248 * Coefficient</a>, "<code>n choose k</code>", the number of
249 * <code>k</code>-element subsets that can be selected from an
250 * <code>n</code>-element set.
251 * <p>
252 * <Strong>Preconditions</strong>:
253 * <ul>
254 * <li> <code>0 <= k <= n </code> (otherwise
255 * <code>IllegalArgumentException</code> is thrown)</li>
256 * <li> The result is small enough to fit into a <code>double</code>. The
257 * largest value of <code>n</code> for which all coefficients are <
258 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
259 * Double.POSITIVE_INFINITY is returned</li>
260 * </ul></p>
261 *
262 * @param n the size of the set
263 * @param k the size of the subsets to be counted
264 * @return <code>n choose k</code>
265 * @throws IllegalArgumentException if preconditions are not met.
266 */
267 public static double binomialCoefficientDouble(final int n, final int k) {
268 checkBinomial(n, k);
269 if ((n == k) || (k == 0)) {
270 return 1d;
271 }
272 if ((k == 1) || (k == n - 1)) {
273 return n;
274 }
275 if (k > n/2) {
276 return binomialCoefficientDouble(n, n - k);
277 }
278 if (n < 67) {
279 return binomialCoefficient(n,k);
280 }
281
282 double result = 1d;
283 for (int i = 1; i <= k; i++) {
284 result *= (double)(n - k + i) / (double)i;
285 }
286
287 return Math.floor(result + 0.5);
288 }
289
290 /**
291 * Returns the natural <code>log</code> of the <a
292 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
293 * Coefficient</a>, "<code>n choose k</code>", the number of
294 * <code>k</code>-element subsets that can be selected from an
295 * <code>n</code>-element set.
296 * <p>
297 * <Strong>Preconditions</strong>:
298 * <ul>
299 * <li> <code>0 <= k <= n </code> (otherwise
300 * <code>IllegalArgumentException</code> is thrown)</li>
301 * </ul></p>
302 *
303 * @param n the size of the set
304 * @param k the size of the subsets to be counted
305 * @return <code>n choose k</code>
306 * @throws IllegalArgumentException if preconditions are not met.
307 */
308 public static double binomialCoefficientLog(final int n, final int k) {
309 checkBinomial(n, k);
310 if ((n == k) || (k == 0)) {
311 return 0;
312 }
313 if ((k == 1) || (k == n - 1)) {
314 return Math.log(n);
315 }
316
317 /*
318 * For values small enough to do exact integer computation,
319 * return the log of the exact value
320 */
321 if (n < 67) {
322 return Math.log(binomialCoefficient(n,k));
323 }
324
325 /*
326 * Return the log of binomialCoefficientDouble for values that will not
327 * overflow binomialCoefficientDouble
328 */
329 if (n < 1030) {
330 return Math.log(binomialCoefficientDouble(n, k));
331 }
332
333 if (k > n / 2) {
334 return binomialCoefficientLog(n, n - k);
335 }
336
337 /*
338 * Sum logs for values that could overflow
339 */
340 double logSum = 0;
341
342 // n!/(n-k)!
343 for (int i = n - k + 1; i <= n; i++) {
344 logSum += Math.log(i);
345 }
346
347 // divide by k!
348 for (int i = 2; i <= k; i++) {
349 logSum -= Math.log(i);
350 }
351
352 return logSum;
353 }
354
355 /**
356 * Check binomial preconditions.
357 * @param n the size of the set
358 * @param k the size of the subsets to be counted
359 * @exception IllegalArgumentException if preconditions are not met.
360 */
361 private static void checkBinomial(final int n, final int k)
362 throws IllegalArgumentException {
363 if (n < k) {
364 throw MathRuntimeException.createIllegalArgumentException(
365 "must have n >= k for binomial coefficient (n,k), got n = {0}, k = {1}",
366 n, k);
367 }
368 if (n < 0) {
369 throw MathRuntimeException.createIllegalArgumentException(
370 "must have n >= 0 for binomial coefficient (n,k), got n = {0}",
371 n);
372 }
373 }
374
375 /**
376 * Compares two numbers given some amount of allowed error.
377 *
378 * @param x the first number
379 * @param y the second number
380 * @param eps the amount of error to allow when checking for equality
381 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li>
382 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li>
383 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul>
384 */
385 public static int compareTo(double x, double y, double eps) {
386 if (equals(x, y, eps)) {
387 return 0;
388 } else if (x < y) {
389 return -1;
390 }
391 return 1;
392 }
393
394 /**
395 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
396 * hyperbolic cosine</a> of x.
397 *
398 * @param x double value for which to find the hyperbolic cosine
399 * @return hyperbolic cosine of x
400 */
401 public static double cosh(double x) {
402 return (Math.exp(x) + Math.exp(-x)) / 2.0;
403 }
404
405 /**
406 * Returns true iff both arguments are NaN or neither is NaN and they are
407 * equal
408 *
409 * @param x first value
410 * @param y second value
411 * @return true if the values are equal or both are NaN
412 */
413 public static boolean equals(double x, double y) {
414 return (Double.isNaN(x) && Double.isNaN(y)) || x == y;
415 }
416
417 /**
418 * Returns true iff both arguments are equal or within the range of allowed
419 * error (inclusive).
420 * <p>
421 * Two NaNs are considered equals, as are two infinities with same sign.
422 * </p>
423 *
424 * @param x first value
425 * @param y second value
426 * @param eps the amount of absolute error to allow
427 * @return true if the values are equal or within range of each other
428 */
429 public static boolean equals(double x, double y, double eps) {
430 return equals(x, y) || (Math.abs(y - x) <= eps);
431 }
432
433 /**
434 * Returns true iff both arguments are equal or within the range of allowed
435 * error (inclusive).
436 * Adapted from <a
437 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
438 * Bruce Dawson</a>
439 *
440 * @param x first value
441 * @param y second value
442 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
443 * values between {@code x} and {@code y}.
444 * @return {@code true} if there are less than {@code maxUlps} floating
445 * point values between {@code x} and {@code y}
446 */
447 public static boolean equals(double x, double y, int maxUlps) {
448 // Check that "maxUlps" is non-negative and small enough so that the
449 // default NAN won't compare as equal to anything.
450 assert maxUlps > 0 && maxUlps < NAN_GAP;
451
452 long xInt = Double.doubleToLongBits(x);
453 long yInt = Double.doubleToLongBits(y);
454
455 // Make lexicographically ordered as a two's-complement integer.
456 if (xInt < 0) {
457 xInt = SGN_MASK - xInt;
458 }
459 if (yInt < 0) {
460 yInt = SGN_MASK - yInt;
461 }
462
463 return Math.abs(xInt - yInt) <= maxUlps;
464 }
465
466 /**
467 * Returns true iff both arguments are null or have same dimensions
468 * and all their elements are {@link #equals(double,double) equals}
469 *
470 * @param x first array
471 * @param y second array
472 * @return true if the values are both null or have same dimension
473 * and equal elements
474 * @since 1.2
475 */
476 public static boolean equals(double[] x, double[] y) {
477 if ((x == null) || (y == null)) {
478 return !((x == null) ^ (y == null));
479 }
480 if (x.length != y.length) {
481 return false;
482 }
483 for (int i = 0; i < x.length; ++i) {
484 if (!equals(x[i], y[i])) {
485 return false;
486 }
487 }
488 return true;
489 }
490
491 /**
492 * Returns n!. Shorthand for <code>n</code> <a
493 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
494 * product of the numbers <code>1,...,n</code>.
495 * <p>
496 * <Strong>Preconditions</strong>:
497 * <ul>
498 * <li> <code>n >= 0</code> (otherwise
499 * <code>IllegalArgumentException</code> is thrown)</li>
500 * <li> The result is small enough to fit into a <code>long</code>. The
501 * largest value of <code>n</code> for which <code>n!</code> <
502 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
503 * an <code>ArithMeticException </code> is thrown.</li>
504 * </ul>
505 * </p>
506 *
507 * @param n argument
508 * @return <code>n!</code>
509 * @throws ArithmeticException if the result is too large to be represented
510 * by a long integer.
511 * @throws IllegalArgumentException if n < 0
512 */
513 public static long factorial(final int n) {
514 if (n < 0) {
515 throw MathRuntimeException.createIllegalArgumentException(
516 "must have n >= 0 for n!, got n = {0}",
517 n);
518 }
519 if (n > 20) {
520 throw new ArithmeticException(
521 "factorial value is too large to fit in a long");
522 }
523 return FACTORIALS[n];
524 }
525
526 /**
527 * Returns n!. Shorthand for <code>n</code> <a
528 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
529 * product of the numbers <code>1,...,n</code> as a <code>double</code>.
530 * <p>
531 * <Strong>Preconditions</strong>:
532 * <ul>
533 * <li> <code>n >= 0</code> (otherwise
534 * <code>IllegalArgumentException</code> is thrown)</li>
535 * <li> The result is small enough to fit into a <code>double</code>. The
536 * largest value of <code>n</code> for which <code>n!</code> <
537 * Double.MAX_VALUE</code> is 170. If the computed value exceeds
538 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
539 * </ul>
540 * </p>
541 *
542 * @param n argument
543 * @return <code>n!</code>
544 * @throws IllegalArgumentException if n < 0
545 */
546 public static double factorialDouble(final int n) {
547 if (n < 0) {
548 throw MathRuntimeException.createIllegalArgumentException(
549 "must have n >= 0 for n!, got n = {0}",
550 n);
551 }
552 if (n < 21) {
553 return factorial(n);
554 }
555 return Math.floor(Math.exp(factorialLog(n)) + 0.5);
556 }
557
558 /**
559 * Returns the natural logarithm of n!.
560 * <p>
561 * <Strong>Preconditions</strong>:
562 * <ul>
563 * <li> <code>n >= 0</code> (otherwise
564 * <code>IllegalArgumentException</code> is thrown)</li>
565 * </ul></p>
566 *
567 * @param n argument
568 * @return <code>n!</code>
569 * @throws IllegalArgumentException if preconditions are not met.
570 */
571 public static double factorialLog(final int n) {
572 if (n < 0) {
573 throw MathRuntimeException.createIllegalArgumentException(
574 "must have n >= 0 for n!, got n = {0}",
575 n);
576 }
577 if (n < 21) {
578 return Math.log(factorial(n));
579 }
580 double logSum = 0;
581 for (int i = 2; i <= n; i++) {
582 logSum += Math.log(i);
583 }
584 return logSum;
585 }
586
587 /**
588 * <p>
589 * Gets the greatest common divisor of the absolute value of two numbers,
590 * using the "binary gcd" method which avoids division and modulo
591 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
592 * Stein (1961).
593 * </p>
594 * Special cases:
595 * <ul>
596 * <li>The invocations
597 * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>,
598 * <code>gcd(Integer.MIN_VALUE, 0)</code> and
599 * <code>gcd(0, Integer.MIN_VALUE)</code> throw an
600 * <code>ArithmeticException</code>, because the result would be 2^31, which
601 * is too large for an int value.</li>
602 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and
603 * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except
604 * for the special cases above.
605 * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
606 * <code>0</code>.</li>
607 * </ul>
608 *
609 * @param p any number
610 * @param q any number
611 * @return the greatest common divisor, never negative
612 * @throws ArithmeticException if the result cannot be represented as a
613 * nonnegative int value
614 * @since 1.1
615 */
616 public static int gcd(final int p, final int q) {
617 int u = p;
618 int v = q;
619 if ((u == 0) || (v == 0)) {
620 if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
621 throw MathRuntimeException.createArithmeticException(
622 "overflow: gcd({0}, {1}) is 2^31",
623 p, q);
624 }
625 return Math.abs(u) + Math.abs(v);
626 }
627 // keep u and v negative, as negative integers range down to
628 // -2^31, while positive numbers can only be as large as 2^31-1
629 // (i.e. we can't necessarily negate a negative number without
630 // overflow)
631 /* assert u!=0 && v!=0; */
632 if (u > 0) {
633 u = -u;
634 } // make u negative
635 if (v > 0) {
636 v = -v;
637 } // make v negative
638 // B1. [Find power of 2]
639 int k = 0;
640 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are
641 // both even...
642 u /= 2;
643 v /= 2;
644 k++; // cast out twos.
645 }
646 if (k == 31) {
647 throw MathRuntimeException.createArithmeticException(
648 "overflow: gcd({0}, {1}) is 2^31",
649 p, q);
650 }
651 // B2. Initialize: u and v have been divided by 2^k and at least
652 // one is odd.
653 int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
654 // t negative: u was odd, v may be even (t replaces v)
655 // t positive: u was even, v is odd (t replaces u)
656 do {
657 /* assert u<0 && v<0; */
658 // B4/B3: cast out twos from t.
659 while ((t & 1) == 0) { // while t is even..
660 t /= 2; // cast out twos
661 }
662 // B5 [reset max(u,v)]
663 if (t > 0) {
664 u = -t;
665 } else {
666 v = t;
667 }
668 // B6/B3. at this point both u and v should be odd.
669 t = (v - u) / 2;
670 // |u| larger: t positive (replace u)
671 // |v| larger: t negative (replace v)
672 } while (t != 0);
673 return -u * (1 << k); // gcd is u*2^k
674 }
675
676 /**
677 * <p>
678 * Gets the greatest common divisor of the absolute value of two numbers,
679 * using the "binary gcd" method which avoids division and modulo
680 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
681 * Stein (1961).
682 * </p>
683 * Special cases:
684 * <ul>
685 * <li>The invocations
686 * <code>gcd(Long.MIN_VALUE, Long.MIN_VALUE)</code>,
687 * <code>gcd(Long.MIN_VALUE, 0L)</code> and
688 * <code>gcd(0L, Long.MIN_VALUE)</code> throw an
689 * <code>ArithmeticException</code>, because the result would be 2^63, which
690 * is too large for a long value.</li>
691 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0L, x)</code> and
692 * <code>gcd(x, 0L)</code> is the absolute value of <code>x</code>, except
693 * for the special cases above.
694 * <li>The invocation <code>gcd(0L, 0L)</code> is the only one which returns
695 * <code>0L</code>.</li>
696 * </ul>
697 *
698 * @param p any number
699 * @param q any number
700 * @return the greatest common divisor, never negative
701 * @throws ArithmeticException if the result cannot be represented as a nonnegative long
702 * value
703 * @since 2.1
704 */
705 public static long gcd(final long p, final long q) {
706 long u = p;
707 long v = q;
708 if ((u == 0) || (v == 0)) {
709 if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){
710 throw MathRuntimeException.createArithmeticException(
711 "overflow: gcd({0}, {1}) is 2^63",
712 p, q);
713 }
714 return Math.abs(u) + Math.abs(v);
715 }
716 // keep u and v negative, as negative integers range down to
717 // -2^63, while positive numbers can only be as large as 2^63-1
718 // (i.e. we can't necessarily negate a negative number without
719 // overflow)
720 /* assert u!=0 && v!=0; */
721 if (u > 0) {
722 u = -u;
723 } // make u negative
724 if (v > 0) {
725 v = -v;
726 } // make v negative
727 // B1. [Find power of 2]
728 int k = 0;
729 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
730 // both even...
731 u /= 2;
732 v /= 2;
733 k++; // cast out twos.
734 }
735 if (k == 63) {
736 throw MathRuntimeException.createArithmeticException(
737 "overflow: gcd({0}, {1}) is 2^63",
738 p, q);
739 }
740 // B2. Initialize: u and v have been divided by 2^k and at least
741 // one is odd.
742 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
743 // t negative: u was odd, v may be even (t replaces v)
744 // t positive: u was even, v is odd (t replaces u)
745 do {
746 /* assert u<0 && v<0; */
747 // B4/B3: cast out twos from t.
748 while ((t & 1) == 0) { // while t is even..
749 t /= 2; // cast out twos
750 }
751 // B5 [reset max(u,v)]
752 if (t > 0) {
753 u = -t;
754 } else {
755 v = t;
756 }
757 // B6/B3. at this point both u and v should be odd.
758 t = (v - u) / 2;
759 // |u| larger: t positive (replace u)
760 // |v| larger: t negative (replace v)
761 } while (t != 0);
762 return -u * (1L << k); // gcd is u*2^k
763 }
764
765 /**
766 * Returns an integer hash code representing the given double value.
767 *
768 * @param value the value to be hashed
769 * @return the hash code
770 */
771 public static int hash(double value) {
772 return new Double(value).hashCode();
773 }
774
775 /**
776 * Returns an integer hash code representing the given double array.
777 *
778 * @param value the value to be hashed (may be null)
779 * @return the hash code
780 * @since 1.2
781 */
782 public static int hash(double[] value) {
783 return Arrays.hashCode(value);
784 }
785
786 /**
787 * For a byte value x, this method returns (byte)(+1) if x >= 0 and
788 * (byte)(-1) if x < 0.
789 *
790 * @param x the value, a byte
791 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
792 */
793 public static byte indicator(final byte x) {
794 return (x >= ZB) ? PB : NB;
795 }
796
797 /**
798 * For a double precision value x, this method returns +1.0 if x >= 0 and
799 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
800 * <code>NaN</code>.
801 *
802 * @param x the value, a double
803 * @return +1.0 or -1.0, depending on the sign of x
804 */
805 public static double indicator(final double x) {
806 if (Double.isNaN(x)) {
807 return Double.NaN;
808 }
809 return (x >= 0.0) ? 1.0 : -1.0;
810 }
811
812 /**
813 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
814 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
815 *
816 * @param x the value, a float
817 * @return +1.0F or -1.0F, depending on the sign of x
818 */
819 public static float indicator(final float x) {
820 if (Float.isNaN(x)) {
821 return Float.NaN;
822 }
823 return (x >= 0.0F) ? 1.0F : -1.0F;
824 }
825
826 /**
827 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
828 *
829 * @param x the value, an int
830 * @return +1 or -1, depending on the sign of x
831 */
832 public static int indicator(final int x) {
833 return (x >= 0) ? 1 : -1;
834 }
835
836 /**
837 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
838 *
839 * @param x the value, a long
840 * @return +1L or -1L, depending on the sign of x
841 */
842 public static long indicator(final long x) {
843 return (x >= 0L) ? 1L : -1L;
844 }
845
846 /**
847 * For a short value x, this method returns (short)(+1) if x >= 0 and
848 * (short)(-1) if x < 0.
849 *
850 * @param x the value, a short
851 * @return (short)(+1) or (short)(-1), depending on the sign of x
852 */
853 public static short indicator(final short x) {
854 return (x >= ZS) ? PS : NS;
855 }
856
857 /**
858 * <p>
859 * Returns the least common multiple of the absolute value of two numbers,
860 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
861 * </p>
862 * Special cases:
863 * <ul>
864 * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and
865 * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a
866 * power of 2, throw an <code>ArithmeticException</code>, because the result
867 * would be 2^31, which is too large for an int value.</li>
868 * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is
869 * <code>0</code> for any <code>x</code>.
870 * </ul>
871 *
872 * @param a any number
873 * @param b any number
874 * @return the least common multiple, never negative
875 * @throws ArithmeticException
876 * if the result cannot be represented as a nonnegative int
877 * value
878 * @since 1.1
879 */
880 public static int lcm(int a, int b) {
881 if (a==0 || b==0){
882 return 0;
883 }
884 int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b));
885 if (lcm == Integer.MIN_VALUE) {
886 throw MathRuntimeException.createArithmeticException(
887 "overflow: lcm({0}, {1}) is 2^31",
888 a, b);
889 }
890 return lcm;
891 }
892
893 /**
894 * <p>
895 * Returns the least common multiple of the absolute value of two numbers,
896 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
897 * </p>
898 * Special cases:
899 * <ul>
900 * <li>The invocations <code>lcm(Long.MIN_VALUE, n)</code> and
901 * <code>lcm(n, Long.MIN_VALUE)</code>, where <code>abs(n)</code> is a
902 * power of 2, throw an <code>ArithmeticException</code>, because the result
903 * would be 2^63, which is too large for an int value.</li>
904 * <li>The result of <code>lcm(0L, x)</code> and <code>lcm(x, 0L)</code> is
905 * <code>0L</code> for any <code>x</code>.
906 * </ul>
907 *
908 * @param a any number
909 * @param b any number
910 * @return the least common multiple, never negative
911 * @throws ArithmeticException if the result cannot be represented as a nonnegative long
912 * value
913 * @since 2.1
914 */
915 public static long lcm(long a, long b) {
916 if (a==0 || b==0){
917 return 0;
918 }
919 long lcm = Math.abs(mulAndCheck(a / gcd(a, b), b));
920 if (lcm == Long.MIN_VALUE){
921 throw MathRuntimeException.createArithmeticException(
922 "overflow: lcm({0}, {1}) is 2^63",
923 a, b);
924 }
925 return lcm;
926 }
927
928 /**
929 * <p>Returns the
930 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
931 * for base <code>b</code> of <code>x</code>.
932 * </p>
933 * <p>Returns <code>NaN<code> if either argument is negative. If
934 * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
935 * If <code>base</code> is positive and <code>x</code> is 0,
936 * <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments
937 * are 0, the result is <code>NaN</code>.</p>
938 *
939 * @param base the base of the logarithm, must be greater than 0
940 * @param x argument, must be greater than 0
941 * @return the value of the logarithm - the number y such that base^y = x.
942 * @since 1.2
943 */
944 public static double log(double base, double x) {
945 return Math.log(x)/Math.log(base);
946 }
947
948 /**
949 * Multiply two integers, checking for overflow.
950 *
951 * @param x a factor
952 * @param y a factor
953 * @return the product <code>x*y</code>
954 * @throws ArithmeticException if the result can not be represented as an
955 * int
956 * @since 1.1
957 */
958 public static int mulAndCheck(int x, int y) {
959 long m = ((long)x) * ((long)y);
960 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
961 throw new ArithmeticException("overflow: mul");
962 }
963 return (int)m;
964 }
965
966 /**
967 * Multiply two long integers, checking for overflow.
968 *
969 * @param a first value
970 * @param b second value
971 * @return the product <code>a * b</code>
972 * @throws ArithmeticException if the result can not be represented as an
973 * long
974 * @since 1.2
975 */
976 public static long mulAndCheck(long a, long b) {
977 long ret;
978 String msg = "overflow: multiply";
979 if (a > b) {
980 // use symmetry to reduce boundary cases
981 ret = mulAndCheck(b, a);
982 } else {
983 if (a < 0) {
984 if (b < 0) {
985 // check for positive overflow with negative a, negative b
986 if (a >= Long.MAX_VALUE / b) {
987 ret = a * b;
988 } else {
989 throw new ArithmeticException(msg);
990 }
991 } else if (b > 0) {
992 // check for negative overflow with negative a, positive b
993 if (Long.MIN_VALUE / b <= a) {
994 ret = a * b;
995 } else {
996 throw new ArithmeticException(msg);
997
998 }
999 } else {
1000 // assert b == 0
1001 ret = 0;
1002 }
1003 } else if (a > 0) {
1004 // assert a > 0
1005 // assert b > 0
1006
1007 // check for positive overflow with positive a, positive b
1008 if (a <= Long.MAX_VALUE / b) {
1009 ret = a * b;
1010 } else {
1011 throw new ArithmeticException(msg);
1012 }
1013 } else {
1014 // assert a == 0
1015 ret = 0;
1016 }
1017 }
1018 return ret;
1019 }
1020
1021 /**
1022 * Get the next machine representable number after a number, moving
1023 * in the direction of another number.
1024 * <p>
1025 * If <code>direction</code> is greater than or equal to<code>d</code>,
1026 * the smallest machine representable number strictly greater than
1027 * <code>d</code> is returned; otherwise the largest representable number
1028 * strictly less than <code>d</code> is returned.</p>
1029 * <p>
1030 * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
1031 *
1032 * @param d base number
1033 * @param direction (the only important thing is whether
1034 * direction is greater or smaller than d)
1035 * @return the next machine representable number in the specified direction
1036 * @since 1.2
1037 */
1038 public static double nextAfter(double d, double direction) {
1039
1040 // handling of some important special cases
1041 if (Double.isNaN(d) || Double.isInfinite(d)) {
1042 return d;
1043 } else if (d == 0) {
1044 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
1045 }
1046 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
1047 // are handled just as normal numbers
1048
1049 // split the double in raw components
1050 long bits = Double.doubleToLongBits(d);
1051 long sign = bits & 0x8000000000000000L;
1052 long exponent = bits & 0x7ff0000000000000L;
1053 long mantissa = bits & 0x000fffffffffffffL;
1054
1055 if (d * (direction - d) >= 0) {
1056 // we should increase the mantissa
1057 if (mantissa == 0x000fffffffffffffL) {
1058 return Double.longBitsToDouble(sign |
1059 (exponent + 0x0010000000000000L));
1060 } else {
1061 return Double.longBitsToDouble(sign |
1062 exponent | (mantissa + 1));
1063 }
1064 } else {
1065 // we should decrease the mantissa
1066 if (mantissa == 0L) {
1067 return Double.longBitsToDouble(sign |
1068 (exponent - 0x0010000000000000L) |
1069 0x000fffffffffffffL);
1070 } else {
1071 return Double.longBitsToDouble(sign |
1072 exponent | (mantissa - 1));
1073 }
1074 }
1075
1076 }
1077
1078 /**
1079 * Scale a number by 2<sup>scaleFactor</sup>.
1080 * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p>
1081 *
1082 * @param d base number
1083 * @param scaleFactor power of two by which d sould be multiplied
1084 * @return d × 2<sup>scaleFactor</sup>
1085 * @since 2.0
1086 */
1087 public static double scalb(final double d, final int scaleFactor) {
1088
1089 // handling of some important special cases
1090 if ((d == 0) || Double.isNaN(d) || Double.isInfinite(d)) {
1091 return d;
1092 }
1093
1094 // split the double in raw components
1095 final long bits = Double.doubleToLongBits(d);
1096 final long exponent = bits & 0x7ff0000000000000L;
1097 final long rest = bits & 0x800fffffffffffffL;
1098
1099 // shift the exponent
1100 final long newBits = rest | (exponent + (((long) scaleFactor) << 52));
1101 return Double.longBitsToDouble(newBits);
1102
1103 }
1104
1105 /**
1106 * Normalize an angle in a 2&pi wide interval around a center value.
1107 * <p>This method has three main uses:</p>
1108 * <ul>
1109 * <li>normalize an angle between 0 and 2π:<br/>
1110 * <code>a = MathUtils.normalizeAngle(a, Math.PI);</code></li>
1111 * <li>normalize an angle between -π and +π<br/>
1112 * <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li>
1113 * <li>compute the angle between two defining angular positions:<br>
1114 * <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li>
1115 * </ul>
1116 * <p>Note that due to numerical accuracy and since π cannot be represented
1117 * exactly, the result interval is <em>closed</em>, it cannot be half-closed
1118 * as would be more satisfactory in a purely mathematical view.</p>
1119 * @param a angle to normalize
1120 * @param center center of the desired 2π interval for the result
1121 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π
1122 * @since 1.2
1123 */
1124 public static double normalizeAngle(double a, double center) {
1125 return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI);
1126 }
1127
1128 /**
1129 * <p>Normalizes an array to make it sum to a specified value.
1130 * Returns the result of the transformation <pre>
1131 * x |-> x * normalizedSum / sum
1132 * </pre>
1133 * applied to each non-NaN element x of the input array, where sum is the
1134 * sum of the non-NaN entries in the input array.</p>
1135 *
1136 * <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite
1137 * or NaN and ArithmeticException if the input array contains any infinite elements
1138 * or sums to 0</p>
1139 *
1140 * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
1141 *
1142 * @param values input array to be normalized
1143 * @param normalizedSum target sum for the normalized array
1144 * @return normalized array
1145 * @throws ArithmeticException if the input array contains infinite elements or sums to zero
1146 * @throws IllegalArgumentException if the target sum is infinite or NaN
1147 * @since 2.1
1148 */
1149 public static double[] normalizeArray(double[] values, double normalizedSum)
1150 throws ArithmeticException, IllegalArgumentException {
1151 if (Double.isInfinite(normalizedSum)) {
1152 throw MathRuntimeException.createIllegalArgumentException(
1153 "Cannot normalize to an infinite value");
1154 }
1155 if (Double.isNaN(normalizedSum)) {
1156 throw MathRuntimeException.createIllegalArgumentException(
1157 "Cannot normalize to NaN");
1158 }
1159 double sum = 0d;
1160 final int len = values.length;
1161 double[] out = new double[len];
1162 for (int i = 0; i < len; i++) {
1163 if (Double.isInfinite(values[i])) {
1164 throw MathRuntimeException.createArithmeticException(
1165 "Array contains an infinite element, {0} at index {1}", values[i], i);
1166 }
1167 if (!Double.isNaN(values[i])) {
1168 sum += values[i];
1169 }
1170 }
1171 if (sum == 0) {
1172 throw MathRuntimeException.createArithmeticException(
1173 "Array sums to zero");
1174 }
1175 for (int i = 0; i < len; i++) {
1176 if (Double.isNaN(values[i])) {
1177 out[i] = Double.NaN;
1178 } else {
1179 out[i] = values[i] * normalizedSum / sum;
1180 }
1181 }
1182 return out;
1183 }
1184
1185 /**
1186 * Round the given value to the specified number of decimal places. The
1187 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
1188 *
1189 * @param x the value to round.
1190 * @param scale the number of digits to the right of the decimal point.
1191 * @return the rounded value.
1192 * @since 1.1
1193 */
1194 public static double round(double x, int scale) {
1195 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1196 }
1197
1198 /**
1199 * Round the given value to the specified number of decimal places. The
1200 * value is rounded using the given method which is any method defined in
1201 * {@link BigDecimal}.
1202 *
1203 * @param x the value to round.
1204 * @param scale the number of digits to the right of the decimal point.
1205 * @param roundingMethod the rounding method as defined in
1206 * {@link BigDecimal}.
1207 * @return the rounded value.
1208 * @since 1.1
1209 */
1210 public static double round(double x, int scale, int roundingMethod) {
1211 try {
1212 return (new BigDecimal
1213 (Double.toString(x))
1214 .setScale(scale, roundingMethod))
1215 .doubleValue();
1216 } catch (NumberFormatException ex) {
1217 if (Double.isInfinite(x)) {
1218 return x;
1219 } else {
1220 return Double.NaN;
1221 }
1222 }
1223 }
1224
1225 /**
1226 * Round the given value to the specified number of decimal places. The
1227 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
1228 *
1229 * @param x the value to round.
1230 * @param scale the number of digits to the right of the decimal point.
1231 * @return the rounded value.
1232 * @since 1.1
1233 */
1234 public static float round(float x, int scale) {
1235 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1236 }
1237
1238 /**
1239 * Round the given value to the specified number of decimal places. The
1240 * value is rounded using the given method which is any method defined in
1241 * {@link BigDecimal}.
1242 *
1243 * @param x the value to round.
1244 * @param scale the number of digits to the right of the decimal point.
1245 * @param roundingMethod the rounding method as defined in
1246 * {@link BigDecimal}.
1247 * @return the rounded value.
1248 * @since 1.1
1249 */
1250 public static float round(float x, int scale, int roundingMethod) {
1251 float sign = indicator(x);
1252 float factor = (float)Math.pow(10.0f, scale) * sign;
1253 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
1254 }
1255
1256 /**
1257 * Round the given non-negative, value to the "nearest" integer. Nearest is
1258 * determined by the rounding method specified. Rounding methods are defined
1259 * in {@link BigDecimal}.
1260 *
1261 * @param unscaled the value to round.
1262 * @param sign the sign of the original, scaled value.
1263 * @param roundingMethod the rounding method as defined in
1264 * {@link BigDecimal}.
1265 * @return the rounded value.
1266 * @since 1.1
1267 */
1268 private static double roundUnscaled(double unscaled, double sign,
1269 int roundingMethod) {
1270 switch (roundingMethod) {
1271 case BigDecimal.ROUND_CEILING :
1272 if (sign == -1) {
1273 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1274 } else {
1275 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1276 }
1277 break;
1278 case BigDecimal.ROUND_DOWN :
1279 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1280 break;
1281 case BigDecimal.ROUND_FLOOR :
1282 if (sign == -1) {
1283 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1284 } else {
1285 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1286 }
1287 break;
1288 case BigDecimal.ROUND_HALF_DOWN : {
1289 unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY);
1290 double fraction = unscaled - Math.floor(unscaled);
1291 if (fraction > 0.5) {
1292 unscaled = Math.ceil(unscaled);
1293 } else {
1294 unscaled = Math.floor(unscaled);
1295 }
1296 break;
1297 }
1298 case BigDecimal.ROUND_HALF_EVEN : {
1299 double fraction = unscaled - Math.floor(unscaled);
1300 if (fraction > 0.5) {
1301 unscaled = Math.ceil(unscaled);
1302 } else if (fraction < 0.5) {
1303 unscaled = Math.floor(unscaled);
1304 } else {
1305 // The following equality test is intentional and needed for rounding purposes
1306 if (Math.floor(unscaled) / 2.0 == Math.floor(Math
1307 .floor(unscaled) / 2.0)) { // even
1308 unscaled = Math.floor(unscaled);
1309 } else { // odd
1310 unscaled = Math.ceil(unscaled);
1311 }
1312 }
1313 break;
1314 }
1315 case BigDecimal.ROUND_HALF_UP : {
1316 unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY);
1317 double fraction = unscaled - Math.floor(unscaled);
1318 if (fraction >= 0.5) {
1319 unscaled = Math.ceil(unscaled);
1320 } else {
1321 unscaled = Math.floor(unscaled);
1322 }
1323 break;
1324 }
1325 case BigDecimal.ROUND_UNNECESSARY :
1326 if (unscaled != Math.floor(unscaled)) {
1327 throw new ArithmeticException("Inexact result from rounding");
1328 }
1329 break;
1330 case BigDecimal.ROUND_UP :
1331 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1332 break;
1333 default :
1334 throw MathRuntimeException.createIllegalArgumentException(
1335 "invalid rounding method {0}, valid methods: {1} ({2}), {3} ({4})," +
1336 " {5} ({6}), {7} ({8}), {9} ({10}), {11} ({12}), {13} ({14}), {15} ({16})",
1337 roundingMethod,
1338 "ROUND_CEILING", BigDecimal.ROUND_CEILING,
1339 "ROUND_DOWN", BigDecimal.ROUND_DOWN,
1340 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
1341 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
1342 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
1343 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
1344 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
1345 "ROUND_UP", BigDecimal.ROUND_UP);
1346 }
1347 return unscaled;
1348 }
1349
1350 /**
1351 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1352 * for byte value <code>x</code>.
1353 * <p>
1354 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
1355 * x = 0, and (byte)(-1) if x < 0.</p>
1356 *
1357 * @param x the value, a byte
1358 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
1359 */
1360 public static byte sign(final byte x) {
1361 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
1362 }
1363
1364 /**
1365 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1366 * for double precision <code>x</code>.
1367 * <p>
1368 * For a double value <code>x</code>, this method returns
1369 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
1370 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
1371 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
1372 *
1373 * @param x the value, a double
1374 * @return +1.0, 0.0, or -1.0, depending on the sign of x
1375 */
1376 public static double sign(final double x) {
1377 if (Double.isNaN(x)) {
1378 return Double.NaN;
1379 }
1380 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
1381 }
1382
1383 /**
1384 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1385 * for float value <code>x</code>.
1386 * <p>
1387 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
1388 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
1389 * is <code>NaN</code>.</p>
1390 *
1391 * @param x the value, a float
1392 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
1393 */
1394 public static float sign(final float x) {
1395 if (Float.isNaN(x)) {
1396 return Float.NaN;
1397 }
1398 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
1399 }
1400
1401 /**
1402 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1403 * for int value <code>x</code>.
1404 * <p>
1405 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
1406 * if x < 0.</p>
1407 *
1408 * @param x the value, an int
1409 * @return +1, 0, or -1, depending on the sign of x
1410 */
1411 public static int sign(final int x) {
1412 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
1413 }
1414
1415 /**
1416 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1417 * for long value <code>x</code>.
1418 * <p>
1419 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
1420 * -1L if x < 0.</p>
1421 *
1422 * @param x the value, a long
1423 * @return +1L, 0L, or -1L, depending on the sign of x
1424 */
1425 public static long sign(final long x) {
1426 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
1427 }
1428
1429 /**
1430 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1431 * for short value <code>x</code>.
1432 * <p>
1433 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
1434 * if x = 0, and (short)(-1) if x < 0.</p>
1435 *
1436 * @param x the value, a short
1437 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
1438 * x
1439 */
1440 public static short sign(final short x) {
1441 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
1442 }
1443
1444 /**
1445 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
1446 * hyperbolic sine</a> of x.
1447 *
1448 * @param x double value for which to find the hyperbolic sine
1449 * @return hyperbolic sine of x
1450 */
1451 public static double sinh(double x) {
1452 return (Math.exp(x) - Math.exp(-x)) / 2.0;
1453 }
1454
1455 /**
1456 * Subtract two integers, checking for overflow.
1457 *
1458 * @param x the minuend
1459 * @param y the subtrahend
1460 * @return the difference <code>x-y</code>
1461 * @throws ArithmeticException if the result can not be represented as an
1462 * int
1463 * @since 1.1
1464 */
1465 public static int subAndCheck(int x, int y) {
1466 long s = (long)x - (long)y;
1467 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
1468 throw new ArithmeticException("overflow: subtract");
1469 }
1470 return (int)s;
1471 }
1472
1473 /**
1474 * Subtract two long integers, checking for overflow.
1475 *
1476 * @param a first value
1477 * @param b second value
1478 * @return the difference <code>a-b</code>
1479 * @throws ArithmeticException if the result can not be represented as an
1480 * long
1481 * @since 1.2
1482 */
1483 public static long subAndCheck(long a, long b) {
1484 long ret;
1485 String msg = "overflow: subtract";
1486 if (b == Long.MIN_VALUE) {
1487 if (a < 0) {
1488 ret = a - b;
1489 } else {
1490 throw new ArithmeticException(msg);
1491 }
1492 } else {
1493 // use additive inverse
1494 ret = addAndCheck(a, -b, msg);
1495 }
1496 return ret;
1497 }
1498
1499 /**
1500 * Raise an int to an int power.
1501 * @param k number to raise
1502 * @param e exponent (must be positive or null)
1503 * @return k<sup>e</sup>
1504 * @exception IllegalArgumentException if e is negative
1505 */
1506 public static int pow(final int k, int e)
1507 throws IllegalArgumentException {
1508
1509 if (e < 0) {
1510 throw MathRuntimeException.createIllegalArgumentException(
1511 "cannot raise an integral value to a negative power ({0}^{1})",
1512 k, e);
1513 }
1514
1515 int result = 1;
1516 int k2p = k;
1517 while (e != 0) {
1518 if ((e & 0x1) != 0) {
1519 result *= k2p;
1520 }
1521 k2p *= k2p;
1522 e = e >> 1;
1523 }
1524
1525 return result;
1526
1527 }
1528
1529 /**
1530 * Raise an int to a long power.
1531 * @param k number to raise
1532 * @param e exponent (must be positive or null)
1533 * @return k<sup>e</sup>
1534 * @exception IllegalArgumentException if e is negative
1535 */
1536 public static int pow(final int k, long e)
1537 throws IllegalArgumentException {
1538
1539 if (e < 0) {
1540 throw MathRuntimeException.createIllegalArgumentException(
1541 "cannot raise an integral value to a negative power ({0}^{1})",
1542 k, e);
1543 }
1544
1545 int result = 1;
1546 int k2p = k;
1547 while (e != 0) {
1548 if ((e & 0x1) != 0) {
1549 result *= k2p;
1550 }
1551 k2p *= k2p;
1552 e = e >> 1;
1553 }
1554
1555 return result;
1556
1557 }
1558
1559 /**
1560 * Raise a long to an int power.
1561 * @param k number to raise
1562 * @param e exponent (must be positive or null)
1563 * @return k<sup>e</sup>
1564 * @exception IllegalArgumentException if e is negative
1565 */
1566 public static long pow(final long k, int e)
1567 throws IllegalArgumentException {
1568
1569 if (e < 0) {
1570 throw MathRuntimeException.createIllegalArgumentException(
1571 "cannot raise an integral value to a negative power ({0}^{1})",
1572 k, e);
1573 }
1574
1575 long result = 1l;
1576 long k2p = k;
1577 while (e != 0) {
1578 if ((e & 0x1) != 0) {
1579 result *= k2p;
1580 }
1581 k2p *= k2p;
1582 e = e >> 1;
1583 }
1584
1585 return result;
1586
1587 }
1588
1589 /**
1590 * Raise a long to a long power.
1591 * @param k number to raise
1592 * @param e exponent (must be positive or null)
1593 * @return k<sup>e</sup>
1594 * @exception IllegalArgumentException if e is negative
1595 */
1596 public static long pow(final long k, long e)
1597 throws IllegalArgumentException {
1598
1599 if (e < 0) {
1600 throw MathRuntimeException.createIllegalArgumentException(
1601 "cannot raise an integral value to a negative power ({0}^{1})",
1602 k, e);
1603 }
1604
1605 long result = 1l;
1606 long k2p = k;
1607 while (e != 0) {
1608 if ((e & 0x1) != 0) {
1609 result *= k2p;
1610 }
1611 k2p *= k2p;
1612 e = e >> 1;
1613 }
1614
1615 return result;
1616
1617 }
1618
1619 /**
1620 * Raise a BigInteger to an int power.
1621 * @param k number to raise
1622 * @param e exponent (must be positive or null)
1623 * @return k<sup>e</sup>
1624 * @exception IllegalArgumentException if e is negative
1625 */
1626 public static BigInteger pow(final BigInteger k, int e)
1627 throws IllegalArgumentException {
1628
1629 if (e < 0) {
1630 throw MathRuntimeException.createIllegalArgumentException(
1631 "cannot raise an integral value to a negative power ({0}^{1})",
1632 k, e);
1633 }
1634
1635 return k.pow(e);
1636
1637 }
1638
1639 /**
1640 * Raise a BigInteger to a long power.
1641 * @param k number to raise
1642 * @param e exponent (must be positive or null)
1643 * @return k<sup>e</sup>
1644 * @exception IllegalArgumentException if e is negative
1645 */
1646 public static BigInteger pow(final BigInteger k, long e)
1647 throws IllegalArgumentException {
1648
1649 if (e < 0) {
1650 throw MathRuntimeException.createIllegalArgumentException(
1651 "cannot raise an integral value to a negative power ({0}^{1})",
1652 k, e);
1653 }
1654
1655 BigInteger result = BigInteger.ONE;
1656 BigInteger k2p = k;
1657 while (e != 0) {
1658 if ((e & 0x1) != 0) {
1659 result = result.multiply(k2p);
1660 }
1661 k2p = k2p.multiply(k2p);
1662 e = e >> 1;
1663 }
1664
1665 return result;
1666
1667 }
1668
1669 /**
1670 * Raise a BigInteger to a BigInteger power.
1671 * @param k number to raise
1672 * @param e exponent (must be positive or null)
1673 * @return k<sup>e</sup>
1674 * @exception IllegalArgumentException if e is negative
1675 */
1676 public static BigInteger pow(final BigInteger k, BigInteger e)
1677 throws IllegalArgumentException {
1678
1679 if (e.compareTo(BigInteger.ZERO) < 0) {
1680 throw MathRuntimeException.createIllegalArgumentException(
1681 "cannot raise an integral value to a negative power ({0}^{1})",
1682 k, e);
1683 }
1684
1685 BigInteger result = BigInteger.ONE;
1686 BigInteger k2p = k;
1687 while (!BigInteger.ZERO.equals(e)) {
1688 if (e.testBit(0)) {
1689 result = result.multiply(k2p);
1690 }
1691 k2p = k2p.multiply(k2p);
1692 e = e.shiftRight(1);
1693 }
1694
1695 return result;
1696
1697 }
1698
1699 /**
1700 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1701 *
1702 * @param p1 the first point
1703 * @param p2 the second point
1704 * @return the L<sub>1</sub> distance between the two points
1705 */
1706 public static double distance1(double[] p1, double[] p2) {
1707 double sum = 0;
1708 for (int i = 0; i < p1.length; i++) {
1709 sum += Math.abs(p1[i] - p2[i]);
1710 }
1711 return sum;
1712 }
1713
1714 /**
1715 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1716 *
1717 * @param p1 the first point
1718 * @param p2 the second point
1719 * @return the L<sub>1</sub> distance between the two points
1720 */
1721 public static int distance1(int[] p1, int[] p2) {
1722 int sum = 0;
1723 for (int i = 0; i < p1.length; i++) {
1724 sum += Math.abs(p1[i] - p2[i]);
1725 }
1726 return sum;
1727 }
1728
1729 /**
1730 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1731 *
1732 * @param p1 the first point
1733 * @param p2 the second point
1734 * @return the L<sub>2</sub> distance between the two points
1735 */
1736 public static double distance(double[] p1, double[] p2) {
1737 double sum = 0;
1738 for (int i = 0; i < p1.length; i++) {
1739 final double dp = p1[i] - p2[i];
1740 sum += dp * dp;
1741 }
1742 return Math.sqrt(sum);
1743 }
1744
1745 /**
1746 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1747 *
1748 * @param p1 the first point
1749 * @param p2 the second point
1750 * @return the L<sub>2</sub> distance between the two points
1751 */
1752 public static double distance(int[] p1, int[] p2) {
1753 double sum = 0;
1754 for (int i = 0; i < p1.length; i++) {
1755 final double dp = p1[i] - p2[i];
1756 sum += dp * dp;
1757 }
1758 return Math.sqrt(sum);
1759 }
1760
1761 /**
1762 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1763 *
1764 * @param p1 the first point
1765 * @param p2 the second point
1766 * @return the L<sub>∞</sub> distance between the two points
1767 */
1768 public static double distanceInf(double[] p1, double[] p2) {
1769 double max = 0;
1770 for (int i = 0; i < p1.length; i++) {
1771 max = Math.max(max, Math.abs(p1[i] - p2[i]));
1772 }
1773 return max;
1774 }
1775
1776 /**
1777 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1778 *
1779 * @param p1 the first point
1780 * @param p2 the second point
1781 * @return the L<sub>∞</sub> distance between the two points
1782 */
1783 public static int distanceInf(int[] p1, int[] p2) {
1784 int max = 0;
1785 for (int i = 0; i < p1.length; i++) {
1786 max = Math.max(max, Math.abs(p1[i] - p2[i]));
1787 }
1788 return max;
1789 }
1790
1791 /**
1792 * Checks that the given array is sorted.
1793 *
1794 * @param val Values
1795 * @param dir Order direction (-1 for decreasing, 1 for increasing)
1796 * @param strict Whether the order should be strict
1797 * @throws IllegalArgumentException if the array is not sorted.
1798 */
1799 public static void checkOrder(double[] val, int dir, boolean strict) {
1800 double previous = val[0];
1801
1802 int max = val.length;
1803 for (int i = 1; i < max; i++) {
1804 if (dir > 0) {
1805 if (strict) {
1806 if (val[i] <= previous) {
1807 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not strictly increasing ({2} >= {3})",
1808 i - 1, i, previous, val[i]);
1809 }
1810 } else {
1811 if (val[i] < previous) {
1812 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not increasing ({2} > {3})",
1813 i - 1, i, previous, val[i]);
1814 }
1815 }
1816 } else {
1817 if (strict) {
1818 if (val[i] >= previous) {
1819 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not strictly decreasing ({2} <= {3})",
1820 i - 1, i, previous, val[i]);
1821 }
1822 } else {
1823 if (val[i] > previous) {
1824 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not decreasing ({2} < {3})",
1825 i - 1, i, previous, val[i]);
1826 }
1827 }
1828 }
1829
1830 previous = val[i];
1831 }
1832 }
1833 }