MathUtils


1   /*
2    * Copyright 2003-2005 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  import java.math.BigDecimal  ;
20  
21  /**
22   * Some useful additions to the built-in functions in {@link Math}.
23   *
24   * @version $Revision$ $Date: 2005-06-25 18:19:16 -0700 (Sat, 25 Jun 2005) $
25   */
26  public final class MathUtils {
27      
28      /** 0.0 cast as a byte. */
29      private static final byte ZB = (byte) 0;
30      
31      /** -1.0 cast as a byte. */
32      private static final byte NB = (byte) -1;
33      
34      /** 1.0 cast as a byte. */
35      private static final byte PB = (byte) 1;
36      
37      /** 0.0 cast as a short. */
38      private static final short ZS = (short) 0;
39      
40      /** -1.0 cast as a short. */
41      private static final short NS = (short) -1;
42      
43      /** 1.0 cast as a short. */
44      private static final short PS = (short) 1;
45      
46      /**
47       * Private Constructor
48       */
49      private MathUtils() {
50      }
51      
52      /**
53       * Round the given value to the specified number of decimal places.  The
54       * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
55       * @param x the value to round.
56       * @param scale the number of digits to the right of the decimal point.
57       * @return the rounded value.
58       * @since 1.1
59       */
60      public static double round(double x, int scale) {
61          return round(x, scale, BigDecimal.ROUND_HALF_UP);
62      }
63  
64      /**
65       * Round the given value to the specified number of decimal places.  The
66       * value is rounded using the given method which is any method defined in
67       * {@link BigDecimal}.
68       * @param x the value to round.
69       * @param scale the number of digits to the right of the decimal point.
70       * @param roundingMethod the rounding method as defined in
71       *        {@link BigDecimal}. 
72       * @return the rounded value.
73       * @since 1.1
74       */
75      public static double round(
76          double x, int scale, int roundingMethod)
77      {
78          return (new BigDecimal  (x).setScale(scale, roundingMethod))
79              .doubleValue();
80      }
81      
82      /**
83       * Round the given value to the specified number of decimal places.  The
84       * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
85       * @param x the value to round.
86       * @param scale the number of digits to the right of the decimal point.
87       * @return the rounded value.
88       * @since 1.1
89       */
90      public static float round(float x, int scale) {
91          return round(x, scale, BigDecimal.ROUND_HALF_UP);
92      }
93  
94      /**
95       * Round the given value to the specified number of decimal places.  The
96       * value is rounded using the given method which is any method defined in
97       * {@link BigDecimal}.
98       * @param x the value to round.
99       * @param scale the number of digits to the right of the decimal point.
100      * @param roundingMethod the rounding method as defined in
101      *        {@link BigDecimal}. 
102      * @return the rounded value.
103      * @since 1.1
104      */
105     public static float round(float x, int scale, int roundingMethod) {
106         return (new BigDecimal  (x).setScale(scale, roundingMethod)).floatValue();
107     }
108     
109     /**
110      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
111      * sign</a> for double precision <code>x</code>.
112      *
113      * <p>
114      * For a double value <code>x</code>, this method returns <code>+1.0</code>
115      * if <code>x > 0</code>, <code>0.0</code> if <code>x = 0.0</code>,
116      * and <code>-1.0</code> if <code>x < 0</code>.  Returns <code>NaN</code> 
117      * if <code>x</code> is <code>NaN</code>.
118      *
119      * @param x the value, a double
120      * @return +1.0, 0.0, or -1.0, depending on the sign of x
121      */
122     public static double sign(final double x) {
123         if (Double.isNaN(x)) {
124             return Double.NaN;
125         }
126         return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
127     }
128     
129     /**
130      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
131      * sign</a> for float value <code>x</code>.
132      *
133      * <p>
134      * For a float value x, this method returns +1.0F if x > 0, 0.0F if
135      * x = 0.0F, and -1.0F if x < 0.  Returns <code>NaN</code> 
136      * if <code>x</code> is <code>NaN</code>.
137      *
138      * @param x the value, a float
139      * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
140      */
141     public static float sign(final float x) {
142         if (Float.isNaN(x)) {
143             return Float.NaN;
144         }
145         return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
146     }
147     
148     /**
149      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
150      * sign</a> for byte value <code>x</code>.
151      *
152      * <p>
153      * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0)
154      * if x = 0, and (byte)(-1) if x < 0.
155      *
156      * @param x the value, a byte
157      * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
158      */
159     public static byte sign(final byte x) {
160         return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
161     }
162     
163     /**
164      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
165      * sign</a> for short value <code>x</code>.
166      *
167      * <p>
168      * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
169      * if x = 0, and (short)(-1) if x < 0.
170      *
171      * @param x the value, a short
172      * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign
173      * of x
174      */
175     public static short sign(final short x) {
176         return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
177     }
178     
179     /**
180      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
181      * sign</a> for int value <code>x</code>.
182      *
183      * <p>
184      * For an int value x, this method returns +1 if x > 0, 0 if x = 0,
185      * and -1 if x < 0.
186      *
187      * @param x the value, an int
188      * @return +1, 0, or -1, depending on the sign of x
189      */
190     public static int sign(final int x) {
191         return (x == 0) ? 0 : (x > 0) ? 1 : -1;
192     }
193     
194     /**
195      * Returns the <a HREF="http://mathworld.wolfram.com/Sign.html">
196      * sign</a> for long value <code>x</code>.
197      *
198      * <p>
199      * For a long value x, this method returns +1L if x > 0, 0L if x = 0,
200      * and -1L if x < 0.
201      *
202      * @param x the value, a long
203      * @return +1L, 0L, or -1L, depending on the sign of x
204      */
205     public static long sign(final long x) {
206         return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
207     }
208     
209     /**
210      * For a double precision value x, this method returns +1.0 if x >= 0
211      * and -1.0 if x < 0.   Returns <code>NaN</code> 
212      * if <code>x</code> is <code>NaN</code>.
213      *
214      * @param x the value, a double
215      * @return +1.0 or -1.0, depending on the sign of x
216      */
217     public static double indicator(final double x) {
218         if (Double.isNaN(x)) {
219             return Double.NaN;
220         }
221         return (x >= 0.0) ? 1.0 : -1.0;
222     }
223     
224     /**
225      * For a float value x, this method returns +1.0F if x >= 0
226      * and -1.0F if x < 0.   Returns <code>NaN</code> 
227      * if <code>x</code> is <code>NaN</code>.
228      *
229      * @param x the value, a float
230      * @return +1.0F or -1.0F, depending on the sign of x
231      */
232     public static float indicator(final float x) {
233         if (Float.isNaN(x)) {
234             return Float.NaN;
235         }
236         return (x >= 0.0F) ? 1.0F : -1.0F;
237     }
238     
239     /**
240      * For a byte value x, this method returns (byte)(+1) if x >= 0
241      * and (byte)(-1) if x < 0.
242      *
243      * @param x the value, a byte
244      * @return (byte)(+1) or (byte)(-1), depending on the sign of x
245      */
246     public static byte indicator(final byte x) {
247         return (x >= ZB) ? PB : NB;
248     }
249     
250     /**
251      * For a short value x, this method returns (short)(+1) if x >= 0
252      * and (short)(-1) if x < 0.
253      *
254      * @param x the value, a short
255      * @return (short)(+1) or (short)(-1), depending on the sign of x
256      */
257     public static short indicator(final short x) {
258         return (x >= ZS) ? PS : NS;
259     }
260     
261     /**
262      * For an int value x, this method returns +1 if x >= 0
263      * and -1 if x < 0.
264      *
265      * @param x the value, an int
266      * @return +1 or -1, depending on the sign of x
267      */
268     public static int indicator(final int x) {
269         return (x >= 0) ? 1 : -1;
270     }
271     
272     /**
273      * For a long value x, this method returns +1L if x >= 0
274      * and -1L if x < 0.
275      *
276      * @param x the value, a long
277      * @return +1L or -1L, depending on the sign of x
278      */
279     public static long indicator(final long x) {
280         return (x >= 0L) ? 1L : -1L;
281     }
282     
283     /**
284      * Returns an exact representation of the
285      * <a HREF="http://mathworld.wolfram.com/BinomialCoefficient.html">
286      * Binomial Coefficient</a>,  "<code>n choose k</code>",
287      * the number of <code>k</code>-element subsets that can be selected from
288      * an <code>n</code>-element set.
289      * <p>
290      * <Strong>Preconditions</strong>:<ul>
291      * <li> <code>0 <= k <= n </code> (otherwise
292      *      <code>IllegalArgumentException</code> is thrown)</li>
293      * <li> The result is small enough to fit into a <code>long</code>.  The
294      *      largest value of <code>n</code> for which all coefficients are
295      *      <code> < Long.MAX_VALUE</code> is 66.  If the computed value
296      *      exceeds <code>Long.MAX_VALUE</code> an <code>ArithMeticException
297      *      </code> is thrown.</li>
298      * </ul>
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      * @throws ArithmeticException if the result is too large to be represented
305      *         by a long integer.
306      */
307     public static long binomialCoefficient(final int n, final int k) {
308         if (n < k) {
309             throw new IllegalArgumentException  (
310             "must have n >= k for binomial coefficient (n,k)");
311         }
312         if (n < 0) {
313             throw new IllegalArgumentException  (
314             "must have n >= 0 for binomial coefficient (n,k)");
315         }
316         if ((n == k) || (k == 0)) {
317             return 1;
318         }
319         if ((k == 1) || (k == n - 1)) {
320             return n;
321         }
322         
323         long result = Math.round(binomialCoefficientDouble(n, k));
324         if (result == Long.MAX_VALUE) {
325             throw new ArithmeticException  (
326             "result too large to represent in a long integer");
327         }
328         return result;
329     }
330     
331     /**
332      * Returns a <code>double</code> representation of the
333      * <a HREF="http://mathworld.wolfram.com/BinomialCoefficient.html">
334      * Binomial Coefficient</a>,  "<code>n choose k</code>",
335      * the number of <code>k</code>-element subsets that can be selected from
336      * an <code>n</code>-element set.
337      * <p>
338      * <Strong>Preconditions</strong>:<ul>
339      * <li> <code>0 <= k <= n </code> (otherwise
340      *      <code>IllegalArgumentException</code> is thrown)</li>
341      * <li> The result is small enough to fit into a <code>double</code>.
342      *      The largest value of <code>n</code> for which all coefficients are
343      *      < Double.MAX_VALUE is 1029.  If the computed value exceeds
344      *      Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
345      * </ul>
346      *
347      * @param n the size of the set
348      * @param k the size of the subsets to be counted
349      * @return <code>n choose k</code>
350      * @throws IllegalArgumentException if preconditions are not met.
351      */
352     public static double binomialCoefficientDouble(final int n, final int k) {
353         return Math.floor(Math.exp(binomialCoefficientLog(n, k)) + 0.5);
354     }
355     
356     /**
357      * Returns the natural <code>log</code> of the
358      * <a HREF="http://mathworld.wolfram.com/BinomialCoefficient.html">
359      * Binomial Coefficient</a>,  "<code>n choose k</code>",
360      * the number of <code>k</code>-element subsets that can be selected from
361      * an <code>n</code>-element set.
362      * <p>
363      * <Strong>Preconditions</strong>:<ul>
364      * <li> <code>0 <= k <= n </code> (otherwise
365      *      <code>IllegalArgumentException</code> is thrown)</li>
366      * </ul>
367      *
368      * @param n the size of the set
369      * @param k the size of the subsets to be counted
370      * @return <code>n choose k</code>
371      * @throws IllegalArgumentException if preconditions are not met.
372      */
373     public static double binomialCoefficientLog(final int n, final int k) {
374         if (n < k) {
375             throw new IllegalArgumentException  (
376             "must have n >= k for binomial coefficient (n,k)");
377         }
378         if (n < 0) {
379             throw new IllegalArgumentException  (
380             "must have n >= 0 for binomial coefficient (n,k)");
381         }
382         if ((n == k) || (k == 0)) {
383             return 0;
384         }
385         if ((k == 1) || (k == n - 1)) {
386             return Math.log((double) n);
387         }
388         double logSum = 0;
389         
390         // n!/k!
391         for (int i = k + 1; i <= n; i++) {
392             logSum += Math.log((double) i);
393         }
394         
395         // divide by (n-k)!
396         for (int i = 2; i <= n - k; i++) {
397             logSum -= Math.log((double) i);
398         }
399         
400         return logSum;
401     }
402     
403     /**
404      * Returns n!.  Shorthand for <code>n</code>
405      * <a HREF="http://mathworld.wolfram.com/Factorial.html">
406      * Factorial</a>, the product of the numbers <code>1,...,n</code>.
407      *
408      * <p>
409      * <Strong>Preconditions</strong>:<ul>
410      * <li> <code>n >= 0</code> (otherwise
411      *      <code>IllegalArgumentException</code> is thrown)</li>
412      * <li> The result is small enough to fit into a <code>long</code>.  The
413      *      largest value of <code>n</code> for which <code>n!</code>
414      *      < Long.MAX_VALUE</code> is 20.  If the computed value
415      *      exceeds <code>Long.MAX_VALUE</code> an <code>ArithMeticException
416      *      </code> is thrown.</li>
417      * </ul>
418      * </p>
419      *
420      * @param n argument
421      * @return <code>n!</code>
422      * @throws ArithmeticException if the result is too large to be represented
423      *         by a long integer.
424      * @throws IllegalArgumentException if n < 0
425      */
426     public static long factorial(final int n) {
427         long result = Math.round(factorialDouble(n));
428         if (result == Long.MAX_VALUE) {
429             throw new ArithmeticException  (
430             "result too large to represent in a long integer");
431         }
432         return result;
433     }
434     
435     /**
436      * Returns n!.  Shorthand for <code>n</code>
437      * <a HREF="http://mathworld.wolfram.com/Factorial.html">
438      * Factorial</a>, the product of the numbers <code>1,...,n</code> as a
439      * <code>double</code>.
440      *
441      * <p>
442      * <Strong>Preconditions</strong>:<ul>
443      * <li> <code>n >= 0</code> (otherwise
444      *      <code>IllegalArgumentException</code> is thrown)</li>
445      * <li> The result is small enough to fit into a <code>double</code>.  The
446      *      largest value of <code>n</code> for which <code>n!</code>
447      *      < Double.MAX_VALUE</code> is 170.  If the computed value exceeds
448      *      Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
449      * </ul>
450      * </p>
451      *
452      * @param n argument
453      * @return <code>n!</code>
454      * @throws IllegalArgumentException if n < 0
455      */
456     public static double factorialDouble(final int n) {
457         if (n < 0) {
458             throw new IllegalArgumentException  ("must have n >= 0 for n!");
459         }
460         return Math.floor(Math.exp(factorialLog(n)) + 0.5);
461     }
462     
463     /**
464      * Returns the natural logarithm of n!.
465      * <p>
466      * <Strong>Preconditions</strong>:<ul>
467      * <li> <code>n >= 0</code> (otherwise
468      *      <code>IllegalArgumentException</code> is thrown)</li>
469      * </ul>
470      *
471      * @param n argument
472      * @return <code>n!</code>
473      * @throws IllegalArgumentException if preconditions are not met.
474      */
475     public static double factorialLog(final int n) {
476         if (n < 0) {
477             throw new IllegalArgumentException  ("must have n > 0 for n!");
478         }
479         double logSum = 0;
480         for (int i = 2; i <= n; i++) {
481             logSum += Math.log((double) i);
482         }
483         return logSum;
484     }
485     
486     /**
487      * Returns the <a HREF="http://mathworld.wolfram.com/HyperbolicCosine.html">
488      * hyperbolic cosine</a> of x.
489      *
490      * @param x double value for which to find the hyperbolic cosine
491      * @return hyperbolic cosine of x
492      */
493     public static double cosh(double x) {
494         return (Math.exp(x) + Math.exp(-x)) / 2.0;
495     }
496     
497     /**
498      * Returns the <a HREF="http://mathworld.wolfram.com/HyperbolicSine.html">
499      * hyperbolic sine</a> of x.
500      *
501      * @param x double value for which to find the hyperbolic sine
502      * @return hyperbolic sine of x
503      */
504     public static double sinh(double x) {
505         return (Math.exp(x) - Math.exp(-x)) / 2.0;
506     }
507     
508     /**
509      * Returns an integer hash code representing the given double value.
510      *
511      * @param value  the value to be hashed
512      * @return the hash code
513      */
514     public static int hash(double value) {
515         long bits = Double.doubleToLongBits(value);
516         return (int)(bits ^ (bits >>> 32));
517     }
518     
519     /**
520      * Returns true iff both arguments are NaN or
521      * neither is NaN and they are equal
522      *
523      * @param x first value
524      * @param y second value
525      * @return true if the values are equal or both are NaN
526      */
527     public static boolean equals(double x, double y) {
528         return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
529     }
530 
531     /**
532      * Returns the least common multiple between two integer values.
533      * 
534      * @param a the first integer value.
535      * @param b the second integer value.
536      * @return the least common multiple between a and b.
537      * @throws ArithmeticException if the lcm is too large to store as an int
538      * @since 1.1
539      */
540     public static int lcm(int a, int b) {
541         return Math.abs(mulAndCheck(a / gcd(a, b) , b));
542     }
543 
544     /**
545      * <p>Gets the greatest common divisor of the absolute value of
546      * two numbers, using the "binary gcd" method which avoids
547      * division and modulo operations.  See Knuth 4.5.2 algorithm B.
548      * This algorithm is due to Josef Stein (1961).</p>
549      *
550      * @param u  a non-zero number
551      * @param v  a non-zero number
552      * @return the greatest common divisor, never zero
553      * @since 1.1
554      */
555     public static int gcd(int u, int v) {
556         if (u * v == 0) {
557             return (Math.abs(u) + Math.abs(v));
558         }
559         // keep u and v negative, as negative integers range down to
560         // -2^31, while positive numbers can only be as large as 2^31-1
561         // (i.e. we can't necessarily negate a negative number without
562         // overflow)
563         /* assert u!=0 && v!=0; */
564         if (u>0) { u=-u; } // make u negative
565         if (v>0) { v=-v; } // make v negative
566         // B1. [Find power of 2]
567         int k=0;
568         while ((u&1)==0 && (v&1)==0 && k<31) { // while u and v are both even...
569             u/=2; v/=2; k++; // cast out twos.
570         }
571         if (k==31) {
572             throw new ArithmeticException  ("overflow: gcd is 2^31");
573         }
574         // B2. Initialize: u and v have been divided by 2^k and at least
575         //     one is odd.
576         int t = ((u&1)==1) ? v : -(u/2)/*B3*/;
577         // t negative: u was odd, v may be even (t replaces v)
578         // t positive: u was even, v is odd (t replaces u)
579         do {
580             /* assert u<0 && v<0; */
581             // B4/B3: cast out twos from t.
582             while ((t&1)==0) { // while t is even..
583                 t/=2; // cast out twos
584             }
585             // B5 [reset max(u,v)]
586             if (t>0) {
587                 u = -t;
588             } else {
589                 v = t;
590             }
591             // B6/B3. at this point both u and v should be odd.
592             t = (v - u)/2;
593             // |u| larger: t positive (replace u)
594             // |v| larger: t negative (replace v)
595         } while (t!=0);
596         return -u*(1<<k); // gcd is u*2^k
597     }
598 
599     /** 
600      * Multiply two integers, checking for overflow.
601      * 
602      * @param x a factor
603      * @param y a factor
604      * @return the product <code>x*y</code>
605      * @throws ArithmeticException if the result can not be represented as
606      *                             an int
607      * @since 1.1
608      */
609     public static int mulAndCheck(int x, int y) {
610         long m = ((long)x)*((long)y);
611         if (m < Integer.MIN_VALUE ||
612                 m > Integer.MAX_VALUE) {
613             throw new ArithmeticException  ("overflow: mul");
614         }
615         return (int)m;
616     }
617     
618     /** 
619      * Add two integers, checking for overflow.
620      * 
621      * @param x an addend
622      * @param y an addend
623      * @return the sum <code>x+y</code>
624      * @throws ArithmeticException if the result can not be represented as
625      * an int
626      * @since 1.1
627      */
628     public static int addAndCheck(int x, int y) {
629         long s = (long)x+(long)y;
630         if (s < Integer.MIN_VALUE ||
631                 s > Integer.MAX_VALUE) {
632             throw new ArithmeticException  ("overflow: add");
633         }
634         return (int)s;
635     }
636     
637     /** 
638      * Subtract two integers, checking for overflow.
639      * 
640      * @param x the minuend
641      * @param y the subtrahend
642      * @return the difference <code>x-y</code>
643      * @throws ArithmeticException if the result can not be represented as
644      * an int
645      * @since 1.1
646      */
647     public static int subAndCheck(int x, int y) {
648         long s = (long)x-(long)y;
649         if (s < Integer.MIN_VALUE ||
650                 s > Integer.MAX_VALUE) {
651             throw new ArithmeticException  ("overflow: add");
652         }
653         return (int)s;
654     }
655 }
656
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags