LCOV - code coverage report
Current view: directory - src/vxl/core/vnl - vnl_gamma.cxx (source / functions) Found Hit Coverage
Test: vnl.info Lines: 69 45 65.2 %
Date: 2011-09-02 Functions: 9 6 66.7 %

       1                 : // This is core/vnl/vnl_gamma.cxx
       2                 : #include "vnl_gamma.h"
       3                 : //:
       4                 : // \file
       5                 : // \brief Complete and incomplete gamma function approximations
       6                 : // \author Tim Cootes
       7                 : 
       8                 : #include <vcl_iostream.h>
       9                 : #include <vcl_cassert.h>
      10                 : 
      11                 : #if defined(__INTEL_COMPILER)
      12                 : # pragma warning (disable:279) /* controlling expression is constant */
      13                 : #endif
      14                 : 
      15                 : //: Approximate gamma function
      16                 : //  Uses 6 parameter Lanczos approximation as described by Viktor Toth
      17                 : //  (http://www.rskey.org/gamma.htm)
      18                 : //  Accurate to about 3e-11.
      19              75 : double vnl_log_gamma(double x)
      20                 : {
      21              75 :   double zp = 2.50662827563479526904;
      22              75 :   zp += 225.525584619175212544/x;
      23              75 :   zp -= 268.295973841304927459/(x+1.0);
      24              75 :   zp += 80.9030806934622512966/(x+2.0);
      25              75 :   zp -= 5.00757863970517583837/(x+3.0);
      26              75 :   zp += 0.0114684895434781459556/(x+4.0);
      27                 : 
      28              75 :   double x1 = x+4.65;
      29                 : 
      30              75 :   return vcl_log(zp)+(x-0.5)*vcl_log(x1)-x1;
      31                 : }
      32                 : 
      33                 : const int MAX_ITS = 100;
      34                 : const double MaxRelError = 3.0e-7;
      35                 : const double vnl_very_small = 1.0e-30;
      36                 : 
      37                 : //: Use series expansion of incomplete gamma function
      38              47 : static double vnl_gamma_series(double a, double x)
      39                 : {
      40              47 :   if (x>0)
      41                 :   {
      42              43 :     double a_i=a;
      43              43 :     double term_i=1.0/a;
      44              43 :     double sum = term_i;
      45             154 :     for (int i=1;i<=MAX_ITS;++i)
      46                 :     {
      47             154 :       a_i+=1;
      48             154 :       term_i *= x/a_i;
      49             154 :       sum += term_i;
      50             154 :       if (vcl_fabs(term_i) < vcl_fabs(sum)*MaxRelError)
      51              43 :         return sum*vcl_exp(-x+a*vcl_log(x)-vnl_log_gamma(a));
      52                 :     }
      53                 :     vcl_cerr<<"vnl_gamma_series : Failed to converge in "<<MAX_ITS<<" steps\n"
      54               0 :             <<"a = "<<a<<"   x= "<< x <<"\nReturning best guess.\n";
      55               0 :     return sum*vcl_exp(-x+a*vcl_log(x)-vnl_log_gamma(a));
      56                 :   }
      57               4 :   else if (x < 0.0)
      58               0 :     assert(!"vnl_gamma_series - x less than 0");
      59                 : 
      60               4 :   return 0.0;
      61                 : }
      62                 : 
      63                 : //: Incomplete gamma using continued fraction representation
      64                 : // Use Lentz's algorithm
      65                 : // Continued fraction with terms a_i/b_i
      66                 : // a_i = i*(a-i), b_i = (x+a-2i)
      67              32 : static double vnl_gamma_cont_frac(double a, double x)
      68                 : {
      69              32 :   double b_i=x+1.0-a;
      70              32 :   double c=1.0/vnl_very_small;
      71              32 :   double d=1.0/b_i;
      72              32 :   double cf=d;
      73             167 :   for (int i=1;i<=MAX_ITS;i++)
      74                 :   {
      75             167 :     double a_i = i*(a-i);
      76             167 :     b_i += 2.0;
      77             167 :     d=a_i*d+b_i;
      78             167 :     if (vcl_fabs(d) < vnl_very_small) d=vnl_very_small;
      79             167 :     c=b_i+a_i/c;
      80             167 :     if (vcl_fabs(c) < vnl_very_small) c=vnl_very_small;
      81             167 :     d=1.0/d;
      82             167 :     double delta=d*c;
      83             167 :     cf *= delta;
      84             167 :     if (vcl_fabs(delta-1.0) < MaxRelError)
      85              32 :       return vcl_exp(-x+a*vcl_log(x)-vnl_log_gamma(a))*cf;
      86                 :   }
      87                 : 
      88                 :   vcl_cerr<<"vnl_gamma_cont_frac : Failed to converge in "<<MAX_ITS<<" steps\n"
      89               0 :           <<"a = "<<a<<"   x= "<<x<<vcl_endl;
      90               0 :   return vcl_exp(-x+a*vcl_log(x)-vnl_log_gamma(a))*cf;
      91                 : }
      92                 : 
      93              79 : double vnl_gamma_p(double a, double x)
      94                 : {
      95              79 :   if (x < 0.0 || a <= 0.0)
      96               0 :     assert(!"vnl_gamma_p - Invalid arguments.");
      97                 : 
      98              79 :   if (x < a+1.0)
      99              47 :     return vnl_gamma_series(a,x); // Use series representation
     100                 :   else
     101              32 :     return 1.0 - vnl_gamma_cont_frac(a,x); // Use continued fraction representation
     102                 : }
     103                 : 
     104               0 : double vnl_gamma_q(double a, double x)
     105                 : {
     106               0 :   if (x < 0.0 || a <= 0.0)
     107               0 :     assert(!"vnl_gamma_q - Invalid arguments.");
     108                 : 
     109               0 :   if (x < a+1.0)
     110               0 :     return 1.0-vnl_gamma_series(a,x); // Use series representation
     111                 :   else
     112               0 :     return vnl_gamma_cont_frac(a,x); // Use continued fraction representation
     113                 : }
     114                 : 
     115               0 : double vnl_digamma(double z)
     116                 : {
     117               0 :   double t0 = (z-0.5)/(z+4.65)-1.0;
     118               0 :   double tlg = vcl_log(4.65+z);
     119               0 :   double tc = 2.50662827563479526904;
     120               0 :   double t1 = 225.525584619175212544/z;
     121               0 :   double t2 = -268.295973841304927459/(1+z);
     122               0 :   double t3 = +80.9030806934622512966/(2+z);
     123               0 :   double t4 = -5.00757863970517583837/(3+z);
     124               0 :   double t5 = 0.0114684895434781459556/(4+z);
     125               0 :   double neu = t1/z + t2/(1+z) + t3/(2+z) + t4/(3+z) + t5/(4+z);
     126               0 :   double den = tc + t1 + t2 + t3 + t4 + t5;
     127               0 :   return (t0 -(neu/den) + tlg);
     128            6315 : }

Generated by: LCOV version 1.7