Cause of a stack overflow in this method (floating-point math)

Question

I occasionally 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 might be getting it ? Any steps I can take to rectify the error ?

Answer 1

Your problem is when x is Not a Number. NAN >= 0.0 is false, -NAN >= 0.0 is also false.

You could check against NAN specially, as others have 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 smarter assumptions about its behaviour. Note that I've marked the "internal" function static (which will limit its scope to the current compilation unit). You could also use unnamed namespaces. (edit: actually Timothy Shields has a simpler way of removing the recursion, which keeps everything in one function)