feat(wing): analytic wing-seed for |k_log| in [2.95, 8.0]

Add an analytic wing-seed for the deep-wing regime, derived clean-room
from Schadner's IG-quantile and the 2-term Mills asymptotic. Replaces a
v1.0.0 path that extrapolated the Chebyshev seed beyond its fit domain
(SEED_K_HI = 3.0) and produced 3.27e-1 catastrophic errors on a
wing-saturated v x Delta stress grid.

Algorithm. For an OTM option with h = |k_log| and IG-survival
q_surv = (1 - c_*)/m:

  z1_0 = -Phi^-1(q_surv);   u_0 = -z1_0 + sqrt(z1_0^2 + 2h)   (W0 lead)

  for n in 0..N_PICARD:
    z2  = h/u + u/2
    Q   = 1/z2 - 1/z2^3                                       (Mills)
    phi = (2*pi)^(-1/2) * exp(-z1^2 / 2)
    delta = phi * Q                                           (= e^h * Phi(-z2), no overflow)
    z1  = -Phi^-1(q_surv + delta)
    u   = -z1 + sqrt(z1^2 + 2h)

  v = u

The collapse e^h * Phi(-z2) === phi(z1) * Q(z2) (with Q(z) :=
sqrt(pi/2) * erfcx(z/sqrt(2)) ~ 1/z - 1/z^3) eliminates the exp(h)
overflow on the wing. N_PICARD = 1 followed by HH3 polish lands at the
f64 conditioning floor.

Dispatch gates (chosen by Wren G's analysis):
  - h >= K_HI_BAILOUT = 2.95   (wing regime starts at |k_log| ~ 3)
  - q_surv < WING_Q_MAX = 0.30 (deep-OTM, the wing's natural support)
  - h < WING_H_MAX = 8.0       (above this the Mills truncation breaks
                                 down; route to Jaeckel rational fallback)

q-convention fix at the wing dispatch site. The IG kernel's q is the IG
CDF; the wing analytic uses IG survival. Without the conversion
`q_surv = 1 - q_kernel`, the wing seed receives the wrong probability
and the dispatch is silently wrong on a band of the wing grid. Caught
during integration, verifier-confirmed.

Cost on the wing v x Delta stress grid (360 cases, T = 1, F = 1, median
of 7 runs of 5000 reps each, single AVX-512 core, taskset -c 0):

  pre-wing  (v1.0.0):  Cheb extrapolation, 3.27e-1 max abs error
  post-wing (v1.0.1):  analytic seed,      8.30e-12 max abs error
  throughput:          81.3 ns/option (+8 ns vs Chebyshev path)

Schadner cold (1M synthetic options, the canonical grid) is unchanged:
73 ns / 3.42e-11 / 0 NaN. The wing path costs ~8 ns when it fires;
the Chebyshev path's 73 ns dominates the synthetic distribution.

Two pre-existing NaN on the wing v x Delta grid at (v=0.01, Delta in
{0.30, 0.70}): tiny-sigma near-ATM puts at the f64 BS price floor
(< 1e-7) where the f64 inverse is not meaningful. Pinned by
tests/wing_seed.rs::volfi_wing_grid_nan_set_bounded_to_two so a future
inner-iteration edit cannot silently expand the NaN set.

Files:
  - src/schadner_fast.rs:  +wing_seed_simd, +dispatch gates, clamps
  - src/otm_context.rs:    +dispatch in solve_with_ctx_simd
  - tests/wing_seed.rs:    +9 tests (Wren G corner, mpmath-200-bit
                            reference at h in {3..8} x q in {0.01..0.30},
                            boundary finiteness, SIMD lane independence,
                            kernel sigma recovery, Cheb non-regression,
                            context-API routing, NaN-set regression pin)
  - bench/wing_grid.rs:    +volfi-style v x Delta stress harness
  - scripts/wing_ref_gen.py:  +regenerates WING_REF from mpmath 200-bit
  - Cargo.toml:            +wing_grid binary entry

Author: RyanJamesStewart

Cargo.toml

@@ -55,3 +55,7 @@ lto = "thin"
 [[bin]]
 name = "kernel_nan"
 path = "bench/kernel_nan.rs"
+
+[[bin]]
+name = "wing_grid"
+path = "bench/wing_grid.rs"

bench/wing_grid.rs

@@ -0,0 +1,175 @@
+//! `bench/wing_grid.rs` — volfi-style v×Δ fixed-grid throughput benchmark.
+//!
+//! The volfi paper's evaluation grid:
+//!   v ∈ {0.01, 0.05, 0.10, 0.15, ..., 2.00}  (41 values)
+//!   Δ ∈ {0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.99} (13)
+//!   533 (v, Δ) points; repeated REPS times for a measurable timing window.
+//!
+//! Each (v, Δ) → (h, c_*) via:
+//!   d1 = Φ⁻¹(Δ);  d2 = d1 - v;  h = |k_log| = |d1·v - 0.5·v²| (need sign).
+//!   c_* = Φ(d1) - K/F · Φ(d2)  with K/F = exp(-h) for OTM call (k>0).
+//!
+//! We translate to voltic's (S, K, T, r, price) by fixing S=1, T=1, r=0,
+//! K = exp(k_log), price = c_* · S.
+//!
+//! Reports: ns/option for implied_vol_fast on the wing-dominated grid,
+//! and max abs err vs σ_true (= v / sqrt(T) = v with T=1).
+
+use std::time::Instant;
+use voltic::{implied_vol_fast, OptionKind};
+
+const V_GRID_START: f64 = 0.01;
+const V_GRID_STEP: f64 = 0.05;
+const V_GRID_N: usize = 41; // 0.01, 0.05, 0.10, ..., 2.00 (overstep is fine)
+const D_GRID: &[f64] = &[
+    0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.99,
+];
+const REPS: usize = 5000;
+
+fn phi_inv(p: f64) -> f64 {
+    // Acklam's algorithm, scalar; accurate enough for grid construction.
+    // Adapted from public domain Wichura AS241 polynomial.
+    let q = p - 0.5;
+    if q.abs() <= 0.425 {
+        let r = q * q;
+        q * ((((-39.69683028665376 * r + 220.9460984245205) * r - 275.9285104469687) * r
+            + 138.357751867269) * r - 30.66479806614716) * r + 2.506628277459239 / (((((-54.47609879822406 * r + 161.5858368580409) * r - 155.6989798598866) * r + 66.80131188771972) * r - 13.28068155288572) * r + 1.0)
+    } else {
+        // Tail. Use scipy-style approximation for ~6 digit accuracy; sufficient.
+        let r = if q < 0.0 { p } else { 1.0 - p };
+        let lr = (-r.ln()).sqrt();
+        let z = (((((2.938163982698783 * lr + 4.374664141464968) * lr - 2.549732539343734) * lr
+            - 2.400758277161838) * lr - 0.3223964580411365) * lr - 0.007784894002430293)
+            / ((((3.754408661907416 * lr + 2.445134137142996) * lr + 0.3224671290700398) * lr + 0.007784695709041462) * lr + 1.0);
+        if q < 0.0 { -z } else { z }
+    }
+}
+
+fn build_grid() -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<OptionKind>, Vec<f64>) {
+    let mut s = Vec::new();
+    let mut k = Vec::new();
+    let mut t = Vec::new();
+    let mut r = Vec::new();
+    let mut price = Vec::new();
+    let mut kind = Vec::new();
+    let mut sigma_true = Vec::new();
+    let spot: f64 = 1.0;
+    let tte: f64 = 1.0;
+    let rate: f64 = 0.0;
+    let sqrt_t: f64 = tte.sqrt();
+    for vi in 0..V_GRID_N {
+        let v = V_GRID_START + V_GRID_STEP * vi as f64;
+        if v < V_GRID_START || v > 2.0 {
+            continue;
+        }
+        for &dlt in D_GRID.iter() {
+            // d1 = Φ⁻¹(Δ); d2 = d1 − v.
+            // k_log = ln(K/F) = −(d1·v − 0.5·v²) = 0.5 v² − d1·v.
+            // For OTM call (positive k), Δ ∈ (0, 0.5): d1 < 0, k > 0.
+            let d1 = phi_inv(dlt);
+            let _d2 = d1 - v;
+            let k_log = 0.5 * v * v - d1 * v;
+            // We want OTM CALLS (positive k_log, Δ < 0.5) and OTM PUTS
+            // (negative k_log, |k_log| via Δ > 0.5). Both feed the wing.
+            let strike = (k_log).exp();
+            let sigma = v / sqrt_t;
+            // Use OTM-leg pricing: if k_log > 0, call; else put.
+            let is_call = k_log >= 0.0;
+            let opt_kind = if is_call { OptionKind::Call } else { OptionKind::Put };
+            // Price via voltic's bs_price so the round-trip is internally consistent.
+            let p = voltic::bs_price(&[spot], &[strike], &[tte], &[rate], &[sigma], &[opt_kind])[0];
+            // Filter: price must be above f64 noise so the inverse is meaningful.
+            if !(p.is_finite() && p > 1e-15) {
+                continue;
+            }
+            s.push(spot);
+            k.push(strike);
+            t.push(tte);
+            r.push(rate);
+            price.push(p);
+            kind.push(opt_kind);
+            sigma_true.push(sigma);
+        }
+    }
+    (s, k, t, r, price, kind, sigma_true)
+}
+
+fn main() {
+    let (s, k, t, r, price, kind, sigma_true) = build_grid();
+    let n = s.len();
+    eprintln!("built grid: {n} cases");
+
+    // Sanity: solve once, report max abs err.
+    let solved = implied_vol_fast(&s, &k, &t, &r, &price, &kind);
+    let mut max_abs = 0.0_f64;
+    let mut nan_count = 0;
+    for i in 0..n {
+        if solved[i].is_nan() {
+            nan_count += 1;
+            continue;
+        }
+        let e = (solved[i] - sigma_true[i]).abs();
+        if e > max_abs {
+            max_abs = e;
+        }
+    }
+    println!("=== volfi v×Δ grid benchmark ===");
+    println!("cases: {n}");
+    println!("max |σ_solved − σ_true|: {max_abs:.3e}");
+    println!("NaN count: {nan_count}");
+
+    // Warmup.
+    let w = implied_vol_fast(&s, &k, &t, &r, &price, &kind);
+    std::hint::black_box(w);
+
+    // Timed loop: REPS passes.
+    let mut samples = Vec::new();
+    for _ in 0..7 {
+        let t0 = Instant::now();
+        let mut total = 0usize;
+        for _ in 0..REPS {
+            let v = implied_vol_fast(&s, &k, &t, &r, &price, &kind);
+            total += v.len();
+            std::hint::black_box(&v);
+        }
+        let dt = t0.elapsed();
+        std::hint::black_box(total);
+        samples.push(dt.as_secs_f64() / (n as f64 * REPS as f64) * 1e9);
+    }
+    samples.sort_by(|a, b| a.partial_cmp(b).unwrap());
+    let median = samples[samples.len() / 2];
+    println!("median ns/option (median of 7, {} reps each): {median:.1}", REPS);
+    println!("options/sec: {:.3e}", 1e9 / median);
+
+    // Stratified report by whether the lane goes through the wing predicate.
+    let mut wing_lanes = 0usize;
+    let mut wing_max_err = 0.0_f64;
+    let mut cheb_lanes = 0usize;
+    let mut cheb_max_err = 0.0_f64;
+    for i in 0..n {
+        let k_log = (k[i] / s[i]).ln() - r[i] * t[i];
+        let h = k_log.abs();
+        let q_cdf = {
+            // Approximation: q = (1-c)/m; computed via the kernel formula.
+            let ek = k_log.exp();
+            let m = if k_log > 0.0 { 1.0 } else { ek };
+            let xn = price[i] / s[i];
+            let c = if matches!(kind[i], OptionKind::Call) { xn } else { xn + 1.0 - ek };
+            (1.0 - c) / m
+        };
+        let q_surv = 1.0 - q_cdf;
+        let in_wing = h >= 2.95 && h < 8.0 && q_surv > 0.0 && q_surv < 0.30;
+        if !solved[i].is_nan() {
+            let e = (solved[i] - sigma_true[i]).abs();
+            if in_wing {
+                wing_lanes += 1;
+                if e > wing_max_err { wing_max_err = e; }
+            } else {
+                cheb_lanes += 1;
+                if e > cheb_max_err { cheb_max_err = e; }
+            }
+        }
+    }
+    println!("wing lanes: {wing_lanes}  max abs err: {wing_max_err:.3e}");
+    println!("cheb lanes: {cheb_lanes}  max abs err: {cheb_max_err:.3e}");
+}

scripts/wing_ref_gen.py

@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+"""Regenerate the WING_REF table in tests/wing_seed.rs from mpmath at 200 bits.
+
+The wing-seed kernel maps (h = |k_log|, q = IG survival) to v = σ√T. The
+mpmath-200-bit reference table pins the seed against an oracle that has no
+dependency on the Rust kernel itself, so test agreement is independent
+evidence of correctness.
+
+Domain: h ∈ {3, 4, 5, 6, 7, 8}, q ∈ {0.01, 0.05, 0.10, 0.20, 0.30}.
+
+Reference equation: the IG-inverse-Gaussian survival relation that maps the
+wing parameter pair (h, q) back to v. Schadner's parameterisation is
+  q = F_IG(4/v² | 2/h, 1)  (the IG CDF with mean 2/h and shape 1, at 4/v²)
+which is equivalent to
+  q_surv = 1 - F_IG(4/v² | 2/h, 1).
+
+We invert by mp.findroot.
+
+Run:
+  python3 scripts/wing_ref_gen.py > /tmp/wing_ref.txt
+then paste the table into tests/wing_seed.rs above the existing WING_REF entries.
+
+Not wired into CI; here for reproducibility.
+"""
+import mpmath as mp
+
+mp.mp.prec = 200  # 200 bits ≈ 60 decimals
+
+
+def ig_cdf(x, mu, lam):
+    """Inverse-Gaussian CDF F(x | mu, lam) at 200 bits.
+
+        F(x) = Φ( √(λ/x) · (x/μ − 1) ) + exp(2λ/μ) · Φ( −√(λ/x) · (x/μ + 1) )
+
+    Φ via mp.ncdf.
+    """
+    a = mp.sqrt(lam / x) * (x / mu - 1)
+    b = -mp.sqrt(lam / x) * (x / mu + 1)
+    return mp.ncdf(a) + mp.exp(2 * lam / mu) * mp.ncdf(b)
+
+
+def solve_v_for(h, q_surv):
+    """Find v such that 1 - ig_cdf(4/v², 2/h, 1) = q_surv."""
+    mu = mp.mpf(2) / mp.mpf(h)
+    lam = mp.mpf(1)
+    target = mp.mpf(q_surv)
+
+    def f(v):
+        v = mp.fabs(v)
+        if v < mp.mpf("1e-12"):
+            v = mp.mpf("1e-12")
+        x = mp.mpf(4) / (v * v)
+        cdf = ig_cdf(x, mu, lam)
+        return (mp.mpf(1) - cdf) - target
+
+    # Initial guess: leading-order W0(h, q).
+    z1 = -mp.mpf(mp.qfunc(target).real if False else 0)  # placeholder; use a coarse seed
+    # Use a robust bracket-style: scan h-dependent range.
+    # In practice v ∈ (0.1, 10.0) covers the (h, q) cells in the table.
+    lo = mp.mpf("0.1")
+    hi = mp.mpf("10.0")
+    root = mp.findroot(f, (lo, hi), solver="anderson", tol=mp.mpf("1e-50"), maxsteps=200)
+    return mp.fabs(root)
+
+
+def main():
+    H_VALUES = [3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+    Q_VALUES = [0.01, 0.05, 0.10, 0.20, 0.30]
+    print("// Regenerated by scripts/wing_ref_gen.py at mpmath 200-bit precision.")
+    print("const WING_REF: &[(f64, f64, f64)] = &[")
+    for h in H_VALUES:
+        for q in Q_VALUES:
+            v = solve_v_for(h, q)
+            print(f"    ({h:.1f}, {q:.4f}, {float(v):.17e}),")
+    print("];")
+
+
+if __name__ == "__main__":
+    main()

src/otm_context.rs

@@ -39,7 +39,7 @@ use std::simd::StdFloat;
 use crate::schadner_fast::{
     cheb_qside_basis, cheb_qside_basis_simd, cheb_seed_from_basis_simd,
     cheb_seed_from_kside_basis_scalar, householder3_step_simd, ig_kt_prelude_scalar,
-    SEED_DEG_PUB,
+    wing_seed_simd, SEED_DEG_PUB, WING_H_MAX, WING_K_LO, WING_Q_MAX,
 };
 use crate::{LANES, M, V, VOL_MAX, VOL_MIN};
 
@@ -481,7 +481,33 @@ fn solve_with_ctx_simd(ctx: &OtmContextSimd, c_v: V) -> V {
     let q_for_seed = q_v.simd_max(p_lo_v).simd_min(p_hi_v);
 
     let cheb_tw = cheb_qside_basis_simd(q_for_seed);
-    let seed_v = cheb_seed_from_basis_simd(&ctx.cheb_tu, &cheb_tw);
+    let cheb_v = cheb_seed_from_basis_simd(&ctx.cheb_tu, &cheb_tw);
+
+    // Wing dispatch: lanes in the deep-OTM wing regime get the analytic
+    // wing seed instead of the Chebyshev seed. ctx.ak carries |k_log|. The
+    // wing seed expects IG SURVIVAL (= c_*); kernel q is the IG CDF
+    // (= 1 - c_*) — convert here.
+    let q_surv = V::splat(1.0) - q_v;
+    let use_wing = ctx.ak.simd_ge(V::splat(WING_K_LO))
+        & q_surv.simd_lt(V::splat(WING_Q_MAX))
+        & q_surv.simd_gt(V::splat(0.0))
+        & ctx.ak.simd_lt(V::splat(WING_H_MAX));
+    // Chunk-level bailout: skip the expensive wing_seed_simd call if no
+    // lane in this chunk needs it. Preserves cold-grid throughput.
+    let seed_v = if use_wing.any() {
+        let q_wing_clamped = q_surv
+            .simd_max(V::splat(1e-300))
+            .simd_min(V::splat(WING_Q_MAX));
+        let h_wing_clamped = ctx
+            .ak
+            .simd_max(V::splat(WING_K_LO))
+            .simd_min(V::splat(WING_H_MAX));
+        let wing_v = wing_seed_simd(h_wing_clamped, q_wing_clamped);
+        use_wing.select(wing_v, cheb_v)
+    } else {
+        cheb_v
+    };
+
     let mut v_iter = seed_v.simd_max(ctx.v_lo).simd_min(ctx.v_hi);
     let mut j = 0;
     while j < HOUSEHOLDER3_STEPS_CTX {