Save This Page
Home » mahout-collections-1.0-src » org.apache.mahout.collections » [javadoc | source]
    1   /**
    2    * Licensed to the Apache Software Foundation (ASF) under one
    3    * or more contributor license agreements. See the NOTICE file
    4    * distributed with this work for additional information
    5    * regarding copyright ownership. The ASF licenses this file
    6    * to you under the Apache License, Version 2.0 (the
    7    * "License"); you may not use this file except in compliance
    8    * with the License. You may obtain a copy of the License at
    9    *
   10    * http://www.apache.org/licenses/LICENSE-2.0
   11    *
   12    * Unless required by applicable law or agreed to in writing,
   13    * software distributed under the License is distributed on an
   14    * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
   15    * KIND, either express or implied. See the License for the
   16    * specific language governing permissions and limitations
   17    * under the License.
   18    */
   19   /*
   20   Copyright 1999 CERN - European Organization for Nuclear Research.
   21   Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
   22   is hereby granted without fee, provided that the above copyright notice appear in all copies and 
   23   that both that copyright notice and this permission notice appear in supporting documentation. 
   24   CERN makes no representations about the suitability of this software for any purpose. 
   25   It is provided "as is" without expressed or implied warranty.
   26   */
   27   package org.apache.mahout.collections;
   28   
   29   /**
   30    * Arithmetic functions.
   31    */
   32   public class Arithmetic extends Constants {
   33     // for method stirlingCorrection(...)
   34     private static final double[] stirlingCorrection = {
   35         0.0,
   36         8.106146679532726e-02, 4.134069595540929e-02,
   37         2.767792568499834e-02, 2.079067210376509e-02,
   38         1.664469118982119e-02, 1.387612882307075e-02,
   39         1.189670994589177e-02, 1.041126526197209e-02,
   40         9.255462182712733e-03, 8.330563433362871e-03,
   41         7.573675487951841e-03, 6.942840107209530e-03,
   42         6.408994188004207e-03, 5.951370112758848e-03,
   43         5.554733551962801e-03, 5.207655919609640e-03,
   44         4.901395948434738e-03, 4.629153749334029e-03,
   45         4.385560249232324e-03, 4.166319691996922e-03,
   46         3.967954218640860e-03, 3.787618068444430e-03,
   47         3.622960224683090e-03, 3.472021382978770e-03,
   48         3.333155636728090e-03, 3.204970228055040e-03,
   49         3.086278682608780e-03, 2.976063983550410e-03,
   50         2.873449362352470e-03, 2.777674929752690e-03,
   51     };
   52   
   53     // for method logFactorial(...)
   54     // log(k!) for k = 0, ..., 29
   55     private static final double[] logFactorials = {
   56         0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
   57         1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
   58         6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
   59         12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
   60         19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
   61         27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
   62         36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
   63         45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
   64         54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
   65         64.55753862700633106, 67.88974313718153498, 71.25703896716800901
   66     };
   67   
   68     // k! for k = 0, ..., 20
   69     private static final long[] longFactorials = {
   70         1L,
   71         1L,
   72         2L,
   73         6L,
   74         24L,
   75         120L,
   76         720L,
   77         5040L,
   78         40320L,
   79         362880L,
   80         3628800L,
   81         39916800L,
   82         479001600L,
   83         6227020800L,
   84         87178291200L,
   85         1307674368000L,
   86         20922789888000L,
   87         355687428096000L,
   88         6402373705728000L,
   89         121645100408832000L,
   90         2432902008176640000L
   91     };
   92   
   93     // k! for k = 21, ..., 170
   94     private static final double[] doubleFactorials = {
   95         5.109094217170944E19,
   96         1.1240007277776077E21,
   97         2.585201673888498E22,
   98         6.204484017332394E23,
   99         1.5511210043330984E25,
  100         4.032914611266057E26,
  101         1.0888869450418352E28,
  102         3.048883446117138E29,
  103         8.841761993739701E30,
  104         2.652528598121911E32,
  105         8.222838654177924E33,
  106         2.6313083693369355E35,
  107         8.68331761881189E36,
  108         2.952327990396041E38,
  109         1.0333147966386144E40,
  110         3.719933267899013E41,
  111         1.3763753091226346E43,
  112         5.23022617466601E44,
  113         2.0397882081197447E46,
  114         8.15915283247898E47,
  115         3.34525266131638E49,
  116         1.4050061177528801E51,
  117         6.041526306337384E52,
  118         2.6582715747884495E54,
  119         1.196222208654802E56,
  120         5.502622159812089E57,
  121         2.5862324151116827E59,
  122         1.2413915592536068E61,
  123         6.082818640342679E62,
  124         3.0414093201713376E64,
  125         1.5511187532873816E66,
  126         8.06581751709439E67,
  127         4.274883284060024E69,
  128         2.308436973392413E71,
  129         1.2696403353658264E73,
  130         7.109985878048632E74,
  131         4.052691950487723E76,
  132         2.350561331282879E78,
  133         1.386831185456898E80,
  134         8.32098711274139E81,
  135         5.075802138772246E83,
  136         3.146997326038794E85,
  137         1.9826083154044396E87,
  138         1.2688693218588414E89,
  139         8.247650592082472E90,
  140         5.443449390774432E92,
  141         3.6471110918188705E94,
  142         2.48003554243683E96,
  143         1.7112245242814127E98,
  144         1.1978571669969892E100,
  145         8.504785885678624E101,
  146         6.123445837688612E103,
  147         4.470115461512686E105,
  148         3.307885441519387E107,
  149         2.4809140811395404E109,
  150         1.8854947016660506E111,
  151         1.451830920282859E113,
  152         1.1324281178206295E115,
  153         8.94618213078298E116,
  154         7.15694570462638E118,
  155         5.797126020747369E120,
  156         4.7536433370128435E122,
  157         3.94552396972066E124,
  158         3.314240134565354E126,
  159         2.8171041143805494E128,
  160         2.4227095383672744E130,
  161         2.107757298379527E132,
  162         1.854826422573984E134,
  163         1.6507955160908465E136,
  164         1.4857159644817605E138,
  165         1.3520015276784033E140,
  166         1.2438414054641305E142,
  167         1.156772507081641E144,
  168         1.0873661566567426E146,
  169         1.0329978488239061E148,
  170         9.916779348709491E149,
  171         9.619275968248216E151,
  172         9.426890448883248E153,
  173         9.332621544394415E155,
  174         9.332621544394418E157,
  175         9.42594775983836E159,
  176         9.614466715035125E161,
  177         9.902900716486178E163,
  178         1.0299016745145631E166,
  179         1.0813967582402912E168,
  180         1.1462805637347086E170,
  181         1.2265202031961373E172,
  182         1.324641819451829E174,
  183         1.4438595832024942E176,
  184         1.5882455415227423E178,
  185         1.7629525510902457E180,
  186         1.974506857221075E182,
  187         2.2311927486598138E184,
  188         2.543559733472186E186,
  189         2.925093693493014E188,
  190         3.393108684451899E190,
  191         3.96993716080872E192,
  192         4.6845258497542896E194,
  193         5.574585761207606E196,
  194         6.689502913449135E198,
  195         8.094298525273444E200,
  196         9.875044200833601E202,
  197         1.2146304367025332E205,
  198         1.506141741511141E207,
  199         1.882677176888926E209,
  200         2.3721732428800483E211,
  201         3.0126600184576624E213,
  202         3.856204823625808E215,
  203         4.974504222477287E217,
  204         6.466855489220473E219,
  205         8.471580690878813E221,
  206         1.1182486511960037E224,
  207         1.4872707060906847E226,
  208         1.99294274616152E228,
  209         2.690472707318049E230,
  210         3.6590428819525483E232,
  211         5.0128887482749884E234,
  212         6.917786472619482E236,
  213         9.615723196941089E238,
  214         1.3462012475717523E241,
  215         1.8981437590761713E243,
  216         2.6953641378881633E245,
  217         3.8543707171800694E247,
  218         5.550293832739308E249,
  219         8.047926057471989E251,
  220         1.1749972043909107E254,
  221         1.72724589045464E256,
  222         2.5563239178728637E258,
  223         3.8089226376305687E260,
  224         5.7133839564458575E262,
  225         8.627209774233244E264,
  226         1.3113358856834527E267,
  227         2.0063439050956838E269,
  228         3.0897696138473515E271,
  229         4.789142901463393E273,
  230         7.471062926282892E275,
  231         1.1729568794264134E278,
  232         1.8532718694937346E280,
  233         2.946702272495036E282,
  234         4.714723635992061E284,
  235         7.590705053947223E286,
  236         1.2296942187394494E289,
  237         2.0044015765453032E291,
  238         3.287218585534299E293,
  239         5.423910666131583E295,
  240         9.003691705778434E297,
  241         1.5036165148649983E300,
  242         2.5260757449731988E302,
  243         4.2690680090047056E304,
  244         7.257415615308004E306
  245     };
  246   
  247     /** Makes this class non instantiable, but still let's others inherit from it. */
  248     protected Arithmetic() {
  249     }
  250   
  251     /**
  252      * Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k". The binomial
  253      * coefficient is defined as <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )</tt>. <ul> <li>k<0<tt>: <tt>0</tt>.
  254      * <li>k==0<tt>: <tt>1</tt>. <li>k==1<tt>: <tt>n</tt>. <li>else: <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k
  255      * )</tt>. </ul>
  256      *
  257      * @return the binomial coefficient.
  258      */
  259     public static double binomial(double n, long k) {
  260       if (k < 0) {
  261         return 0;
  262       }
  263       if (k == 0) {
  264         return 1;
  265       }
  266       if (k == 1) {
  267         return n;
  268       }
  269   
  270       // binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
  271       double a = n - k + 1;
  272       double b = 1;
  273       double binomial = 1;
  274       for (long i = k; i-- > 0;) {
  275         binomial *= (a++) / (b++);
  276       }
  277       return binomial;
  278     }
  279   
  280     /**
  281      * Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k". The binomial
  282      * coefficient is defined as <ul> <li>k<0<tt>: <tt>0</tt>. <li>k==0 || k==n<tt>: <tt>1</tt>. <li>k==1 || k==n-1<tt>:
  283      * <tt>n</tt>. <li>else: <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )</tt>. </ul>
  284      *
  285      * @return the binomial coefficient.
  286      */
  287     public static double binomial(long n, long k) {
  288       if (k < 0) {
  289         return 0;
  290       }
  291       if (k == 0 || k == n) {
  292         return 1;
  293       }
  294       if (k == 1 || k == n - 1) {
  295         return n;
  296       }
  297   
  298       // try quick version and see whether we get numeric overflows.
  299       // factorial(..) is O(1); requires no loop; only a table lookup.
  300       if (n > k) {
  301         int max = longFactorials.length + doubleFactorials.length;
  302         if (n < max) { // if (n! < inf && k! < inf)
  303           double n_fac = factorial((int) n);
  304           double k_fac = factorial((int) k);
  305           double n_minus_k_fac = factorial((int) (n - k));
  306           double nk = n_minus_k_fac * k_fac;
  307           if (nk != Double.POSITIVE_INFINITY) { // no numeric overflow?
  308             // now this is completely safe and accurate
  309             return n_fac / nk;
  310           }
  311         }
  312         if (k > n / 2) {
  313           k = n - k;
  314         } // quicker
  315       }
  316   
  317       // binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
  318       long a = n - k + 1;
  319       long b = 1;
  320       double binomial = 1;
  321       for (long i = k; i-- > 0;) {
  322         binomial *= ((double) (a++)) / (b++);
  323       }
  324       return binomial;
  325     }
  326   
  327     /**
  328      * Returns the smallest <code>long &gt;= value</code>. <dt>Examples: <code>1.0 -> 1, 1.2 -> 2, 1.9 -> 2</code>. This
  329      * method is safer than using (long) Math.ceil(value), because of possible rounding error.
  330      */
  331     public static long ceil(double value) {
  332       return Math.round(Math.ceil(value));
  333     }
  334   
  335     /**
  336      * Evaluates the series of Chebyshev polynomials Ti at argument x/2. The series is given by
  337      * <pre>
  338      *        N-1
  339      *         - '
  340      *  y  =   >   coef[i] T (x/2)
  341      *         -            i
  342      *        i=0
  343      * </pre>
  344      * Coefficients are stored in reverse order, i.e. the zero order term is last in the array.  Note N is the number of
  345      * coefficients, not the order. <p> If coefficients are for the interval a to b, x must have been transformed to x ->
  346      * 2(2x - b - a)/(b-a) before entering the routine.  This maps x from (a, b) to (-1, 1), over which the Chebyshev
  347      * polynomials are defined. <p> If the coefficients are for the inverted interval, in which (a, b) is mapped to (1/b,
  348      * 1/a), the transformation required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity, this becomes x -> 4a/x - 1.
  349      * <p> SPEED: <p> Taking advantage of the recurrence properties of the Chebyshev polynomials, the routine requires one
  350      * more addition per loop than evaluating a nested polynomial of the same degree.
  351      *
  352      * @param x    argument to the polynomial.
  353      * @param coef the coefficients of the polynomial.
  354      * @param N    the number of coefficients.
  355      */
  356     public static double chbevl(double x, double[] coef, int N) throws ArithmeticException {
  357   
  358       int p = 0;
  359   
  360       double b0 = coef[p++];
  361       double b1 = 0.0;
  362       int i = N - 1;
  363   
  364       double b2;
  365       do {
  366         b2 = b1;
  367         b1 = b0;
  368         b0 = x * b1 - b2 + coef[p++];
  369       } while (--i > 0);
  370   
  371       return (0.5 * (b0 - b2));
  372     }
  373   
  374     /**
  375      * Instantly returns the factorial <tt>k!</tt>.
  376      *
  377      * @param k must hold <tt>k &gt;= 0</tt>.
  378      */
  379     private static double factorial(int k) {
  380       if (k < 0) {
  381         throw new IllegalArgumentException();
  382       }
  383   
  384       int length1 = longFactorials.length;
  385       if (k < length1) {
  386         return longFactorials[k];
  387       }
  388   
  389       int length2 = doubleFactorials.length;
  390       if (k < length1 + length2) {
  391         return doubleFactorials[k - length1];
  392       } else {
  393         return Double.POSITIVE_INFINITY;
  394       }
  395     }
  396   
  397     /**
  398      * Returns the largest <code>long &lt;= value</code>. <dt>Examples: <code> 1.0 -> 1, 1.2 -> 1, 1.9 -> 1 <dt> 2.0 -> 2,
  399      * 2.2 -> 2, 2.9 -> 2 </code><dt> This method is safer than using (long) Math.floor(value), because of possible
  400      * rounding error.
  401      */
  402     public static long floor(double value) {
  403       return Math.round(Math.floor(value));
  404     }
  405   
  406     /** Returns <tt>log<sub>base</sub>value</tt>. */
  407     public static double log(double base, double value) {
  408       return Math.log(value) / Math.log(base);
  409     }
  410   
  411     /** Returns <tt>log<sub>10</sub>value</tt>. */
  412     public static double log10(double value) {
  413       // 1.0 / Math.log(10) == 0.43429448190325176
  414       return Math.log(value) * 0.43429448190325176;
  415     }
  416   
  417     /** Returns <tt>log<sub>2</sub>value</tt>. */
  418     public static double log2(double value) {
  419       // 1.0 / Math.log(2) == 1.4426950408889634
  420       return Math.log(value) * 1.4426950408889634;
  421     }
  422   
  423     /**
  424      * Returns <tt>log(k!)</tt>. Tries to avoid overflows. For <tt>k<30</tt> simply looks up a table in O(1). For
  425      * <tt>k>=30</tt> uses stirlings approximation.
  426      *
  427      * @param k must hold <tt>k &gt;= 0</tt>.
  428      */
  429     public static double logFactorial(int k) {
  430       if (k >= 30) {
  431   
  432         double r = 1.0 / (double) k;
  433         double rr = r * r;
  434         double C7 = -5.95238095238095238e-04;
  435         double C5 = 7.93650793650793651e-04;
  436         double C3 = -2.77777777777777778e-03;
  437         double C1 = 8.33333333333333333e-02;
  438         double C0 = 9.18938533204672742e-01;
  439         return (k + 0.5) * Math.log(k) - k + C0 + r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
  440       } else {
  441         return logFactorials[k];
  442       }
  443     }
  444   
  445     /**
  446      * Instantly returns the factorial <tt>k!</tt>.
  447      *
  448      * @param k must hold <tt>k &gt;= 0 && k &lt; 21</tt>.
  449      */
  450     public static long longFactorial(int k) throws IllegalArgumentException {
  451       if (k < 0) {
  452         throw new IllegalArgumentException("Negative k");
  453       }
  454   
  455       if (k < longFactorials.length) {
  456         return longFactorials[k];
  457       }
  458       throw new IllegalArgumentException("Overflow");
  459     }
  460   
  461     /**
  462      * Returns the StirlingCorrection. <p> Correction term of the Stirling approximation for <tt>log(k!)</tt> (series in
  463      * 1/k, or table values for small k) with int parameter k. <p> <tt> log k! = (k + 1/2)log(k + 1) - (k + 1) +
  464      * (1/2)log(2Pi) + stirlingCorrection(k + 1) <p> log k! = (k + 1/2)log(k)     -  k      + (1/2)log(2Pi) +
  465      * stirlingCorrection(k) </tt>
  466      */
  467     public static double stirlingCorrection(int k) {
  468   
  469       if (k > 30) {
  470         double r = 1.0 / (double) k;
  471         double rr = r * r;
  472         double C7 = -5.95238095238095238e-04;     //  -1/1680
  473         double C5 = 7.93650793650793651e-04;     //  +1/1260
  474         double C3 = -2.77777777777777778e-03;     //  -1/360
  475         double C1 = 8.33333333333333333e-02;     //  +1/12
  476         return r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
  477       } else {
  478         return stirlingCorrection[k];
  479       }
  480     }
  481   
  482   }

Save This Page
Home » mahout-collections-1.0-src » org.apache.mahout.collections » [javadoc | source]