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