View Javadoc

1   /*
2    * Copyright 2003-2004 The Apache Software Foundation.
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *      http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
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         // n!/k!
332         for (int i = k + 1; i <= n; i++) {
333             logSum += Math.log((double) i);
334         }
335         
336         // divide by (n-k)!
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 }