[Quantlib-dev] subtle bug in InverseCumulativeNormal

[Kakhkhor A. <kab...@gm...>]

Bug Report:

Where:            In "normaldistribution.cpp".
What:              Undefined behavior or floating point exception.
Severity:          Subtle.
Burden:           Easy to fix. See below.


Description.
If the argument x is close to 0.0 or 1.0, the implementation sets it
exactly to 0.0 or 1.0.
This leads to an attempt to evaluate std::log(0.0) or std::log(1.0).
I propose that we cut it at x=1E-12 and x = 1 - 1E-12, which covers
all draws between
plus/minus 7 standard deviations. Any argument value beyond that range
should be considered
as either erroneous or astronomically improbable.

Regards,
Kakhkhor Abdijalilov.

=============================================================
// current implementation
Real InverseCumulativeNormal::operator()(Real x) const {
    if (x < 0.0 || x > 1.0) {
        // try to recover if due to numerical error
        if (close_enough(x, 1.0)) {
            x = 1.0;
        } else if (std::fabs(x) < QL_EPSILON) {
            x = 0.0;
        } else {
            QL_FAIL("InverseCumulativeNormal(" << x
                    << ") undefined: must be 0 < x < 1");
        }
    }
.........................................
}

// new implementation
Real InverseCumulativeNormal::operator()(Real x) const {
    if (x < 1e-12 || x > (1.0 - 1e-12)) {
        // try to recover if due to numerical error
        if (close_enough(x, 1.0)) {
            x = 1.0 - 1e-12;
        } else if (std::fabs(x) < QL_EPSILON) {
            x = 1e-12;
        } else {
            QL_FAIL("InverseCumulativeNormal(" << x
                    << ") undefined: must be 0 < x < 1");
        }
    }
.........................................
}