1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.util;
18
19 /***
20 * Some useful additions to the built-in functions in {@link Math}.
21 *
22 * @version $Revision: 1.20 $ $Date: 2004/10/14 04:01:04 $
23 */
24 public final class MathUtils {
25
26 /*** 0.0 cast as a byte. */
27 private static final byte ZB = (byte) 0;
28
29 /*** -1.0 cast as a byte. */
30 private static final byte NB = (byte) -1;
31
32 /*** 1.0 cast as a byte. */
33 private static final byte PB = (byte) 1;
34
35 /*** 0.0 cast as a short. */
36 private static final short ZS = (short) 0;
37
38 /*** -1.0 cast as a short. */
39 private static final short NS = (short) -1;
40
41 /*** 1.0 cast as a short. */
42 private static final short PS = (short) 1;
43
44 /***
45 * Private Constructor
46 */
47 private MathUtils() {
48 }
49
50 /***
51 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
52 * sign</a> for double precision <code>x</code>.
53 *
54 * <p>
55 * For a double value <code>x</code>, this method returns <code>+1.0</code>
56 * if <code>x > 0</code>, <code>0.0</code> if <code>x = 0.0</code>,
57 * and <code>-1.0</code> if <code>x < 0</code>. Returns <code>NaN</code>
58 * if <code>x</code> is <code>NaN</code>.
59 *
60 * @param x the value, a double
61 * @return +1.0, 0.0, or -1.0, depending on the sign of x
62 */
63 public static double sign(final double x) {
64 if (Double.isNaN(x)) {
65 return Double.NaN;
66 }
67 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
68 }
69
70 /***
71 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
72 * sign</a> for float value <code>x</code>.
73 *
74 * <p>
75 * For a float value x, this method returns +1.0F if x > 0, 0.0F if
76 * x = 0.0F, and -1.0F if x < 0. Returns <code>NaN</code>
77 * if <code>x</code> is <code>NaN</code>.
78 *
79 * @param x the value, a float
80 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
81 */
82 public static float sign(final float x) {
83 if (Float.isNaN(x)) {
84 return Float.NaN;
85 }
86 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
87 }
88
89 /***
90 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
91 * sign</a> for byte value <code>x</code>.
92 *
93 * <p>
94 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0)
95 * if x = 0, and (byte)(-1) if x < 0.
96 *
97 * @param x the value, a byte
98 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
99 */
100 public static byte sign(final byte x) {
101 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
102 }
103
104 /***
105 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
106 * sign</a> for short value <code>x</code>.
107 *
108 * <p>
109 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
110 * if x = 0, and (short)(-1) if x < 0.
111 *
112 * @param x the value, a short
113 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign
114 * of x
115 */
116 public static short sign(final short x) {
117 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
118 }
119
120 /***
121 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
122 * sign</a> for int value <code>x</code>.
123 *
124 * <p>
125 * For an int value x, this method returns +1 if x > 0, 0 if x = 0,
126 * and -1 if x < 0.
127 *
128 * @param x the value, an int
129 * @return +1, 0, or -1, depending on the sign of x
130 */
131 public static int sign(final int x) {
132 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
133 }
134
135 /***
136 * Returns the <a href="http://mathworld.wolfram.com/Sign.html">
137 * sign</a> for long value <code>x</code>.
138 *
139 * <p>
140 * For a long value x, this method returns +1L if x > 0, 0L if x = 0,
141 * and -1L if x < 0.
142 *
143 * @param x the value, a long
144 * @return +1L, 0L, or -1L, depending on the sign of x
145 */
146 public static long sign(final long x) {
147 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
148 }
149
150 /***
151 * For a double precision value x, this method returns +1.0 if x >= 0
152 * and -1.0 if x < 0. Returns <code>NaN</code>
153 * if <code>x</code> is <code>NaN</code>.
154 *
155 * @param x the value, a double
156 * @return +1.0 or -1.0, depending on the sign of x
157 */
158 public static double indicator(final double x) {
159 if (Double.isNaN(x)) {
160 return Double.NaN;
161 }
162 return (x >= 0.0) ? 1.0 : -1.0;
163 }
164
165 /***
166 * For a float value x, this method returns +1.0F if x >= 0
167 * and -1.0F if x < 0. Returns <code>NaN</code>
168 * if <code>x</code> is <code>NaN</code>.
169 *
170 * @param x the value, a float
171 * @return +1.0F or -1.0F, depending on the sign of x
172 */
173 public static float indicator(final float x) {
174 if (Float.isNaN(x)) {
175 return Float.NaN;
176 }
177 return (x >= 0.0F) ? 1.0F : -1.0F;
178 }
179
180 /***
181 * For a byte value x, this method returns (byte)(+1) if x >= 0
182 * and (byte)(-1) if x < 0.
183 *
184 * @param x the value, a byte
185 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
186 */
187 public static byte indicator(final byte x) {
188 return (x >= ZB) ? PB : NB;
189 }
190
191 /***
192 * For a short value x, this method returns (short)(+1) if x >= 0
193 * and (short)(-1) if x < 0.
194 *
195 * @param x the value, a short
196 * @return (short)(+1) or (short)(-1), depending on the sign of x
197 */
198 public static short indicator(final short x) {
199 return (x >= ZS) ? PS : NS;
200 }
201
202 /***
203 * For an int value x, this method returns +1 if x >= 0
204 * and -1 if x < 0.
205 *
206 * @param x the value, an int
207 * @return +1 or -1, depending on the sign of x
208 */
209 public static int indicator(final int x) {
210 return (x >= 0) ? 1 : -1;
211 }
212
213 /***
214 * For a long value x, this method returns +1L if x >= 0
215 * and -1L if x < 0.
216 *
217 * @param x the value, a long
218 * @return +1L or -1L, depending on the sign of x
219 */
220 public static long indicator(final long x) {
221 return (x >= 0L) ? 1L : -1L;
222 }
223
224 /***
225 * Returns an exact representation of the
226 * <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
227 * Binomial Coefficient</a>, "<code>n choose k</code>",
228 * the number of <code>k</code>-element subsets that can be selected from
229 * an <code>n</code>-element set.
230 * <p>
231 * <Strong>Preconditions</strong>:<ul>
232 * <li> <code>0 <= k <= n </code> (otherwise
233 * <code>IllegalArgumentException</code> is thrown)</li>
234 * <li> The result is small enough to fit into a <code>long</code>. The
235 * largest value of <code>n</code> for which all coefficients are
236 * <code> < Long.MAX_VALUE</code> is 66. If the computed value
237 * exceeds <code>Long.MAX_VALUE</code> an <code>ArithMeticException
238 * </code> is thrown.</li>
239 * </ul>
240 *
241 * @param n the size of the set
242 * @param k the size of the subsets to be counted
243 * @return <code>n choose k</code>
244 * @throws IllegalArgumentException if preconditions are not met.
245 * @throws ArithmeticException if the result is too large to be represented
246 * by a long integer.
247 */
248 public static long binomialCoefficient(final int n, final int k) {
249 if (n < k) {
250 throw new IllegalArgumentException(
251 "must have n >= k for binomial coefficient (n,k)");
252 }
253 if (n < 0) {
254 throw new IllegalArgumentException(
255 "must have n >= 0 for binomial coefficient (n,k)");
256 }
257 if ((n == k) || (k == 0)) {
258 return 1;
259 }
260 if ((k == 1) || (k == n - 1)) {
261 return n;
262 }
263
264 long result = Math.round(binomialCoefficientDouble(n, k));
265 if (result == Long.MAX_VALUE) {
266 throw new ArithmeticException(
267 "result too large to represent in a long integer");
268 }
269 return result;
270 }
271
272 /***
273 * Returns a <code>double</code> representation of the
274 * <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
275 * Binomial Coefficient</a>, "<code>n choose k</code>",
276 * the number of <code>k</code>-element subsets that can be selected from
277 * an <code>n</code>-element set.
278 * <p>
279 * <Strong>Preconditions</strong>:<ul>
280 * <li> <code>0 <= k <= n </code> (otherwise
281 * <code>IllegalArgumentException</code> is thrown)</li>
282 * <li> The result is small enough to fit into a <code>double</code>.
283 * The largest value of <code>n</code> for which all coefficients are
284 * < Double.MAX_VALUE is 1029. If the computed value exceeds
285 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
286 * </ul>
287 *
288 * @param n the size of the set
289 * @param k the size of the subsets to be counted
290 * @return <code>n choose k</code>
291 * @throws IllegalArgumentException if preconditions are not met.
292 */
293 public static double binomialCoefficientDouble(final int n, final int k) {
294 return Math.floor(Math.exp(binomialCoefficientLog(n, k)) + 0.5);
295 }
296
297 /***
298 * Returns the natural <code>log</code> of the
299 * <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
300 * Binomial Coefficient</a>, "<code>n choose k</code>",
301 * the number of <code>k</code>-element subsets that can be selected from
302 * an <code>n</code>-element set.
303 * <p>
304 * <Strong>Preconditions</strong>:<ul>
305 * <li> <code>0 <= k <= n </code> (otherwise
306 * <code>IllegalArgumentException</code> is thrown)</li>
307 * </ul>
308 *
309 * @param n the size of the set
310 * @param k the size of the subsets to be counted
311 * @return <code>n choose k</code>
312 * @throws IllegalArgumentException if preconditions are not met.
313 */
314 public static double binomialCoefficientLog(final int n, final int k) {
315 if (n < k) {
316 throw new IllegalArgumentException(
317 "must have n >= k for binomial coefficient (n,k)");
318 }
319 if (n < 0) {
320 throw new IllegalArgumentException(
321 "must have n >= 0 for binomial coefficient (n,k)");
322 }
323 if ((n == k) || (k == 0)) {
324 return 0;
325 }
326 if ((k == 1) || (k == n - 1)) {
327 return Math.log((double) n);
328 }
329 double logSum = 0;
330
331
332 for (int i = k + 1; i <= n; i++) {
333 logSum += Math.log((double) i);
334 }
335
336
337 for (int i = 2; i <= n - k; i++) {
338 logSum -= Math.log((double) i);
339 }
340
341 return logSum;
342 }
343
344 /***
345 * Returns n!. Shorthand for <code>n</code>
346 * <a href="http://mathworld.wolfram.com/Factorial.html">
347 * Factorial</a>, the product of the numbers <code>1,...,n</code>.
348 *
349 * <p>
350 * <Strong>Preconditions</strong>:<ul>
351 * <li> <code>n >= 0</code> (otherwise
352 * <code>IllegalArgumentException</code> is thrown)</li>
353 * <li> The result is small enough to fit into a <code>long</code>. The
354 * largest value of <code>n</code> for which <code>n!</code>
355 * < Long.MAX_VALUE</code> is 20. If the computed value
356 * exceeds <code>Long.MAX_VALUE</code> an <code>ArithMeticException
357 * </code> is thrown.</li>
358 * </ul>
359 * </p>
360 *
361 * @param n argument
362 * @return <code>n!</code>
363 * @throws ArithmeticException if the result is too large to be represented
364 * by a long integer.
365 * @throws IllegalArgumentException if n < 0
366 */
367 public static long factorial(final int n) {
368 long result = Math.round(factorialDouble(n));
369 if (result == Long.MAX_VALUE) {
370 throw new ArithmeticException(
371 "result too large to represent in a long integer");
372 }
373 return result;
374 }
375
376 /***
377 * Returns n!. Shorthand for <code>n</code>
378 * <a href="http://mathworld.wolfram.com/Factorial.html">
379 * Factorial</a>, the product of the numbers <code>1,...,n</code> as a
380 * <code>double</code>.
381 *
382 * <p>
383 * <Strong>Preconditions</strong>:<ul>
384 * <li> <code>n >= 0</code> (otherwise
385 * <code>IllegalArgumentException</code> is thrown)</li>
386 * <li> The result is small enough to fit into a <code>double</code>. The
387 * largest value of <code>n</code> for which <code>n!</code>
388 * < Double.MAX_VALUE</code> is 170. If the computed value exceeds
389 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
390 * </ul>
391 * </p>
392 *
393 * @param n argument
394 * @return <code>n!</code>
395 * @throws IllegalArgumentException if n < 0
396 */
397 public static double factorialDouble(final int n) {
398 if (n < 0) {
399 throw new IllegalArgumentException("must have n >= 0 for n!");
400 }
401 return Math.floor(Math.exp(factorialLog(n)) + 0.5);
402 }
403
404 /***
405 * Returns the natural logarithm of n!.
406 * <p>
407 * <Strong>Preconditions</strong>:<ul>
408 * <li> <code>n >= 0</code> (otherwise
409 * <code>IllegalArgumentException</code> is thrown)</li>
410 * </ul>
411 *
412 * @param n argument
413 * @return <code>n!</code>
414 * @throws IllegalArgumentException if preconditions are not met.
415 */
416 public static double factorialLog(final int n) {
417 if (n < 0) {
418 throw new IllegalArgumentException("must have n > 0 for n!");
419 }
420 double logSum = 0;
421 for (int i = 2; i <= n; i++) {
422 logSum += Math.log((double) i);
423 }
424 return logSum;
425 }
426
427 /***
428 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
429 * hyperbolic cosine</a> of x.
430 *
431 * @param x double value for which to find the hyperbolic cosine
432 * @return hyperbolic cosine of x
433 */
434 public static double cosh(double x) {
435 return (Math.exp(x) + Math.exp(-x)) / 2.0;
436 }
437
438 /***
439 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
440 * hyperbolic sine</a> of x.
441 *
442 * @param x double value for which to find the hyperbolic sine
443 * @return hyperbolic sine of x
444 */
445 public static double sinh(double x) {
446 return (Math.exp(x) - Math.exp(-x)) / 2.0;
447 }
448
449 /***
450 * Returns an integer hash code representing the given double value.
451 *
452 * @param value the value to be hashed
453 * @return the hash code
454 */
455 public static int hash(double value) {
456 long bits = Double.doubleToLongBits(value);
457 return (int)(bits ^ (bits >>> 32));
458 }
459
460 /***
461 * Returns true iff both arguments are NaN or
462 * neither is NaN and they are equal
463 *
464 * @param x first value
465 * @param y second value
466 * @return true if the values are equal or both are NaN
467 */
468 public static boolean equals(double x, double y) {
469 return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
470 }
471 }