The reason for this method (floating point math)

Question

Sometimes I get a stackoverflow exception in this method.

double norm_cdf(const double x) { double k = 1.0/(1.0 + 0.2316419*x); double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k)))); if (x >= 0.0) { return (1.0 - (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum); } else { return 1.0 - norm_cdf(-x); } } 

Any suggestions on why I can get it? Any steps I can take to fix the error?

Answer 1

This method is probably just straw that breaks the camel back. This feature will only ring more than once, so this is not a problem. (Edit: Or this is indicated by other NAN users, which will result in infinite recursion.)

You can easily make the function non-recursive anyway, which could be simpler.

 double norm_cdf(double x) { bool negative = x < 0; x = abs(x); double k = 1.0/(1.0 + 0.2316419*x); double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k)))); double result = (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum; if (!negative) result = 1.0 - result; return result; } 

Answer 2

Your problem is that x not a number. NAN >= 0.0 is false, -NAN >= 0.0 also false.

You could check against NAN specifically, as others suggested, but I would suggest simplifying things:

 static double norm_cdf_positive(const double x) { double k = 1.0/(1.0 + 0.2316419*x); double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k)))); return (1.0 - (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum); } double norm_cdf(const double x) { if (x >= 0.0) { return norm_cdf_positive(x); } else { return 1.0 - norm_cdf_positive(-x); } } 

This has the advantage that compilers can make more reasonable assumptions about their behavior. Note that I marked the “internal” static function (which will limit the scope to the current compilation unit). You can also use unnamed namespaces. (edit: Timothy Shields actually has an easier way to remove recursion, which saves everything in one function)

Answer 3

When x NaN , the recursion will never end. Add a check: either std::isnan in C ++ 11, or lazy x != x , or a document-and-blame-user. If you decide to process NaN , distributing it may be a smart choice.