Mercurial > hg > octave-lyh

comparison liboctave/lo-specfun.cc @ 14781:e190f6da40f6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Correct comments and use Octave spacing conventions for erfinv. lo-specfun.cc (do_erfinv): Correct comment about Halley's 3rd order method. Use space between function call and open parenthesis.
author Rik <octave@nomad.inbox5.com>
date Tue, 19 Jun 2012 10:05:01 -0700
parents 10ed11922f19
children e70a0c9cada6
comparison
equal deleted inserted replaced
14780:1d83d1539b2b 14781:e190f6da40f6
3150 } 3150 }
3151 3151
3152 // This algorithm is due to P. J. Acklam. 3152 // This algorithm is due to P. J. Acklam.
3153 // See http://home.online.no/~pjacklam/notes/invnorm/ 3153 // See http://home.online.no/~pjacklam/notes/invnorm/
3154 // The rational approximation has relative accuracy 1.15e-9 in the whole region. 3154 // The rational approximation has relative accuracy 1.15e-9 in the whole region.
3155 // For doubles, it is refined by a single step of Higham's 3rd order method. 3155 // For doubles, it is refined by a single step of Halley's 3rd order method.
3156 // For single precision, the accuracy is already OK, so we skip it to get 3156 // For single precision, the accuracy is already OK, so we skip it to get
3157 // faster evaluation. 3157 // faster evaluation.
3158 3158
3159 static double do_erfinv (double x, bool refine) 3159 static double do_erfinv (double x, bool refine)
3160 { 3160 {
3173 3.093354679843505e+00, 2.077595676404383e+00 }; 3173 3.093354679843505e+00, 2.077595676404383e+00 };
3174 static const double d[] = 3174 static const double d[] =
3175 { 7.784695709041462e-03, 3.224671290700398e-01, 3175 { 7.784695709041462e-03, 3.224671290700398e-01,
3176 2.445134137142996e+00, 3.754408661907416e+00 }; 3176 2.445134137142996e+00, 3.754408661907416e+00 };
3177 3177
3178 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. 3178 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
3179 static const double pbreak = 0.95150; 3179 static const double pbreak = 0.95150;
3180 double ax = fabs (x), y; 3180 double ax = fabs (x), y;
3181 3181
3182 // Select case. 3182 // Select case.
3183 if (ax <= pbreak) 3183 if (ax <= pbreak)
3202 return octave_NaN; 3202 return octave_NaN;
3203 3203
3204 if (refine) 3204 if (refine)
3205 { 3205 {
3206 // One iteration of Halley's method gives full precision. 3206 // One iteration of Halley's method gives full precision.
3207 double u = (erf(y) - x) * spi2 * exp (y*y); 3207 double u = (erf (y) - x) * spi2 * exp (y*y);
3208 y -= u / (1 + y*u); 3208 y -= u / (1 + y*u);
3209 } 3209 }
3210 3210
3211 return y; 3211 return y;
3212 } 3212 }