[anveshachuri]

fixes #2969.

C(), D(), E() and F() each form products like pow(H/S, 2 mu) * N(y), and returned 0 whenever N(y) underflowed. at low volatility and large carry pow overflows at the same time, and the product is finite and not small (0.076 for D in the issue's example), so the engine was off by that amount.

this adds a helper that computes power * N(y) as before when both factors are in range, so results don't change anywhere else (bit for bit), and otherwise forms it as exp(log(power) + log N(y)). for N(y) below the smallest normal double, log N(y) uses the log of the same asymptotic expansion (A&S 26.2.12) that CumulativeNormalDistribution already uses for very negative y, with the log of the gaussian density in closed form.

the issue only checked the down-and-out put without rebate, but the down-and-in put is off by the same amount (e.g. 65.5304 vs 65.6064), and with a rebate E() and F() hit the same problem (down-and-out put with rebate 2.0 was 3.6% too high at 1.6% vol). the new test testLowVolatilityLargeCarry covers down-and-out and down-and-in puts, with and without rebate, at 0.5% to 2% volatility (20 cases), against the same closed-form formulas evaluated with 60-digit arithmetic in mpmath (script below). before the change 16 of the 20 cases fail (every vol below 2%), after it all pass.

the tolerance is 1e-10 relative to max(price, 1): for the 0.5% down-and-out the price is 7.6e-6 but it's a difference of terms of the order of the spot, so 1e-14 absolute is the best a double can do there.

BarrierOptionTests (including the existing testLowVolatility at 1e-7 vol), SoftBarrierOptionTests, QuantoOptionTests, PartialTimeBarrierOptionTests, DoubleBarrierOptionTests, BinaryOptionTests and FdHestonTests/testFdmHestonBarrierVsBlackScholes pass (gcc 14.3, cmake + ninja, release).

ai use: I used Claude (Anthropic) to help write the fix, the test and the reference script; I checked the math and ran the tests before and after the change myself.

reference values (mpmath, 60 digits)

import mpmath as mp
import QuantLib as ql
S, K, H = 40.809474692328266, 100.0, 26.788905850262115
r, q = 0.04943364523655579, 0.27314598534530515
days = 747

mp.mp.dps = 60


def ql_npv(sigma, btype, rebate):
    today = ql.Date(15, 10, 2026)
    ql.Settings.instance().evaluationDate = today
    dc = ql.Actual365Fixed()
    p = ql.BlackScholesMertonProcess(
        ql.QuoteHandle(ql.SimpleQuote(S)),
        ql.YieldTermStructureHandle(ql.FlatForward(today, q, dc)),
        ql.YieldTermStructureHandle(ql.FlatForward(today, r, dc)),
        ql.BlackVolTermStructureHandle(ql.BlackConstantVol(today, ql.NullCalendar(), sigma, dc)))
    o = ql.BarrierOption(btype, H, rebate, ql.PlainVanillaPayoff(ql.Option.Put, K),
                         ql.EuropeanExercise(today + days))
    o.setPricingEngine(ql.AnalyticBarrierEngine(p))
    return o.NPV()


def ref(sigma, rebate):
    s, k, h, rb = mp.mpf(S), mp.mpf(K), mp.mpf(H), mp.mpf(rebate)
    rr, qq, vol = mp.mpf(r), mp.mpf(q), mp.mpf(sigma)
    t = mp.mpf(days) / 365
    sd = vol * mp.sqrt(t)
    mu = (rr - qq) / vol**2 - mp.mpf(1) / 2
    lam = mp.sqrt(mu**2 + 2 * rr / vol**2)
    musd = (1 + mu) * sd
    dr, dq = mp.exp(-rr * t), mp.exp(-qq * t)
    N = mp.ncdf
    phi, eta = -1, 1
    hs = h / s
    def AB(x): return phi * (s * dq * N(phi * x) - k * dr * N(phi * (x - sd)))
    def CD(y): return phi * (s * dq * hs ** (2 * (mu + 1)) * N(eta * y) - k * dr * hs ** (2 * mu) * N(eta * (y - sd)))
    x1 = mp.log(s / k) / sd + musd; x2 = mp.log(s / h) / sd + musd
    y1 = mp.log(h * hs / k) / sd + musd; y2 = mp.log(hs) / sd + musd
    A, B, C, D = AB(x1), AB(x2), CD(y1), CD(y2)
    E = rb * dr * (N(eta * (x2 - sd)) - hs ** (2 * mu) * N(eta * (y2 - sd)))
    z = mp.log(hs) / sd + lam * sd
    F = rb * (hs ** (mu + lam) * N(eta * z) + hs ** (mu - lam) * N(eta * (z - 2 * lam * sd)))
    return {"DownOut": A - B + C - D + F, "DownIn": B - C + D + E}


for rebate in (0.0, 2.0):
    for sigma in (0.005, 0.010, 0.012985838614515886, 0.016, 0.020):
        rf = ref(sigma, rebate)
        for name, bt in (("DownOut", ql.Barrier.DownOut), ("DownIn", ql.Barrier.DownIn)):
            a = ql_npv(sigma, bt, rebate)
            print(f"rebate {rebate} sigma {sigma:.6f} {name:7s} ql {a:.15g}  ref {mp.nstr(rf[name], 17)}  rel {float((a - rf[name]) / rf[name]):.2e}")

[boring-cyborg]

Thanks for opening this pull request! It might take a while before we look at it, so don't worry if there seems to be no feedback. We'll get to it.