/[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.11 by bjg, Sat Mar 1 20:36:13 2003 UTC revision 1.1.2.12 by bjg, Sat Mar 1 20:58:33 2003 UTC
# Line 46  Line 46 
46  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
47  #include <gsl/gsl_cdf.h>  #include <gsl/gsl_cdf.h>
48    
49  #ifndef M_1_SQRT2PI  #ifndef M_1_SQRT2PI
50  #define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)  #define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)
51  #endif  #endif
52    
53    #define SQRT32 (4.0 * M_SQRT2)
54    
55  /*  /*
56   * IEEE double precision dependent constants.   * IEEE double precision dependent constants.
57   *   *
# Line 60  Line 62 
62   */   */
63    
64  #define GAUSS_EPSILON  GSL_DBL_EPSILON / 2  #define GAUSS_EPSILON  GSL_DBL_EPSILON / 2
65  #define GAUSS_XUPPER 8.572  #define GAUSS_XUPPER 8.572
66  #define GAUSS_XLOWER -37.519  #define GAUSS_XLOWER -37.519
67    
68  #define GAUSS_SCALE 16.0  #define GAUSS_SCALE 16.0
69    
70  static double get_del ( double x, double rational )  static double
71    get_del (double x, double rational)
72  {  {
73    double xsq = 0.0;    double xsq = 0.0;
74    double del = 0.0;    double del = 0.0;
75    double result = 0.0;    double result = 0.0;
76    
77    xsq = floor( x * GAUSS_SCALE ) / GAUSS_SCALE;    xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE;
78    del = (x - xsq) * (x + xsq);    del = (x - xsq) * (x + xsq);
79    del *= 0.5;    del *= 0.5;
80    
81    result = exp(-0.5 * xsq * xsq ) * exp(-1.0 * del) * rational;    result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational;
82    
83    return result;    return result;
84  }  }
85    
86  /*  /*
87   * Normal cdf for fabs(x) < 0.66291   * Normal cdf for fabs(x) < 0.66291
88   */   */
89  static double gauss_smallx(const double x)  static double
90    gauss_smallx (const double x)
91  {  {
92      unsigned int i;
93    double result = 0.0;    double result = 0.0;
   double temp = 0.0;  
94    double xsq;    double xsq;
95    double xnum;    double xnum;
96    double xden;    double xden;
# Line 103  static double gauss_smallx(const double Line 108  static double gauss_smallx(const double
108      10260.932208618978205,      10260.932208618978205,
109      45507.789335026729956      45507.789335026729956
110    };    };
   int n;  
111    
112    xsq = x * x;    xsq = x * x;
113    xnum = a[4] * xsq;    xnum = a[4] * xsq;
114    xden = xsq;    xden = xsq;
115    for (n = 0; n < 3; ++n) {  
116      xnum = (xnum + a[n]) * xsq;    for (i = 0; i < 3; i++)
117      xden = (xden + b[n]) * xsq;      {
118    }        xnum = (xnum + a[i]) * xsq;
119          xden = (xden + b[i]) * xsq;
120        }
121    
122    result = x * (xnum + a[3]) / (xden + b[3]);    result = x * (xnum + a[3]) / (xden + b[3]);
     
   return result;  
123    
124      return result;
125  }  }
126    
127  /*  /*
128   * Normal cdf for 0.66291 < fabs(x) < sqrt(32).   * Normal cdf for 0.66291 < fabs(x) < sqrt(32).
129   */   */
130  static double gauss_mediumx ( const double x)  static double
131    gauss_mediumx (const double x)
132  {  {
133    int n;    unsigned int i;
134    double temp = 0.0;    double temp = 0.0;
135    double result = 0.0;    double result = 0.0;
136    double xnum;    double xnum;
# Line 151  static double gauss_mediumx ( const doub Line 158  static double gauss_mediumx ( const doub
158      38912.003286093271411,      38912.003286093271411,
159      19685.429676859990727      19685.429676859990727
160    };    };
161    absx = fabs(x);  
162      absx = fabs (x);
163    
164    xnum = c[8] * absx;    xnum = c[8] * absx;
165    xden = absx;    xden = absx;
166    for (n = 0; n < 7; ++n) {  
167      xnum = (xnum + c[n]) * absx;    for (i = 0; i < 7; i++)
168      xden = (xden + d[n]) * absx;      {
169    }        xnum = (xnum + c[i]) * absx;
170          xden = (xden + d[i]) * absx;
171        }
172    
173    temp = (xnum + c[7]) / (xden + d[7]);    temp = (xnum + c[7]) / (xden + d[7]);
174    
175    result = get_del(x, temp);    result = get_del (x, temp);
176    
177    return result;    return result;
178  }  }
# Line 170  static double gauss_mediumx ( const doub Line 181  static double gauss_mediumx ( const doub
181   * Normal cdf for   * Normal cdf for
182   * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.   * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
183   */   */
184  static double gauss_grandex ( const double x )  static double
185    gauss_grandex (const double x)
186  {  {
187    int n;    int i;
188    double result;    double result;
189    double xsq;    double xsq;
190    double temp;    double temp;
191    double xnum;    double xnum;
192    double xden;    double xden;
193    double absx;    double absx;
194      
195    const double p[6] = {    const double p[6] = {
196      0.21589853405795699,      0.21589853405795699,
197      0.1274011611602473639,      0.1274011611602473639,
# Line 196  static double gauss_grandex ( const doub Line 208  static double gauss_grandex ( const doub
208      7.29751555083966205e-5      7.29751555083966205e-5
209    };    };
210    
211    absx = fabs(x);    absx = fabs (x);
212    xsq = 1.0 / (x * x);    xsq = 1.0 / (x * x);
213    xnum = p[5] * xsq;    xnum = p[5] * xsq;
214    xden = xsq;    xden = xsq;
215    for (n = 0; n < 4; ++n) {  
216      xnum = (xnum + p[n]) * xsq;    for (i = 0; i < 4; i++)
217      xden = (xden + q[n]) * xsq;      {
218    }        xnum = (xnum + p[i]) * xsq;
219          xden = (xden + q[i]) * xsq;
220        }
221    
222    temp = xsq * (xnum + p[4]) / (xden + q[4]);    temp = xsq * (xnum + p[4]) / (xden + q[4]);
223    temp = (M_1_SQRT2PI - temp) / absx;    temp = (M_1_SQRT2PI - temp) / absx;
224    
225    result = get_del(x, temp);    result = get_del (x, temp);
226    
227    return result;    return result;
228  }  }
229    
230  double gsl_cdf_gaussian_Q (const double x )  double
231    gsl_cdf_gaussian_P (const double x)
232  {  {
233    double result;    double result;
234    double absx;    double absx = fabs (x);
   double sqrt32 = sqrt(32);  
   
   absx = fabs(x);  
235    
236    if ( absx < 0.66291 )    if (absx < GAUSS_EPSILON)
237      {      {
238        result = gauss_smallx( x);        result = 0.5;
239        if ( x < 0.0 )        return result;
         {  
           result = fabs(result) + 0.5;  
         }  
       else  
         result = 0.5 - result;  
         
240      }      }
241    else if ( absx < sqrt(32))    else if (absx < 0.66291)
242      {      {
243        result = gauss_mediumx(x);        result = 0.5 + gauss_smallx (x);
244        if( x < 0.0 )        return result;
         {  
           result = 1.0 - result;  
         }  
245      }      }
246    else if ( x >= -GAUSS_XUPPER || x <= -GAUSS_XLOWER )    else if (absx < SQRT32)
247      {      {
248        result = gauss_grandex (x);        result = gauss_mediumx (x);
249        if ( x < 0.0 )  
250          if (x > 0.0)
251          {          {
252            result = 1.0 - result;            result = 1.0 - result;
253          }          }
254    
255          return result;
256      }      }
257    else if ( x > -GAUSS_XLOWER )    else if (x > GAUSS_XUPPER)
258      {      {
259        result = 0.0;        result = 1.0;
260          return result;
261      }      }
262    else if ( x < -GAUSS_XUPPER )    else if (x < GAUSS_XLOWER)
263      {      {
264        result = 1.0;        result = 0.0;
265          return result;
266      }      }
267    else if ( x > 0.0 && x < GAUSS_EPSILON )    else
268      {      {
269        result = 0.5;        result = gauss_grandex (x);
270    
271          if (x > 0.0)
272            {
273              result = 1.0 - result;
274            }
275      }      }
276    
277    return result;    return result;
278  }  }
279  double gsl_cdf_gaussian_P(const double x)  
280    double
281    gsl_cdf_gaussian_Q (const double x)
282  {  {
283    double result;    double result;
284    double absx;    double absx = fabs (x);
   double sqrt32 = sqrt(32);  
285    
286    absx = fabs(x);    if (absx < GAUSS_EPSILON)
   
   if ( absx < 0.66291 )  
287      {      {
288        result = 0.5 + gauss_smallx( x);        result = 0.5;
289          return result;
290      }      }
291    else if ( absx < sqrt(32))    else if (absx < 0.66291)
292      {      {
293        result = gauss_mediumx(x);        result = gauss_smallx (x);
294        if ( x > 0.0 )  
295          if (x < 0.0)
296          {          {
297            result = 1.0 - result;            result = fabs (result) + 0.5;
298            }
299          else
300            {
301              result = 0.5 - result;
302          }          }
303    
304          return result;
305      }      }
306    else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER )    else if (absx < SQRT32)
307      {      {
308        result = gauss_grandex (x);        result = gauss_mediumx (x);
309        if( x > 0.0 )  
310          if (x < 0.0)
311          {          {
312            result = 1.0 - result;            result = 1.0 - result;
313          }          }
314    
315          return result;
316      }      }
317    else if ( x > GAUSS_XUPPER )    else if (x > (-GAUSS_XLOWER))
318      {      {
319        result = 1.0;        result = 0.0;
320          return result;
321      }      }
322    else if ( x < GAUSS_XLOWER )    else if (x < (-GAUSS_XUPPER))
323      {      {
324        result = 0.0;        result = 1.0;
325          return result;
326      }      }
327    else if ( x > 0.0 && x < GAUSS_EPSILON )    else
328      {      {
329        result = 0.5;        result = gauss_grandex (x);
330    
331          if (x < 0.0)
332            {
333              result = 1.0 - result;
334            }
335    
336      }      }
337    
338    return result;    return result;

Legend:
Removed from v.1.1.2.11  
changed lines
  Added in v.1.1.2.12

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