regit-svi 2.0.0

// Copyright 2026 Regit.io — Nicolas Koenig
// SPDX-License-Identifier: Apache-2.0

//! Static-arbitrage checks: butterfly (`g(k) >= 0`) and calendar-spread.
//!
//! # Butterfly arbitrage and the g function
//!
//! A slice admits butterfly arbitrage when the risk-neutral density it
//! implies is negative somewhere — a butterfly spread with negative cost.
//! Define
//!
//! ```text
//! g(k) = ( 1 - k*w'(k) / (2*w(k)) )^2
//!      - ( w'(k) / 2 )^2 * ( 1/w(k) + 1/4 )
//!      + w''(k) / 2
//! ```
//!
//! A regular Raw slice needs `g(k) >= 0` for all `k`, the strict right-call
//! boundary `b*(1 + rho) < 2`, and the non-strict left density-factor bound
//! `b*(1 - rho) <= 2`. Strictness matters: equality on the right does not make
//! the call price vanish at infinite strike.
//!
//! `g(k)` is evaluated in closed form from `w, w', w''`. [`assess_raw`] uses
//! the Martini–Mingone normalization, explicit root domains, and a recorded
//! numerical search on compactified half-wings. [`butterfly_scan`] remains a
//! deliberately bounded diagnostic over a caller-selected quoted range.
//!
//! # Calendar-spread arbitrage
//!
//! Two slices at maturities `t_1 < t_2` admit calendar-spread arbitrage when
//! their total-variance curves cross. Absence is the pointwise monotonicity
//! `w(k, t_1) <= w(k, t_2)` for all `k`. The difference
//! `D(k) = w(k, t_2) - w(k, t_1)` is scanned for negativity, with Brent
//! refinement of any crossing.
//!
//! For SSVI, [`crate::surface::ssvi`] exposes the analytic Theorem 4.1 calendar
//! characterization and the analytic sufficient Theorem 4.2 butterfly test.
//!
//! # References
//!
//! - Gatheral, J. & Jacquier, A., "Arbitrage-free SVI volatility surfaces",
//!   *Quantitative Finance* 14(1):59-71 (2014), Section 2.
//! - Roper, M., "Arbitrage free implied volatility surfaces", preprint,
//!   University of Sydney (2010).
//! - Lee, R. W., "The moment formula for implied volatility at extreme
//!   strikes", *Mathematical Finance* 14(3):469-480 (2004).

use crate::no_arb::evidence::{
    ArbitrageAssessment, ArbitrageEvidence, ArbitrageStatus, DiagnosticError, RawSearchStage,
    RootEvidence, RootTermination, ScanConfig, ScanEvidence, SearchEvidence,
};
use crate::numerics::{brent_root_with_evidence, index_to_f64};
use crate::smile::raw::RawSvi;

/// Number of grid points used by the butterfly and calendar scans.
const SCAN_POINTS: usize = 401;
/// Wing margin added on each side of the quoted range when scanning.
const SCAN_MARGIN: f64 = 1.0;
/// Brent tolerance for refining a reported arbitrage boundary.
const REFINE_TOL: f64 = 1e-10;
/// Brent iteration cap for boundary refinement.
const REFINE_MAX_ITER: usize = 200;

/// The butterfly function `g(k)` for a raw SVI slice (MATH.md §7).
///
/// `g` is the sign-controlling factor in the risk-neutral density:
/// `g(k) >= 0` everywhere, together with the relevant call-price boundary,
/// characterizes butterfly freedom for a regular positive slice.
/// This infallible kernel assumes finite `k` and positive finite `w(k)` and
/// propagates IEEE non-finite results when those preconditions fail.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::no_arb::butterfly::g;
///
/// // A gently curved slice has positive g near the money.
/// let svi = RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3)?;
/// assert!(g(&svi, 0.0) > 0.0);
/// # Ok(())
/// # }
/// ```
#[must_use]
pub fn g(svi: &RawSvi, k: f64) -> f64 {
    g_with_scale(svi, k).0
}

/// Evaluates `g` and a conservative scale for cancellation in its three
/// signed terms.
pub(crate) fn g_with_scale(svi: &RawSvi, k: f64) -> (f64, f64) {
    let w = svi.total_variance(k);
    let wp = svi.w_prime(k);
    let wpp = svi.w_double_prime(k);

    // term1 = (1 - k*w'/(2*w))^2
    let t1 = {
        let inner = 1.0 - k * wp / (2.0 * w);
        inner * inner
    };
    // term2 = (w'/2)^2 * (1/w + 1/4)
    let t2 = {
        let half_wp = wp / 2.0;
        half_wp * half_wp * (1.0 / w + 0.25)
    };
    // term3 = w''/2
    let t3 = wpp / 2.0;

    (t1 - t2 + t3, 1.0 + t1.abs() + t2.abs() + t3.abs())
}

/// The result of a butterfly-arbitrage scan over a raw SVI slice.
///
/// A clean diagnostic says only that no violation was observed on its grid.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ButterflyDiagnostic {
    /// Whether a density factor below the recorded negative tolerance was sampled.
    violation_observed: bool,
    /// The smallest value of `g` observed on the scan grid.
    min_g: f64,
    /// The log-moneyness at which [`Self::min_g`] was observed.
    worst_k: f64,
    /// Refined zero boundary near the sampled violation, when bracketed.
    refined_boundary: Option<f64>,
    /// Reproducible evidence for the bounded scan.
    evidence: ScanEvidence,
}

impl ButterflyDiagnostic {
    /// Returns whether a density factor below the recorded negative tolerance was observed.
    #[must_use]
    pub const fn violation_observed(self) -> bool {
        self.violation_observed
    }
    /// Returns the minimum sampled density factor.
    #[must_use]
    pub const fn min_g(self) -> f64 {
        self.min_g
    }
    /// Returns the worst sampled log-moneyness.
    #[must_use]
    pub const fn worst_k(self) -> f64 {
        self.worst_k
    }
    /// Returns a nearby refined `g(k)=0` boundary when one was bracketed.
    #[must_use]
    pub const fn refined_boundary(self) -> Option<f64> {
        self.refined_boundary
    }
    /// Returns the bounded numerical evidence.
    #[must_use]
    pub const fn evidence(self) -> ScanEvidence {
        self.evidence
    }
}

/// Tests the Raw SVI call/density tail bounds (Lee 2004).
///
/// This checks `b*(1 + rho) < 2` on the right and `b*(1 - rho) <= 2` on the
/// left. It is a necessary tail check, not a test of the interior `g` factor.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::no_arb::butterfly::wing_bound_ok;
///
/// let ok = RawSvi::new(0.04, 0.5, -0.3, 0.0, 0.1)?;
/// assert!(wing_bound_ok(&ok));
/// let steep = RawSvi::new(0.04, 3.0, -0.3, 0.0, 0.1)?;
/// assert!(!wing_bound_ok(&steep));
/// # Ok(())
/// # }
/// ```
#[must_use]
#[inline]
pub fn wing_bound_ok(svi: &RawSvi) -> bool {
    svi.b * (1.0 + svi.rho) < 2.0 && svi.b * (1.0 - svi.rho) <= 2.0
}

/// Scans a raw SVI slice for butterfly arbitrage over `[k_lo, k_hi]` plus a
/// wing margin (MATH.md §7).
///
/// Evaluates `g` on a dense grid. A value below the scale-aware negative
/// tolerance recorded in the scan evidence flags arbitrage; values inside the
/// band are unresolved floating-point boundary observations. When an adjacent
/// pair brackets a sign change, [`ButterflyDiagnostic::refined_boundary`]
/// reports a Brent-refined boundary selected nearest the deepest sampled
/// violation. [`ButterflyDiagnostic::worst_k`] remains a sampled witness.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::no_arb::butterfly::butterfly_scan;
///
/// // A well-behaved slice is free of butterfly arbitrage.
/// let svi = RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3)?;
/// let report = butterfly_scan(&svi, -0.5, 0.5)?;
/// assert!(!report.violation_observed());
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`DiagnosticError`] if either bound is non-finite, `k_lo >= k_hi`,
/// or expanding the interval by the built-in margin overflows.
pub fn butterfly_scan(
    svi: &RawSvi,
    k_lo: f64,
    k_hi: f64,
) -> Result<ButterflyDiagnostic, DiagnosticError> {
    if !k_lo.is_finite() || !k_hi.is_finite() {
        return Err(DiagnosticError::NonFiniteBound);
    }
    if k_lo >= k_hi {
        return Err(DiagnosticError::InvalidOrder);
    }
    let lo = k_lo - SCAN_MARGIN;
    let hi = k_hi + SCAN_MARGIN;
    if ScanConfig::new(lo, hi, SCAN_POINTS, 0.0).is_none() {
        return Err(DiagnosticError::DomainOverflow);
    }
    let step = (hi - lo) / index_to_f64(SCAN_POINTS - 1);

    let mut min_g = f64::INFINITY;
    let mut max_evaluation_scale = 1.0_f64;
    let mut worst_k = lo;
    let mut samples = Vec::with_capacity(SCAN_POINTS);

    for i in 0..SCAN_POINTS {
        let k = step.mul_add(index_to_f64(i), lo);
        let (gi, evaluation_scale) = g_with_scale(svi, k);
        if !gi.is_finite() {
            return Err(DiagnosticError::NonFiniteEvaluation);
        }
        max_evaluation_scale = max_evaluation_scale.max(evaluation_scale);
        if gi < min_g {
            min_g = gi;
            worst_k = k;
        }
        samples.push((k, gi));
    }

    let evaluation_tolerance = 128.0 * f64::EPSILON * max_evaluation_scale;
    let Some(config) = ScanConfig::new(lo, hi, SCAN_POINTS, evaluation_tolerance) else {
        return Err(DiagnosticError::DomainOverflow);
    };
    let violation_observed = min_g < -evaluation_tolerance;

    let refinement_bracket = violation_observed
        .then(|| {
            samples
                .windows(2)
                .filter(|pair| {
                    pair[0].1 == 0.0
                        || pair[1].1 == 0.0
                        || pair[0].1.is_sign_negative() != pair[1].1.is_sign_negative()
                })
                .min_by(|left, right| {
                    let left_distance = (0.5 * (left[0].0 + left[1].0) - worst_k).abs();
                    let right_distance = (0.5 * (right[0].0 + right[1].0) - worst_k).abs();
                    left_distance.total_cmp(&right_distance)
                })
                .map(|pair| (pair[0].0, pair[1].0))
        })
        .flatten();
    let refinement = refinement_bracket.and_then(|(lower, upper)| {
        brent_root_with_evidence(|x| g(svi, x), lower, upper, REFINE_TOL, REFINE_MAX_ITER)
    });
    let refined_boundary = refinement.map(crate::numerics::BrentRoot::root);
    let root_evidence = refinement.map(|root| {
        RootEvidence::new(
            root.root(),
            root.lower(),
            root.upper(),
            root.residual(),
            root.evaluations(),
            if root.exact() {
                RootTermination::ExactRoot
            } else {
                RootTermination::BracketTolerance
            },
        )
    });

    Ok(ButterflyDiagnostic {
        violation_observed,
        min_g,
        worst_k,
        refined_boundary,
        evidence: ScanEvidence::new(
            config,
            (k_lo, k_hi),
            SCAN_POINTS,
            SCAN_POINTS,
            refinement_bracket.is_some(),
            root_evidence,
        ),
    })
}

#[derive(Clone, Copy)]
struct NormalizedRaw {
    alpha: f64,
    b: f64,
    rho: f64,
    mu: f64,
}

impl NormalizedRaw {
    fn values(self, ell: f64) -> Option<(f64, f64, f64, f64, f64)> {
        let radius = ell.hypot(1.0);
        let n = self.b.mul_add(self.rho.mul_add(ell, radius), self.alpha);
        let np = self.b * (self.rho + ell / radius);
        let npp = self.b / (radius * radius * radius);
        let g1_plus = 1.0 - np * ((ell + self.mu) / (2.0 * n) + 0.25);
        let g1_minus = 1.0 - np * ((ell + self.mu) / (2.0 * n) - 0.25);
        let g1 = g1_plus * g1_minus;
        let g2 = npp - np * np / (2.0 * n);
        [n, np, npp, g1, g2]
            .iter()
            .all(|x| x.is_finite())
            .then_some((n, np, npp, g1, g2))
    }

    fn h_plus(self, ell: f64) -> Option<f64> {
        let (n, np, npp, _, _) = self.values(ell)?;
        (np * np * (1.0 - np / 2.0) - 2.0 * n * npp)
            .is_finite()
            .then_some(np * np * (1.0 - np / 2.0) - 2.0 * n * npp)
    }

    fn h_minus(self, ell: f64) -> Option<f64> {
        let (n, np, npp, _, _) = self.values(ell)?;
        (np * np * (1.0 + np / 2.0) - 2.0 * n * npp)
            .is_finite()
            .then_some(np * np * (1.0 + np / 2.0) - 2.0 * n * npp)
    }

    fn g2(self, ell: f64) -> Option<f64> {
        self.values(ell).map(|(_, _, _, _, g2)| g2)
    }

    fn objective(self, ell: f64) -> Option<f64> {
        let (_, _, _, g1, g2) = self.values(ell)?;
        let value = -g2 / (2.0 * g1);
        value.is_finite().then_some(value.max(0.0))
    }
}

#[derive(Clone, Copy)]
enum WingBranch {
    B1,
    B2,
    B3,
    B4,
}

fn scale_band(abs_tol: f64, values: &[f64]) -> f64 {
    const REL_TOL: f64 = 64.0 * f64::EPSILON;
    let scale = values
        .iter()
        .fold(1.0_f64, |acc, value| acc.max(value.abs()));
    abs_tol + REL_TOL * scale
}

fn root_with_evidence(
    function: impl Fn(f64) -> Option<f64>,
    mut lower: f64,
    mut upper: f64,
    tolerance: f64,
) -> Option<(f64, RootEvidence)> {
    let mut f_lower = function(lower)?;
    let f_upper = function(upper)?;
    let mut evaluations = 2;
    if f_lower == 0.0 {
        return Some((
            lower,
            RootEvidence::new(
                lower,
                lower,
                lower,
                0.0,
                evaluations,
                RootTermination::ExactRoot,
            ),
        ));
    }
    if f_upper == 0.0 {
        return Some((
            upper,
            RootEvidence::new(
                upper,
                upper,
                upper,
                0.0,
                evaluations,
                RootTermination::ExactRoot,
            ),
        ));
    }
    if f_lower.is_sign_positive() == f_upper.is_sign_positive() {
        return None;
    }
    let mut midpoint = lower + (upper - lower) / 2.0;
    let mut f_midpoint = function(midpoint)?;
    evaluations += 1;
    for _ in 0..160 {
        if (upper - lower).abs() <= tolerance * (1.0 + midpoint.abs()) || f_midpoint == 0.0 {
            break;
        }
        if f_lower.is_sign_positive() == f_midpoint.is_sign_positive() {
            lower = midpoint;
            f_lower = f_midpoint;
        } else {
            upper = midpoint;
        }
        midpoint = lower + (upper - lower) / 2.0;
        f_midpoint = function(midpoint)?;
        evaluations += 1;
    }
    let termination = if f_midpoint == 0.0 {
        RootTermination::ExactRoot
    } else if (upper - lower).abs() <= tolerance * (1.0 + midpoint.abs()) {
        RootTermination::BracketTolerance
    } else {
        return None;
    };
    Some((
        midpoint,
        RootEvidence::new(
            midpoint,
            lower,
            upper,
            f_midpoint.abs(),
            evaluations,
            termination,
        ),
    ))
}

fn bracket_from_stationary(
    normalized: NormalizedRaw,
    left: bool,
    tolerance: f64,
) -> Option<(f64, RootEvidence)> {
    let ell_star = -normalized.rho / (1.0 - normalized.rho * normalized.rho).sqrt();
    let direction = if left { -1.0 } else { 1.0 };
    let function = |ell| {
        if left {
            normalized.h_minus(ell)
        } else {
            normalized.h_plus(ell)
        }
    };
    let mut distance = (1.0 + ell_star.abs()) * 1e-6;
    let inner = ell_star + direction * distance;
    let inner_value = function(inner)?;
    if inner_value >= 0.0 {
        return None;
    }
    for _ in 0..160 {
        distance *= 2.0;
        let outer = ell_star + direction * distance;
        let outer_value = function(outer)?;
        if outer_value > 0.0 {
            return if left {
                root_with_evidence(function, outer, inner, tolerance)
            } else {
                root_with_evidence(function, inner, outer, tolerance)
            };
        }
    }
    None
}

fn bracket_g2(
    normalized: NormalizedRaw,
    left: bool,
    tolerance: f64,
) -> Option<(f64, RootEvidence)> {
    let direction = if left { -1.0 } else { 1.0 };
    let at_zero = normalized.g2(0.0)?;
    if at_zero <= 0.0 {
        return None;
    }
    let mut distance = 1.0;
    for _ in 0..160 {
        let outer = direction * distance;
        let outer_value = normalized.g2(outer)?;
        if outer_value < 0.0 {
            return if left {
                root_with_evidence(|ell| normalized.g2(ell), outer, 0.0, tolerance)
            } else {
                root_with_evidence(|ell| normalized.g2(ell), 0.0, outer, tolerance)
            };
        }
        distance *= 2.0;
    }
    None
}

#[derive(Clone, Copy)]
struct SearchState {
    roots: [Option<RootEvidence>; 4],
    interval: Option<(f64, f64)>,
    domains: [(f64, f64); 2],
    evaluations: usize,
}

#[allow(clippy::too_many_arguments)] // Mirrors the immutable public audit record field-for-field.
fn search_evidence(
    stage: RawSearchStage,
    tolerance: f64,
    trace: SearchState,
    optimization_error: f64,
    subdivisions: usize,
    sigma_star: Option<f64>,
    argmax_ell: Option<f64>,
    terminated: bool,
) -> ArbitrageEvidence {
    ArbitrageEvidence::NumericalSearch(SearchEvidence::new(
        "Martini–Mingone roots with global dyadic compactification",
        stage,
        trace.domains,
        tolerance,
        64.0 * f64::EPSILON,
        optimization_error,
        trace.evaluations,
        subdivisions,
        trace.roots,
        trace.interval,
        sigma_star,
        argmax_ell,
        terminated,
    ))
}

#[allow(clippy::type_complexity)] // Private fixed tuple avoids another one-use state type.
fn compactified_search(
    normalized: NormalizedRaw,
    left_root: f64,
    right_root: f64,
) -> Option<(f64, f64, f64, usize, usize, [(f64, f64); 2])> {
    const LEVEL: usize = 15;
    let intervals = 1_usize << LEVEL;
    let domains = [(1.0 / left_root, 0.0), (0.0, 1.0 / right_root)];
    let mut best = 0.0_f64;
    let mut argmax = left_root;
    let mut previous_best = 0.0_f64;
    let mut evaluations = 0;
    for level in 8..=LEVEL {
        let count = 1_usize << level;
        let mut level_best = 0.0_f64;
        for &(lower, upper) in &domains {
            for i in 0..=count {
                let h = (upper - lower).mul_add(index_to_f64(i) / index_to_f64(count), lower);
                let value = if i == 0 && lower == 0.0 || i == count && upper == 0.0 {
                    0.0
                } else {
                    normalized.objective(1.0 / h)?
                };
                evaluations += 1;
                if value > level_best {
                    level_best = value;
                }
                if value > best {
                    best = value;
                    argmax = 1.0 / h;
                }
            }
        }
        if level < LEVEL {
            previous_best = level_best;
        }
    }
    let max_width = domains
        .iter()
        .map(|(lo, hi)| (hi - lo).abs())
        .fold(0.0_f64, f64::max)
        / index_to_f64(intervals);
    let optimization_error = (best - previous_best).abs() + max_width * (1.0 + best.abs());
    Some((
        best,
        argmax,
        optimization_error,
        evaluations,
        intervals * 2,
        domains,
    ))
}

/// Runs a theorem-guided Raw SVI assessment.
///
/// Flat and isolated-zero cases and failed necessary conditions carry
/// analytic evidence. The regular interior uses the explicit
/// Martini–Mingone reductions followed by bracket-preserving root solves and
/// a global dyadic search on compactified half-wings. It therefore never carries
/// [`ArbitrageEvidence::AnalyticNecessaryAndSufficient`] evidence.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::no_arb::butterfly::assess_raw;
/// use regit_svi::ArbitrageStatus;
///
/// let svi = RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3)?;
/// assert_eq!(assess_raw(&svi, 1e-10).status(), ArbitrageStatus::NoViolationDetected);
/// # Ok(())
/// # }
/// ```
#[must_use]
#[allow(clippy::too_many_lines)] // Kept linear so the published decision stages are auditable in order.
#[allow(clippy::float_cmp)] // Exact equality selects the paper's equality branch and strict right-tail failure.
pub fn assess_raw(svi: &RawSvi, boundary_tolerance: f64) -> ArbitrageAssessment {
    if !boundary_tolerance.is_finite() || boundary_tolerance < 0.0 {
        let state = SearchState {
            roots: [None; 4],
            interval: None,
            domains: [(f64::NAN, f64::NAN); 2],
            evaluations: 0,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::Domain,
                0.0,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    }
    if svi.b == 0.0 {
        let evidence = ArbitrageEvidence::AnalyticNecessaryAndSufficient {
            theorem: "constant positive total-variance special case",
            boundary_tolerance,
        };
        let status = if svi.a > 0.0 {
            ArbitrageStatus::NoViolationDetected
        } else {
            ArbitrageStatus::Indeterminate
        };
        return ArbitrageAssessment::new(status, evidence, svi.a, None);
    }
    if svi.w_min() == 0.0 {
        let evidence = ArbitrageEvidence::AnalyticNecessaryAndSufficient {
            theorem: "Martini–Mingone isolated-zero result",
            boundary_tolerance,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::ViolationDetected,
            evidence,
            -0.0,
            Some(svi.k_min()),
        );
    }
    let alpha = svi.a / svi.sigma;
    let mu = svi.m / svi.sigma;
    if [alpha, mu, svi.b, svi.rho]
        .iter()
        .any(|value| !value.is_finite() || value.is_subnormal())
    {
        let state = SearchState {
            roots: [None; 4],
            interval: None,
            domains: [(f64::NAN, f64::NAN); 2],
            evaluations: 0,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::Domain,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    }
    let normalized = NormalizedRaw {
        alpha,
        b: svi.b,
        rho: svi.rho,
        mu,
    };
    let right_slope = svi.b * (1.0 + svi.rho);
    let left_slope = svi.b * (1.0 - svi.rho);
    let right_margin = 2.0 - right_slope;
    let left_margin = 2.0 - left_slope;
    let tail_margin = right_margin.min(left_margin);
    let wing_band = scale_band(boundary_tolerance, &[right_slope, left_slope, 2.0]);
    if right_slope == 2.0 || right_margin < -wing_band || left_margin < -wing_band {
        let evidence = ArbitrageEvidence::AnalyticNecessary {
            condition: "strict right-call tail and non-strict density-factor wing bounds",
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::ViolationDetected,
            evidence,
            tail_margin,
            None,
        );
    }
    let near_right = right_margin <= wing_band;
    let near_left = left_margin.abs() <= wing_band && left_slope != 2.0;
    if near_right || near_left {
        let state = SearchState {
            roots: [None; 4],
            interval: None,
            domains: [(f64::NAN, f64::NAN); 2],
            evaluations: 0,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::Wings,
                boundary_tolerance,
                state,
                wing_band,
                0,
                None,
                None,
                false,
            ),
            tail_margin,
            None,
        );
    }

    let branch = match (left_slope == 2.0, right_slope == 2.0) {
        (false, false) => WingBranch::B1,
        (true, false) => WingBranch::B2,
        (false, true) => WingBranch::B3,
        (true, true) => WingBranch::B4,
    };
    let root_tolerance = (boundary_tolerance * 0.1).clamp(1e-13, 1e-10);
    let left_fukasawa = if matches!(branch, WingBranch::B1 | WingBranch::B3) {
        bracket_from_stationary(normalized, true, root_tolerance)
    } else {
        None
    };
    let right_fukasawa = if matches!(branch, WingBranch::B1 | WingBranch::B2) {
        bracket_from_stationary(normalized, false, root_tolerance)
    } else {
        None
    };
    if matches!(branch, WingBranch::B1 | WingBranch::B3) && left_fukasawa.is_none()
        || matches!(branch, WingBranch::B1 | WingBranch::B2) && right_fukasawa.is_none()
    {
        let state = SearchState {
            roots: [
                left_fukasawa.map(|x| x.1),
                right_fukasawa.map(|x| x.1),
                None,
                None,
            ],
            interval: None,
            domains: [(f64::NAN, f64::NAN); 2],
            evaluations: 0,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::FukasawaInterval,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    }
    let lower = if let Some((ell, _)) = left_fukasawa {
        normalized
            .values(ell)
            .map(|(n, np, _, _, _)| 2.0 * n * (1.0 / np + 0.25) - ell)
    } else {
        Some(-alpha / 2.0)
    };
    let upper = if let Some((ell, _)) = right_fukasawa {
        normalized
            .values(ell)
            .map(|(n, np, _, _, _)| 2.0 * n * (1.0 / np - 0.25) - ell)
    } else {
        Some(alpha / 2.0)
    };
    let (Some(lower), Some(upper)) = (lower, upper) else {
        let state = SearchState {
            roots: [
                left_fukasawa.map(|x| x.1),
                right_fukasawa.map(|x| x.1),
                None,
                None,
            ],
            interval: None,
            domains: [(f64::NAN, f64::NAN); 2],
            evaluations: 0,
        };
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::FukasawaInterval,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    };
    let interval = (lower, upper);
    let root_evaluations = left_fukasawa.map_or(0, |x| x.1.evaluations())
        + right_fukasawa.map_or(0, |x| x.1.evaluations());
    let mut state = SearchState {
        roots: [
            left_fukasawa.map(|x| x.1),
            right_fukasawa.map(|x| x.1),
            None,
            None,
        ],
        interval: Some(interval),
        domains: [(f64::NAN, f64::NAN); 2],
        evaluations: root_evaluations,
    };
    let interval_band = scale_band(boundary_tolerance, &[lower, upper, mu]);
    if !lower.is_finite() || !upper.is_finite() {
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::FukasawaInterval,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    }
    let interval_margin = upper - lower;
    if interval_margin < -interval_band {
        return ArbitrageAssessment::new(
            ArbitrageStatus::ViolationDetected,
            search_evidence(
                RawSearchStage::FukasawaInterval,
                boundary_tolerance,
                state,
                interval_band,
                0,
                None,
                None,
                true,
            ),
            interval_margin,
            None,
        );
    }
    if interval_margin <= interval_band {
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::FukasawaInterval,
                boundary_tolerance,
                state,
                interval_band,
                0,
                None,
                None,
                true,
            ),
            interval_margin,
            None,
        );
    }
    let location_margin = (mu - lower).min(upper - mu);
    if location_margin < -interval_band {
        return ArbitrageAssessment::new(
            ArbitrageStatus::ViolationDetected,
            search_evidence(
                RawSearchStage::Location,
                boundary_tolerance,
                state,
                interval_band,
                0,
                None,
                None,
                true,
            ),
            location_margin,
            None,
        );
    }
    if location_margin <= interval_band {
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::Location,
                boundary_tolerance,
                state,
                interval_band,
                0,
                None,
                None,
                true,
            ),
            location_margin,
            None,
        );
    }

    let left_g2 = bracket_g2(normalized, true, root_tolerance);
    let right_g2 = bracket_g2(normalized, false, root_tolerance);
    let (Some((left_g2_root, left_g2_evidence)), Some((right_g2_root, right_g2_evidence))) =
        (left_g2, right_g2)
    else {
        state.roots[2] = left_g2.map(|x| x.1);
        state.roots[3] = right_g2.map(|x| x.1);
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::G2Roots,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    };
    state.roots[2] = Some(left_g2_evidence);
    state.roots[3] = Some(right_g2_evidence);
    state.evaluations += left_g2_evidence.evaluations() + right_g2_evidence.evaluations();
    let Some((sigma_star, argmax_ell, optimization_error, evaluations, subdivisions, domains)) =
        compactified_search(normalized, left_g2_root, right_g2_root)
    else {
        return ArbitrageAssessment::new(
            ArbitrageStatus::Indeterminate,
            search_evidence(
                RawSearchStage::SigmaSearch,
                boundary_tolerance,
                state,
                f64::NAN,
                0,
                None,
                None,
                false,
            ),
            f64::NAN,
            None,
        );
    };
    state.domains = domains;
    state.evaluations += evaluations;
    let margin = svi.sigma - sigma_star;
    let uncertainty = scale_band(boundary_tolerance, &[svi.sigma, sigma_star]) + optimization_error;
    let witness_k = svi.sigma.mul_add(argmax_ell, svi.m);
    let witness_g = g(svi, witness_k);
    let evidence = search_evidence(
        RawSearchStage::Complete,
        boundary_tolerance,
        state,
        optimization_error,
        subdivisions,
        Some(sigma_star),
        Some(argmax_ell),
        true,
    );
    let status =
        if margin < -uncertainty && witness_g < -scale_band(boundary_tolerance, &[witness_g]) {
            ArbitrageStatus::ViolationDetected
        } else if margin > uncertainty && witness_g.is_finite() {
            ArbitrageStatus::NoViolationDetected
        } else {
            ArbitrageStatus::Indeterminate
        };
    ArbitrageAssessment::new(status, evidence, margin.min(tail_margin), Some(witness_k))
}

#[cfg(test)]
#[allow(clippy::expect_used)] // Validated fixtures use contextual expectations.
mod tests {
    use super::*;
    use crate::surface::ssvi::{Phi, Ssvi};

    #[test]
    fn g_is_positive_for_benign_slice() {
        let svi =
            RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3).expect("valid test or documentation fixture");
        for &k in &[-1.0, -0.3, 0.0, 0.3, 1.0] {
            assert!(g(&svi, k) > 0.0, "g({k}) should be positive");
        }
    }

    #[test]
    fn g_equals_density_factor_at_atm() {
        // For a flat slice (b = 0), w' = w'' = 0, so g(0) = 1.
        let flat =
            RawSvi::new(0.04, 0.0, 0.0, 0.0, 0.1).expect("valid test or documentation fixture");
        assert!((g(&flat, 0.0) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn wing_bound_accepts_gentle_rejects_steep() {
        assert!(wing_bound_ok(
            &RawSvi::new(0.04, 0.5, -0.3, 0.0, 0.1).expect("valid test or documentation fixture")
        ));
        assert!(!wing_bound_ok(
            &RawSvi::new(0.04, 3.0, -0.3, 0.0, 0.1).expect("valid test or documentation fixture")
        ));
    }

    #[test]
    fn butterfly_scan_passes_benign_slice() {
        let svi =
            RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3).expect("valid test or documentation fixture");
        let report = butterfly_scan(&svi, -0.5, 0.5).expect("valid test or documentation fixture");
        assert!(!report.violation_observed);
        assert!(report.min_g > 0.0);
    }

    #[test]
    fn butterfly_scan_rejects_unordered_bounds() {
        let slice = RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3).expect("valid test fixture");
        assert_eq!(
            butterfly_scan(&slice, 1.0, -1.0),
            Err(DiagnosticError::InvalidOrder)
        );
        assert_eq!(
            butterfly_scan(&slice, 0.0, 0.0),
            Err(DiagnosticError::InvalidOrder)
        );
    }

    #[test]
    fn butterfly_scan_flags_vogt_slice() {
        // The Axel Vogt slice from Gatheral & Jacquier (2014), Section 2.2 —
        // a raw SVI slice with documented butterfly arbitrage.
        // a = -0.0410, b = 0.1331, rho = 0.3060, m = 0.3586, sigma = 0.4153.
        let vogt = RawSvi::new(-0.0410, 0.1331, 0.3060, 0.3586, 0.4153)
            .expect("valid test or documentation fixture");
        let report = butterfly_scan(&vogt, -1.5, 1.5).expect("valid test or documentation fixture");
        assert!(
            report.violation_observed,
            "Vogt slice must be flagged as arbitrageable"
        );
        assert!(report.refined_boundary().is_some());
        let evidence = report.evidence();
        assert_eq!(evidence.requested_domain(), (-1.5, 1.5));
        assert_eq!(evidence.requested_points(), SCAN_POINTS);
        assert!(evidence.config().lower() < evidence.requested_domain().0);
        assert!(evidence.refinement_attempted());
        let refinement = evidence
            .refinement()
            .expect("completed sign-change refinement");
        assert!(
            (refinement.root() - report.refined_boundary().expect("refined root")).abs() < 1e-15
        );
        assert!(refinement.residual() <= REFINE_TOL);
        assert!(refinement.evaluations() >= 2);
        assert!(matches!(
            refinement.termination(),
            RootTermination::ExactRoot | RootTermination::BracketTolerance
        ));
        assert!(report.min_g < 0.0, "min_g = {}", report.min_g);
    }

    #[test]
    fn strict_right_tail_equality_is_an_analytic_violation() {
        let sigma = 3.0 * 3.0_f64.sqrt() / 4.0;
        let slice = RawSvi::new(1.5, 4.0 / 3.0, 0.5, -0.75, sigma)
            .expect("valid test or documentation fixture");
        assert!(g(&slice, 10.0) > 0.0);
        let assessment = assess_raw(&slice, 1e-12);
        assert_eq!(assessment.status(), ArbitrageStatus::ViolationDetected);
        assert!(matches!(
            assessment.evidence(),
            ArbitrageEvidence::AnalyticNecessary { .. }
        ));
    }

    #[test]
    fn fukasawa_interval_roundoff_boundary_is_indeterminate() {
        let alpha = -0.499_567_894_485_558_65;
        let slice = RawSvi::new(alpha * 0.1, 0.5, 0.0, 0.0, 0.1).expect("valid boundary fixture");
        let assessment = assess_raw(&slice, 1e-10);
        assert_eq!(assessment.status(), ArbitrageStatus::Indeterminate);
        assert!(assessment.margin().abs() <= 1e-10);
        assert!(matches!(
            assessment.evidence(),
            ArbitrageEvidence::NumericalSearch(search)
                if search.stage() == RawSearchStage::FukasawaInterval
        ));
    }

    #[test]
    fn isolated_zero_is_not_hidden_by_nan_density_factor() {
        let slice =
            RawSvi::new(-0.125, 0.5, 0.0, 0.0, 0.25).expect("valid test or documentation fixture");
        assert!(slice.w_min().abs() <= f64::EPSILON);
        assert_eq!(
            assess_raw(&slice, 1e-12).status(),
            ArbitrageStatus::ViolationDetected
        );
    }

    #[test]
    fn martini_mingone_symmetric_anchor_matches_reference() {
        let sigma = 0.1;
        let slice = RawSvi::new(0.1 * sigma, 0.5, 0.0, 0.0, sigma)
            .expect("valid test or documentation fixture");
        let assessment = assess_raw(&slice, 1e-12);
        let evidence = match assessment.evidence() {
            ArbitrageEvidence::NumericalSearch(evidence) => evidence,
            other => {
                assert!(matches!(other, ArbitrageEvidence::NumericalSearch(_)));
                return;
            }
        };
        let (lower, upper) = evidence
            .fukasawa_interval()
            .expect("valid test or documentation fixture");
        assert!(
            (lower + 2.797_192_337_566_493_6).abs() < 2e-9,
            "lower={lower}"
        );
        assert!(
            (upper - 2.797_192_337_566_493_6).abs() < 2e-9,
            "upper={upper}"
        );
        let roots = evidence.roots();
        let left_g2 = roots[2].expect("valid test or documentation fixture");
        let right_g2 = roots[3].expect("valid test or documentation fixture");
        assert!((left_g2.lower() + 1.490_910_442_675_791_7).abs() < 2e-9);
        assert!((right_g2.upper() - 1.490_910_442_675_791_7).abs() < 2e-9);
        let sigma_star = evidence
            .sigma_star()
            .expect("valid test or documentation fixture");
        assert!(
            (sigma_star - 0.082_735_513_672_510_96).abs() < 2e-8,
            "sigma_star={sigma_star}"
        );
        assert!(
            (evidence
                .argmax_ell()
                .expect("valid test or documentation fixture")
                .abs()
                - 3.562_806_352_647_235)
                .abs()
                < 2e-4
        );
        assert_eq!(evidence.stage(), RawSearchStage::Complete);
        assert!(evidence.terminated());
        assert_eq!(assessment.status(), ArbitrageStatus::NoViolationDetected);
    }

    #[test]
    fn martini_mingone_pass_violation_and_boundary_regressions() {
        let pass =
            RawSvi::new(0.01, 0.5, 0.0, 0.0, 0.1).expect("valid test or documentation fixture");
        let fail =
            RawSvi::new(0.004, 0.5, 0.0, 0.0, 0.04).expect("valid test or documentation fixture");
        let sigma_star = 0.082_735_513_672_510_96;
        let boundary = RawSvi::new(0.1 * sigma_star, 0.5, 0.0, 0.0, sigma_star)
            .expect("valid test or documentation fixture");
        let pass_assessment = assess_raw(&pass, 1e-12);
        let fail_assessment = assess_raw(&fail, 1e-12);
        let boundary_assessment = assess_raw(&boundary, 1e-12);
        assert_eq!(
            pass_assessment.status(),
            ArbitrageStatus::NoViolationDetected
        );
        assert_eq!(fail_assessment.status(), ArbitrageStatus::ViolationDetected);
        assert_eq!(boundary_assessment.status(), ArbitrageStatus::Indeterminate);
        assert!(
            g(
                &fail,
                fail_assessment
                    .witness()
                    .expect("valid test or documentation fixture")
            ) < 0.0
        );
        assert!(
            (pass_assessment
                .witness()
                .expect("valid test or documentation fixture")
                .abs()
                - 0.356_280_635)
                .abs()
                < 5e-4
        );
        assert!(
            (fail_assessment
                .witness()
                .expect("valid test or documentation fixture")
                .abs()
                - 0.142_512_254)
                .abs()
                < 5e-4
        );
    }

    #[test]
    fn ssvi_slice_passing_theorem_42_is_butterfly_free() {
        // An SSVI slice satisfying Theorem 4.2 should also pass the g-scan.
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.butterfly_assessment_at(0.04).status(),
            ArbitrageStatus::NoViolationDetected
        );
        let raw = ssvi
            .slice_at(0.04)
            .expect("valid test or documentation fixture");
        let report = butterfly_scan(&raw, -1.0, 1.0).expect("valid test or documentation fixture");
        assert!(!report.violation_observed, "min_g = {}", report.min_g);
    }
}