/[gsl]/gsl/cdf/gammainv.c
ViewVC logotype

Diff of /gsl/cdf/gammainv.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.1 by bjg, Sun Jun 29 13:33:36 2003 UTC revision 1.1.2.2 by bjg, Sun Jun 29 17:31:38 2003 UTC
# Line 23  Line 23 
23  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
24  #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
25    
26    #include <stdio.h>
27    
28  double  double
29  gsl_cdf_gamma_Pinv (double p, double a, double b)  gsl_cdf_gamma_Pinv (double p, double a, double b)
30  {  {
31    double x;    double x, xc;
32    
33    if (p == 1.0)    if (p == 1.0)
34      {      {
# Line 37  gsl_cdf_gamma_Pinv (double p, double a, Line 39  gsl_cdf_gamma_Pinv (double p, double a,
39        return 0.0;        return 0.0;
40      }      }
41    
   if (p < 0.05)  
     {  
       double x0 = exp((lgamma(b) + log(p))/a);  
       x = x0;  
     }  
   else if (p > 0.95)  
     {  
       double x0 = exp(-log(1-p)-lgamma(b));  
       x = x0;  
     }  
   else  
     {  
       double xg = gsl_cdf_ugaussian_Pinv (p);  
       double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b;  
       x = x0;  
     }  
   
   {  
     double dp = p - gsl_cdf_gamma_P (x, 1.0, b);  
     double u = log(x);  
   
     while (fabs (dp) > 1e-10 * p)  
       {  
         double dp1, step = dp / exp((b-1.0)*u);  
   
         do  
           {  
             u = u + step;  
             x = exp(u);  
             dp1 = p - gsl_cdf_gamma_P (x, 1.0, b);  
             step *= 2.0;  
           }  
         while((dp > 0 && dp1 > 0) || (dp < 0 && dp1 < 0));  
   
         dp = dp1;  
       }  
   }  
42    
43    return a * x;    return a * x;
44  }  }

Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.1.2.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26