Back to index

lightning-sunbird  0.9+nobinonly
nsIncompleteGamma.h
Go to the documentation of this file.
00001 /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
00002 /***** BEGIN LICENSE BLOCK *****
00003   Version: MPL 1.1/GPL 2.0/LGPL 2.1
00004 
00005   The contents of this file are subject to the Mozilla Public License Version 
00006   1.1 (the "License"); you may not use this file except in compliance with 
00007   the License. You may obtain a copy of the License at 
00008   http://www.mozilla.org/MPL/
00009 
00010   Software distributed under the License is distributed on an "AS IS" basis,
00011   WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00012   for the specific language governing rights and limitations under the
00013   License.
00014 
00015   The Original Code is bayesian spam filter
00016 
00017   The Initial Developer of the Original Code is
00018   Scott MacGregor.
00019   Portions created by the Initial Developer are Copyright (C) 2004
00020   the Initial Developer. All Rights Reserved.
00021 
00022   Contributor(s):
00023     Scott MacGregor <mscott@mozilla.org>
00024     tenthumbs@cybernex.net
00025 
00026   Alternatively, the contents of this file may be used under the terms of
00027   either of the GNU General Public License Version 2 or later (the "GPL"),
00028   or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029   in which case the provisions of the GPL or the LGPL are applicable instead
00030   of those above. If you wish to allow use of your version of this file only
00031   under the terms of either the GPL or the LGPL, and not to allow others to
00032   use your version of this file under the terms of the MPL, indicate your
00033   decision by deleting the provisions above and replace them with the notice
00034   and other provisions required by the GPL or the LGPL. If you do not delete
00035   the provisions above, a recipient may use your version of this file under
00036   the terms of any one of the MPL, the GPL or the LGPL.
00037 
00038 ***** END LICENSE BLOCK *****/
00039 
00040 #ifndef nsIncompleteGamma_h__
00041 #define nsIncompleteGamma_h__
00042 
00043 /* An implementation of the incomplete gamma functions for real
00044    arguments. P is defined as
00045 
00046                           x
00047                          /
00048                   1      [     a - 1  - t
00049     P(a, x) =  --------  I    t      e    dt
00050                Gamma(a)  ]
00051                          /
00052                           0
00053 
00054    and
00055 
00056                           infinity
00057                          /
00058                   1      [     a - 1  - t
00059     Q(a, x) =  --------  I    t      e    dt
00060                Gamma(a)  ]
00061                          /
00062                           x
00063 
00064    so that P(a,x) + Q(a,x) = 1.
00065 
00066    Both a series expansion and a continued fraction exist.  This
00067    implementation uses the more efficient method based on the arguments.
00068 
00069    Either case involves calculating a multiplicative term:
00070       e^(-x)*x^a/Gamma(a).
00071    Here we calculate the log of this term. Most math libraries have a
00072    "lgamma" function but it is not re-entrant. Some libraries have a
00073    "lgamma_r" which is re-entrant. Use it if possible. I have included a
00074    simple replacement but it is certainly not as accurate.
00075 
00076    Relative errors are almost always < 1e-10 and usually < 1e-14. Very
00077    small and very large arguments cause trouble.
00078 
00079    The region where a < 0.5 and x < 0.5 has poor error properties and is
00080    not too stable. Get a better routine if you need results in this
00081    region.
00082 
00083    The error argument will be set negative if there is a domain error or
00084    positive for an internal calculation error, currently lack of
00085    convergence. A value is always returned, though.
00086 
00087  */
00088 
00089 #include <math.h>
00090 #include <float.h>
00091 
00092 // the main routine
00093 static double nsIncompleteGammaP (double a, double x, int *error);
00094 
00095 // nsLnGamma(z): either a wrapper around lgamma_r or the internal function.
00096 // C_m = B[2*m]/(2*m*(2*m-1)) where B is a Bernoulli number
00097 static const double C_1 = 1.0 / 12.0;
00098 static const double C_2 = -1.0 / 360.0;
00099 static const double C_3 = 1.0 / 1260.0;
00100 static const double C_4 = -1.0 / 1680.0;
00101 static const double C_5 = 1.0 / 1188.0;
00102 static const double C_6 = -691.0 / 360360.0;
00103 static const double C_7 = 1.0 / 156.0;
00104 static const double C_8 = -3617.0 / 122400.0;
00105 static const double C_9 = 43867.0 / 244188.0;
00106 static const double C_10 = -174611.0 / 125400.0;
00107 static const double C_11 = 77683.0 / 5796.0;
00108 
00109 // truncated asymptotic series in 1/z
00110 static inline double lngamma_asymp (double z)
00111 {
00112   double w, w2, sum;
00113   w = 1.0 / z;
00114   w2 = w * w;
00115   sum = w * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2
00116         * (C_11 * w2 + C_10) + C_9) + C_8) + C_7) + C_6)
00117         + C_5) + C_4) + C_3) + C_2) + C_1);
00118 
00119   return sum;
00120 }
00121 
00122 struct fact_table_s
00123 {
00124   double fact;
00125   double lnfact;
00126 };
00127 
00128 // for speed and accuracy
00129 static const struct fact_table_s FactTable[] = {
00130     {1.000000000000000, 0.0000000000000000000000e+00},
00131     {1.000000000000000, 0.0000000000000000000000e+00},
00132     {2.000000000000000, 6.9314718055994530942869e-01},
00133     {6.000000000000000, 1.7917594692280550007892e+00},
00134     {24.00000000000000, 3.1780538303479456197550e+00},
00135     {120.0000000000000, 4.7874917427820459941458e+00},
00136     {720.0000000000000, 6.5792512120101009952602e+00},
00137     {5040.000000000000, 8.5251613610654142999881e+00},
00138     {40320.00000000000, 1.0604602902745250228925e+01},
00139     {362880.0000000000, 1.2801827480081469610995e+01},
00140     {3628800.000000000, 1.5104412573075515295248e+01},
00141     {39916800.00000000, 1.7502307845873885839769e+01},
00142     {479001600.0000000, 1.9987214495661886149228e+01},
00143     {6227020800.000000, 2.2552163853123422886104e+01},
00144     {87178291200.00000, 2.5191221182738681499610e+01},
00145     {1307674368000.000, 2.7899271383840891566988e+01},
00146     {20922789888000.00, 3.0671860106080672803835e+01},
00147     {355687428096000.0, 3.3505073450136888885825e+01},
00148     {6402373705728000., 3.6395445208033053576674e+01}
00149 };
00150 #define FactTableLength (int)(sizeof(FactTable)/sizeof(FactTable[0]))
00151 
00152 // for speed
00153 static const double ln_2pi_2 = 0.918938533204672741803; // log(2*PI)/2
00154 
00155 /* A simple lgamma function, not very robust.
00156 
00157    Valid for z_in > 0 ONLY.
00158 
00159    For z_in > 8 precision is quite good, relative errors < 1e-14 and
00160    usually better. For z_in < 8 relative errors increase but are usually
00161    < 1e-10. In two small regions, 1 +/- .001 and 2 +/- .001 errors
00162    increase quickly.
00163 */
00164 static double nsLnGamma (double z_in, int *gsign)
00165 {
00166   double scale, z, sum, result;
00167   *gsign = 1;
00168 
00169   int zi = (int) z_in;
00170   if (z_in == (double) zi)
00171   {
00172     if (0 < zi && zi <= FactTableLength)
00173            return FactTable[zi - 1].lnfact;      // gamma(z) = (z-1)!
00174   }
00175 
00176   for (scale = 1.0, z = z_in; z < 8.0; ++z)
00177     scale *= z;
00178   
00179   sum = lngamma_asymp (z);
00180   result = (z - 0.5) * log (z) - z + ln_2pi_2 - log (scale);
00181   result += sum;
00182   return result;
00183 }
00184 
00185 // log( e^(-x)*x^a/Gamma(a) )
00186 static inline double lnPQfactor (double a, double x)
00187 {
00188   int gsign;                       // ignored because a > 0
00189   return a * log (x) - x - nsLnGamma (a, &gsign);
00190 }
00191 
00192 static double Pseries (double a, double x, int *error)
00193 {
00194   double sum, term;
00195   const double eps = 2.0 * DBL_EPSILON;
00196   const int imax = 5000;
00197   int i;
00198 
00199   sum = term = 1.0 / a;
00200   for (i = 1; i < imax; ++i)
00201   {
00202          term *= x / (a + i);
00203          sum += term;
00204          if (fabs (term) < eps * fabs (sum))
00205       break;
00206   }
00207 
00208   if (i >= imax)
00209          *error = 1;
00210 
00211   return sum;
00212 }
00213 
00214 static double Qcontfrac (double a, double x, int *error)
00215 {
00216   double result, D, C, e, f, term;
00217   const double eps = 2.0 * DBL_EPSILON;
00218   const double small =
00219   DBL_EPSILON * DBL_EPSILON * DBL_EPSILON * DBL_EPSILON;
00220   const int imax = 5000;
00221   int i;
00222 
00223   // modified Lentz method
00224   f = x - a + 1.0;
00225   if (fabs (f) < small)
00226     f = small;
00227   C = f + 1.0 / small;
00228   D = 1.0 / f;
00229   result = D;
00230   for (i = 1; i < imax; ++i)
00231   {
00232     e = i * (a - i);
00233     f += 2.0;
00234     D = f + e * D;
00235     if (fabs (D) < small)
00236       D = small;
00237     D = 1.0 / D;
00238     C = f + e / C;
00239     if (fabs (C) < small)
00240       C = small;
00241     term = C * D;
00242     result *= term;
00243     if (fabs (term - 1.0) < eps)
00244       break;
00245   }
00246 
00247   if (i >= imax)
00248     *error = 1;
00249   return result;
00250 }
00251 
00252 static double nsIncompleteGammaP (double a, double x, int *error)
00253 {
00254   double result, dom, ldom;
00255   //  domain errors. the return values are meaningless but have
00256   //  to return something. 
00257   *error = -1;
00258   if (a <= 0.0)
00259          return 1.0;
00260   if (x < 0.0)
00261          return 0.0;
00262   *error = 0;
00263   if (x == 0.0)
00264     return 0.0;
00265 
00266   ldom = lnPQfactor (a, x);
00267   dom = exp (ldom);
00268   // might need to adjust the crossover point
00269   if (a <= 0.5)
00270   {
00271     if (x < a + 1.0)
00272       result = dom * Pseries (a, x, error);
00273     else
00274       result = 1.0 - dom * Qcontfrac (a, x, error);
00275   }
00276   else
00277   {
00278     if (x < a)
00279       result = dom * Pseries (a, x, error);
00280          else
00281            result = 1.0 - dom * Qcontfrac (a, x, error);
00282   }
00283 
00284   // not clear if this can ever happen
00285   if (result > 1.0)
00286     result = 1.0;
00287   if (result < 0.0)
00288     result = 0.0;
00289   return result;
00290 }
00291 
00292 #endif
00293