56
62
/** 0.0 cast as a short. */
57
63
private static final short ZS = (short)0;
60
private static final double TWO_PI = 2 * Math.PI;
62
65
/** Gap between NaN and regular numbers. */
63
66
private static final int NAN_GAP = 4 * 1024 * 1024;
65
68
/** Offset to order signed double numbers lexicographically. */
66
69
private static final long SGN_MASK = 0x8000000000000000L;
71
/** All long-representable factorials */
72
private static final long[] FACTORIALS = new long[] {
76
362880l, 3628800l, 39916800l,
77
479001600l, 6227020800l, 87178291200l,
78
1307674368000l, 20922789888000l, 355687428096000l,
79
6402373705728000l, 121645100408832000l, 2432902008176640000l };
69
82
* Private Constructor
197
210
// For n <= 61, the naive implementation cannot overflow.
198
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
212
for (int j = 1; j <= k; j++) {
199
213
result = result * i / j;
201
216
} else if (n <= 66) {
202
217
// For n > 61 but n <= 66, the result cannot overflow,
203
218
// but we must take care not to overflow intermediate values.
204
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
220
for (int j = 1; j <= k; j++) {
205
221
// We know that (result * i) is divisible by j,
206
222
// but (result * i) may overflow, so we split j:
207
223
// Filter out the gcd, d, so j/d and i/d are integer.
208
224
// result is divisible by (j/d) because (j/d)
209
225
// is relative prime to (i/d) and is a divisor of
210
226
// result * (i/d).
227
final long d = gcd(i, j);
212
228
result = (result / (j / d)) * (i / d);
215
232
// For n > 66, a result overflow might occur, so we check
216
233
// the multiplication, taking care to not overflow
218
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
220
result = mulAndCheck((result / (j / d)), (i / d));
236
for (int j = 1; j <= k; j++) {
237
final long d = gcd(i, j);
238
result = mulAndCheck(result / (j / d), i / d);
294
313
if ((k == 1) || (k == n - 1)) {
295
314
return Math.log(n);
299
318
* For values small enough to do exact integer computation,
300
* return the log of the exact value
319
* return the log of the exact value
303
322
return Math.log(binomialCoefficient(n,k));
307
326
* Return the log of binomialCoefficientDouble for values that will not
308
327
* overflow binomialCoefficientDouble
311
330
return Math.log(binomialCoefficientDouble(n, k));
315
334
return binomialCoefficientLog(n, n - k);
376
395
* Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
377
396
* hyperbolic cosine</a> of x.
379
398
* @param x double value for which to find the hyperbolic cosine
380
399
* @return hyperbolic cosine of x
382
401
public static double cosh(double x) {
383
402
return (Math.exp(x) + Math.exp(-x)) / 2.0;
387
406
* Returns true iff both arguments are NaN or neither is NaN and they are
390
409
* @param x first value
391
410
* @param y second value
392
411
* @return true if the values are equal or both are NaN
394
413
public static boolean equals(double x, double y) {
395
return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
414
return (Double.isNaN(x) && Double.isNaN(y)) || x == y;
472
/** All long-representable factorials */
473
private static final long[] factorials = new long[]
474
{1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
475
479001600, 6227020800l, 87178291200l, 1307674368000l, 20922789888000l,
476
355687428096000l, 6402373705728000l, 121645100408832000l,
477
2432902008176640000l};
480
492
* Returns n!. Shorthand for <code>n</code> <a
593
605
* <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
594
606
* <code>0</code>.</li>
597
609
* @param p any number
598
610
* @param q any number
599
611
* @return the greatest common divisor, never negative
600
* @throws ArithmeticException
601
* if the result cannot be represented as a nonnegative int
612
* @throws ArithmeticException if the result cannot be represented as a
613
* nonnegative int value
605
616
public static int gcd(final int p, final int q) {
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
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>
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
705
public static long gcd(final long p, final long 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",
714
return Math.abs(u) + Math.abs(v);
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
720
/* assert u!=0 && v!=0; */
727
// B1. [Find power of 2]
729
while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
733
k++; // cast out twos.
736
throw MathRuntimeException.createArithmeticException(
737
"overflow: gcd({0}, {1}) is 2^63",
740
// B2. Initialize: u and v have been divided by 2^k and at least
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)
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
751
// B5 [reset max(u,v)]
757
// B6/B3. at this point both u and v should be odd.
759
// |u| larger: t positive (replace u)
760
// |v| larger: t negative (replace v)
762
return -u * (1L << k); // gcd is u*2^k
666
766
* Returns an integer hash code representing the given double value.
668
768
* @param value the value to be hashed
669
769
* @return the hash code
784
884
int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b));
785
if (lcm == Integer.MIN_VALUE){
786
throw new ArithmeticException("overflow: lcm is 2^31");
885
if (lcm == Integer.MIN_VALUE) {
886
throw MathRuntimeException.createArithmeticException(
887
"overflow: lcm({0}, {1}) is 2^31",
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>.
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>.
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
915
public static long lcm(long a, long b) {
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",
793
930
* <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
794
931
* for base <code>b</code> of <code>x</code>.
796
* <p>Returns <code>NaN<code> if either argument is negative. If
933
* <p>Returns <code>NaN<code> if either argument is negative. If
797
934
* <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
798
* If <code>base</code> is positive and <code>x</code> is 0,
935
* If <code>base</code> is positive and <code>x</code> is 0,
799
936
* <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments
800
937
* are 0, the result is <code>NaN</code>.</p>
802
939
* @param base the base of the logarithm, must be greater than 0
803
940
* @param x argument, must be greater than 0
804
941
* @return the value of the logarithm - the number y such that base^y = x.
807
944
public static double log(double base, double x) {
808
945
return Math.log(x)/Math.log(base);
812
949
* Multiply two integers, checking for overflow.
814
951
* @param x a factor
815
952
* @param y a factor
816
953
* @return the product <code>x*y</code>
988
1125
return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI);
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
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>
1136
* <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite
1137
* or NaN and ArithmeticException if the input array contains any infinite elements
1140
* <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
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
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");
1155
if (Double.isNaN(normalizedSum)) {
1156
throw MathRuntimeException.createIllegalArgumentException(
1157
"Cannot normalize to NaN");
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);
1167
if (!Double.isNaN(values[i])) {
1172
throw MathRuntimeException.createArithmeticException(
1173
"Array sums to zero");
1175
for (int i = 0; i < len; i++) {
1176
if (Double.isNaN(values[i])) {
1177
out[i] = Double.NaN;
1179
out[i] = values[i] * normalizedSum / sum;
992
1186
* Round the given value to the specified number of decimal places. The
993
1187
* value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
995
1189
* @param x the value to round.
996
1190
* @param scale the number of digits to the right of the decimal point.
997
1191
* @return the rounded value.
1160
1354
* For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
1161
1355
* x = 0, and (byte)(-1) if x < 0.</p>
1163
1357
* @param x the value, a byte
1164
1358
* @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
1175
1369
* <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
1176
1370
* <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
1177
1371
* Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
1179
1373
* @param x the value, a double
1180
1374
* @return +1.0, 0.0, or -1.0, depending on the sign of x
1193
1387
* For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
1194
1388
* 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
1195
1389
* is <code>NaN</code>.</p>
1197
1391
* @param x the value, a float
1198
1392
* @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
1239
1433
* For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
1240
1434
* if x = 0, and (short)(-1) if x < 0.</p>
1242
1436
* @param x the value, a short
1243
1437
* @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
1509
1703
* @param p2 the second point
1510
1704
* @return the L<sub>1</sub> distance between the two points
1512
public static final double distance1(double[] p1, double[] p2) {
1706
public static double distance1(double[] p1, double[] p2) {
1513
1707
double sum = 0;
1514
1708
for (int i = 0; i < p1.length; i++) {
1515
1709
sum += Math.abs(p1[i] - p2[i]);
1521
1715
* Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1539
1733
* @param p2 the second point
1540
1734
* @return the L<sub>2</sub> distance between the two points
1542
public static final double distance(double[] p1, double[] p2) {
1736
public static double distance(double[] p1, double[] p2) {
1543
1737
double sum = 0;
1544
1738
for (int i = 0; i < p1.length; i++) {
1545
1739
final double dp = p1[i] - p2[i];
1555
1749
* @param p2 the second point
1556
1750
* @return the L<sub>2</sub> distance between the two points
1558
public static final double distance(int[] p1, int[] p2) {
1752
public static double distance(int[] p1, int[] p2) {
1560
1754
for (int i = 0; i < p1.length; i++) {
1561
final int dp = p1[i] - p2[i];
1755
final double dp = p1[i] - p2[i];
1562
1756
sum += dp * dp;
1564
1758
return Math.sqrt(sum);
1568
1762
* Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1571
1765
* @param p2 the second point
1572
1766
* @return the L<sub>∞</sub> distance between the two points
1574
public static final double distanceInf(double[] p1, double[] p2) {
1768
public static double distanceInf(double[] p1, double[] p2) {
1575
1769
double max = 0;
1576
1770
for (int i = 0; i < p1.length; i++) {
1577
1771
max = Math.max(max, Math.abs(p1[i] - p2[i]));
1583
1777
* Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1586
1780
* @param p2 the second point
1587
1781
* @return the L<sub>∞</sub> distance between the two points
1589
public static final int distanceInf(int[] p1, int[] p2) {
1783
public static int distanceInf(int[] p1, int[] p2) {
1591
1785
for (int i = 0; i < p1.length; i++) {
1592
1786
max = Math.max(max, Math.abs(p1[i] - p2[i]));
1792
* Checks that the given array is sorted.
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.
1799
public static void checkOrder(double[] val, int dir, boolean strict) {
1800
double previous = val[0];
1802
int max = val.length;
1803
for (int i = 1; i < max; i++) {
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]);
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]);
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]);
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]);
added
removed
src/main/java/org/apache/commons/math/util/MathUtils.java
