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

Diff of /gsl/cdf/gauss.c

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

revision 1.1.2.6 by jstover, Mon Feb 17 23:08:54 2003 UTC revision 1.1.2.7 by jstover, Sun Feb 23 02:11:42 2003 UTC
# Line 130  static double gauss_smallx(const double Line 130  static double gauss_smallx(const double
130  /*  /*
131   * Normal cdf for 0.66291 < fabs(x) < sqrt(32).   * Normal cdf for 0.66291 < fabs(x) < sqrt(32).
132   */   */
133  static double gauss_mediumx_Q ( const double x)  static double gauss_mediumx ( const double x)
134  {  {
135    int n;    int n;
136    double temp = 0.0;    double temp = 0.0;
# Line 172  static double gauss_mediumx_Q ( const do Line 172  static double gauss_mediumx_Q ( const do
172    
173    result = get_del(x, temp);    result = get_del(x, temp);
174    
   if ( x < 0.0 )  
     {  
       result = 1.0 - result;  
     }  
   return result;  
 }  
 static double gauss_mediumx_P ( const double x)  
 {  
   int n;  
   double temp = 0.0;  
   double result = 0.0;  
   double xnum;  
   double xden;  
   double absx;  
   
   const double c[9] = {  
     0.39894151208813466764,  
     8.8831497943883759412,  
     93.506656132177855979,  
     597.27027639480026226,  
     2494.5375852903726711,  
     6848.1904505362823326,  
     11602.651437647350124,  
     9842.7148383839780218,  
     1.0765576773720192317e-8  
   };  
   const double d[8] = {  
     22.266688044328115691,  
     235.38790178262499861,  
     1519.377599407554805,  
     6485.558298266760755,  
     18615.571640885098091,  
     34900.952721145977266,  
     38912.003286093271411,  
     19685.429676859990727  
   };  
   absx = fabs(x);  
   
   xnum = c[8] * absx;  
   xden = absx;  
   for (n = 0; n < 7; ++n) {  
     xnum = (xnum + c[n]) * absx;  
     xden = (xden + d[n]) * absx;  
   }  
   temp = (xnum + c[7]) / (xden + d[7]);  
   
   result = get_del(x, temp);  
   
   if( x > 0.0 )  
     {  
       result = 1.0 - result;  
     }  
   
175    return result;    return result;
176  }  }
177    
# Line 232  static double gauss_mediumx_P ( const do Line 179  static double gauss_mediumx_P ( const do
179   * Normal cdf for   * Normal cdf for
180   * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.   * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
181   */   */
182  static double gauss_grandex_Q ( const double x )  static double gauss_grandex ( const double x )
 {  
   int n;  
   double result;  
   double xsq;  
   double temp;  
   double xnum;  
   double xden;  
   double absx;  
     
   const double p[6] = {  
     0.21589853405795699,  
     0.1274011611602473639,  
     0.022235277870649807,  
     0.001421619193227893466,  
     2.9112874951168792e-5,  
     0.02307344176494017303  
   };  
   const double q[5] = {  
     1.28426009614491121,  
     0.468238212480865118,  
     0.0659881378689285515,  
     0.00378239633202758244,  
     7.29751555083966205e-5  
   };  
   
   absx = fabs(x);  
   xsq = 1.0 / (x * x);  
   xnum = p[5] * xsq;  
   xden = xsq;  
   for (n = 0; n < 4; ++n) {  
     xnum = (xnum + p[n]) * xsq;  
     xden = (xden + q[n]) * xsq;  
   }  
   temp = xsq * (xnum + p[4]) / (xden + q[4]);  
   temp = (M_1_SQRT2PI - temp) / absx;  
   
   result = get_del(x, temp);  
   
   if( x < 0.0 )  
     {  
       result = 1.0 - result;  
     }  
   return result;  
 }  
 static double gauss_grandex_P ( const double x )  
183  {  {
184    int n;    int n;
185    double result;    double result;
# Line 316  static double gauss_grandex_P ( const do Line 218  static double gauss_grandex_P ( const do
218    
219    result = get_del(x, temp);    result = get_del(x, temp);
220    
   /* If result is the upper tail area, subtract it from 1.0. */  
   if( x > 0.0 )  
     {  
       result = 1.0 - result;  
     }  
   
221    return result;    return result;
222  }  }
223    
# Line 346  double gsl_cdf_gauss_Q (const double x ) Line 242  double gsl_cdf_gauss_Q (const double x )
242      }      }
243    else if ( absx < sqrt(32))    else if ( absx < sqrt(32))
244      {      {
245        result = gauss_mediumx_Q(x);        result = gauss_mediumx(x);
246          if( x < 0.0 )
247            {
248              result = 1.0 - result;
249            }
250      }      }
251    else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER )    else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER )
252      {      {
253        result = gauss_grandex_Q (x);        result = gauss_grandex (x);
254          if ( x < 0.0 )
255            {
256              result = 1.0 - result;
257            }
258      }      }
259    else if ( x > GAUSS_XUPPER )    else if ( x > GAUSS_XUPPER )
260      {      {
# Line 381  double gsl_cdf_gauss_P(const double x) Line 285  double gsl_cdf_gauss_P(const double x)
285      }      }
286    else if ( absx < sqrt(32))    else if ( absx < sqrt(32))
287      {      {
288        result = gauss_mediumx_P(x);        result = gauss_mediumx(x);
289          if ( x > 0.0 )
290            {
291              result = 1.0 - result;
292            }
293      }      }
294    else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER )    else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER )
295      {      {
296        result = gauss_grandex_P (x);        result = gauss_grandex (x);
297          if( x > 0.0 )
298            {
299              result = 1.0 - result;
300            }
301      }      }
302    else if ( x > GAUSS_XUPPER )    else if ( x > GAUSS_XUPPER )
303      {      {

Legend:
Removed from v.1.1.2.6  
changed lines
  Added in v.1.1.2.7

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