solmath 0.2.0

// European barrier option pricing via Rubinstein-Reiner building blocks.
//
// Uses Haug building blocks A, B, C, D with eta = phi (not barrier direction).
// Verified against QuantLib AnalyticBarrierEngine on 443K vectors.
//
// All arithmetic at HP precision (1e15). Final SBF audit: 270,156 CU average,
// 415,531 max for the legacy/unbreached calculation.

use crate::arithmetic::{fp_div_i, fp_mul_i, isqrt_u128};
use crate::constants::*;
use crate::error::SolMathError;
use crate::hp::{
    black_scholes_price_hp, downscale_hp_to_std, exp_fixed_hp, fp_div_hp_safe, fp_mul_hp_i,
    ln_fixed_hp, norm_cdf_poly_hp, upscale_std_to_hp,
};

/// Barrier option type (single barrier, European exercise).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BarrierType {
    /// Knocked out if spot falls to or below the barrier.
    DownAndOut,
    /// Knocked in (activated) if spot falls to or below the barrier.
    DownAndIn,
    /// Knocked out if spot rises to or above the barrier.
    UpAndOut,
    /// Knocked in (activated) if spot rises to or above the barrier.
    UpAndIn,
}

/// Result of a barrier option pricing computation.
///
/// The returned `price` and `vanilla` are rounded so that paired knock-in/knock-out
/// calls satisfy the exact public identity `in_price + out_price == vanilla`.
#[derive(Debug, Clone, Copy)]
pub struct BarrierResult {
    /// Barrier option price at SCALE.
    pub price: u128,
    /// Vanilla BS price for reference (in + out = vanilla).
    pub vanilla: u128,
}

/// All HP intermediates for barrier pricing.
struct HaugIntermediates {
    s_hp: i128,
    k_disc_hp: i128,
    x1_hp: i128,
    y1_hp: i128,
    d1_hp: i128,
    y_hp: i128,
    sigma_sqrt_t_hp: i128,
    discount_hp: i128,
    pow_2l_hp: i128,
    pow_2lm2_hp: i128,
    phi: i128,
}

/// Compute all HP intermediates for barrier pricing.
/// eta = phi (call/put sign), NOT barrier direction.
#[inline(never)]
fn compute_intermediates(
    s: u128,
    k: u128,
    h: u128,
    r: u128,
    sigma: u128,
    t: u128,
    is_call: bool,
) -> Result<HaugIntermediates, SolMathError> {
    let s_hp = upscale_std_to_hp(s)?;
    let k_hp = upscale_std_to_hp(k)?;
    let h_hp = upscale_std_to_hp(h)?;
    let r_hp = upscale_std_to_hp(r)?;
    let sigma_hp = upscale_std_to_hp(sigma)?;
    let t_hp = upscale_std_to_hp(t)?;

    let sqrt_t_hp = isqrt_u128(
        (t_hp as u128)
            .checked_mul(SCALE_HP_U)
            .ok_or(SolMathError::Overflow)?,
    ) as i128;
    let sigma_sqrt_t_hp = fp_mul_hp_i(sigma_hp, sqrt_t_hp)?;

    let r_t_hp = fp_mul_hp_i(r_hp, t_hp)?;
    let discount_hp = exp_fixed_hp(-r_t_hp)?;
    let k_disc_hp = fp_mul_hp_i(k_hp, discount_hp)?;

    let sigma_sq_hp = fp_mul_hp_i(sigma_hp, sigma_hp)?;
    let drift_rate_hp = r_hp
        .checked_add(sigma_sq_hp / 2)
        .ok_or(SolMathError::Overflow)?;
    let drift_hp = fp_mul_hp_i(drift_rate_hp, t_hp)?;
    let lambda_sst = if sigma_sqrt_t_hp > 0 {
        fp_div_hp_safe(drift_hp, sigma_sqrt_t_hp)?
    } else {
        0
    };

    let ln_sk = ln_fixed_hp(fp_div_hp_safe(s_hp, k_hp)?)?;
    let ln_sh = ln_fixed_hp(fp_div_hp_safe(s_hp, h_hp)?)?;
    let ln_hk = ln_fixed_hp(fp_div_hp_safe(h_hp, k_hp)?)?;

    let mk = |log_val: i128| -> Result<i128, SolMathError> {
        if sigma_sqrt_t_hp > 0 {
            // fp_div_hp_safe result ∈ [-~1e15, ~1e15]; lambda_sst ∈ [-~1e15, ~1e15] (finite-rate drift); sum ≤ ~2e15, fits i128
            fp_div_hp_safe(log_val, sigma_sqrt_t_hp)?
                .checked_add(lambda_sst)
                .ok_or(SolMathError::Overflow)
        } else {
            Ok(0)
        }
    };

    let d1_hp = mk(ln_sk)?;
    let x1_hp = mk(ln_sh)?;
    let y1_hp = mk(-ln_sh)?;
    // -ln_sh ∈ [-~1e15, ~1e15], ln_hk ∈ [-~1e15, ~1e15]; sum ≤ ~2e15, fits i128
    let y_hp = mk(ln_sh
        .checked_neg()
        .and_then(|v| v.checked_add(ln_hk))
        .ok_or(SolMathError::Overflow)?)?;

    // Power terms at HP via exp(2λ·ln(H/S))
    let sigma_sq_std = fp_mul_i(sigma as i128, sigma as i128)?;
    // r as i128 ≤ ~1e12 (rate at SCALE), sigma_sq_std ≤ SCALE (volatility² ≤ 1.0 at SCALE); sum ≤ ~2e12, fits i128
    let lambda_num = (r as i128)
        .checked_add(sigma_sq_std / 2)
        .and_then(|v| v.checked_mul(2))
        .ok_or(SolMathError::Overflow)?;
    let two_lambda_std = fp_div_i(lambda_num, sigma_sq_std)?;
    let two_lambda_hp = upscale_std_to_hp(two_lambda_std as u128)?;
    // two_lambda_hp ≤ ~100·SCALE_HP (lambda is a dimensionless financial ratio, typically ≤ 100); 2·SCALE_HP ≈ 2e15; no underflow for lambda > 1
    let two_lambda_m2_hp = two_lambda_hp - 2 * SCALE_HP;
    let ln_h_over_s_hp = -ln_sh;

    let pow_2l_hp = if ln_h_over_s_hp == 0 {
        SCALE_HP
    } else {
        exp_fixed_hp(fp_mul_hp_i(two_lambda_hp, ln_h_over_s_hp)?)?
    };
    let pow_2lm2_hp = if ln_h_over_s_hp == 0 {
        SCALE_HP
    } else {
        exp_fixed_hp(fp_mul_hp_i(two_lambda_m2_hp, ln_h_over_s_hp)?)?
    };

    Ok(HaugIntermediates {
        s_hp,
        k_disc_hp,
        x1_hp,
        y1_hp,
        d1_hp,
        y_hp,
        sigma_sqrt_t_hp,
        discount_hp,
        pow_2l_hp,
        pow_2lm2_hp,
        phi: if is_call { 1 } else { -1 },
    })
}

/// Compute a single Haug building block at HP.
/// block(z) = φ·[s_eff·N(eta·z) - k_eff·N(eta·(z - σ√T))]
/// where eta = phi.
#[inline(never)]
fn block_hp(phi: i128, z: i128, s_eff: i128, k_eff: i128, sst: i128) -> Result<i128, SolMathError> {
    // eta = phi for all blocks
    // phi ∈ {-1, +1} (literal scalar, not SCALE-valued); z and sst are HP-scale ∈ [-~1e15, ~1e15]
    // z - sst: both ≤ ~1e15; difference ≤ ~2e15, fits i128
    // phi * z, phi * (z - sst): sign flips only, magnitude unchanged, fits i128
    // fp_mul_hp_i outputs ∈ [-~1e15, ~1e15] (price × N(·) where N ∈ [0,1]); difference ≤ ~2e15, fits i128
    // outer phi * (...): sign flip, magnitude unchanged; fits i128
    Ok(phi
        * (fp_mul_hp_i(s_eff, norm_cdf_poly_hp(phi * z)?)?
            - fp_mul_hp_i(k_eff, norm_cdf_poly_hp(phi * (z - sst))?)?))
}

/// Compute all 4 building blocks: (A, B, C, D).
#[inline(never)]
fn all_blocks(im: &HaugIntermediates) -> Result<(i128, i128, i128, i128), SolMathError> {
    let s_pow = fp_mul_hp_i(im.s_hp, im.pow_2l_hp)?;
    let k_pow = fp_mul_hp_i(im.k_disc_hp, im.pow_2lm2_hp)?;

    let a = block_hp(im.phi, im.x1_hp, im.s_hp, im.k_disc_hp, im.sigma_sqrt_t_hp)?;
    let b = block_hp(im.phi, im.d1_hp, im.s_hp, im.k_disc_hp, im.sigma_sqrt_t_hp)?;
    let c = block_hp(im.phi, im.y1_hp, s_pow, k_pow, im.sigma_sqrt_t_hp)?;
    let d = block_hp(im.phi, im.y_hp, s_pow, k_pow, im.sigma_sqrt_t_hp)?;

    Ok((a, b, c, d))
}

/// Single barrier European option price via Rubinstein-Reiner building blocks.
///
/// Prices a European option with a single knock-in or knock-out barrier
/// using Haug's ABCD decomposition, verified against QuantLib on 443K vectors.
///
/// This formula assumes continuous monitoring, zero rebate, no dividends, and
/// that the barrier has **not** been breached before the valuation instant.
/// On-chain callers with persisted path state should use
/// [`barrier_option_with_state`]. Discretely sampled oracle barriers require a
/// separate monitoring correction and must not be priced as continuous.
///
/// # Parameters
/// - `s` -- Spot price at SCALE (u128)
/// - `k` -- Strike price at SCALE (u128)
/// - `h` -- Barrier level at SCALE (u128)
/// - `r` -- Risk-free rate at SCALE (u128, e.g. `50_000_000_000` = 5%)
/// - `sigma` -- Volatility at SCALE (u128, e.g. `250_000_000_000` = 25%)
/// - `t` -- Time to expiry in years at SCALE (u128)
/// - `is_call` -- `true` for call, `false` for put
/// - `barrier_type` -- [`BarrierType`] variant
///
/// # Errors
/// Returns `Err(DomainError)` if `s`, `k`, `h`, `sigma`, or `t` are zero.
///
/// # Accuracy
/// Max 1.7K ULP, P99 33, median 1. Final SBF audit: 270,156 CU
/// average and 415,531 max for this math call.
///
/// Public return values preserve exact in/out conservation after rounding.
///
/// # Example
/// ```
/// use solmath::{barrier_option, BarrierType, SCALE};
/// let result = barrier_option(
///     100 * SCALE, 100 * SCALE, 90 * SCALE,
///     50_000_000_000, 250_000_000_000, SCALE,
///     true, BarrierType::DownAndOut,
/// )?;
/// assert!(result.price > 0);
/// assert!(result.price <= result.vanilla);
/// # Ok::<(), solmath::SolMathError>(())
/// ```
pub fn barrier_option(
    s: u128,
    k: u128,
    h: u128,
    r: u128,
    sigma: u128,
    t: u128,
    is_call: bool,
    barrier_type: BarrierType,
) -> Result<BarrierResult, SolMathError> {
    if s == 0 || k == 0 || sigma == 0 || t == 0 || h == 0 {
        return Err(SolMathError::DomainError);
    }

    let is_down = matches!(
        barrier_type,
        BarrierType::DownAndOut | BarrierType::DownAndIn
    );
    let is_out = matches!(
        barrier_type,
        BarrierType::DownAndOut | BarrierType::UpAndOut
    );

    // Already at or past barrier
    if (is_down && s <= h) || (!is_down && s >= h) {
        let (call, put) = black_scholes_price_hp(s, k, r, sigma, t)?;
        let vanilla = if is_call { call } else { put };
        return Ok(BarrierResult {
            price: if is_out { 0 } else { vanilla },
            vanilla,
        });
    }

    // Impossible payoff: up call K≥H, down put K≤H
    if is_call && !is_down && k >= h {
        let (call, _) = black_scholes_price_hp(s, k, r, sigma, t)?;
        return Ok(BarrierResult {
            price: if is_out { 0 } else { call },
            vanilla: call,
        });
    }
    if !is_call && is_down && k <= h {
        let (_, put) = black_scholes_price_hp(s, k, r, sigma, t)?;
        return Ok(BarrierResult {
            price: if is_out { 0 } else { put },
            vanilla: put,
        });
    }

    let im = compute_intermediates(s, k, h, r, sigma, t, is_call)?;
    let (a, b, c, d) = all_blocks(&im)?;
    let vanilla_hp = b;

    // Select formula based on verified QuantLib match:
    //   Down call K≥H: out = B-D        Down call K<H: out = A-C
    //   Down put  K>H: out = B-A+C-D    Up call   K<H: out = B-A+C-D
    //   Up put   K≤H: out = B-D         Up put    K>H: digital decomposition
    let out_hp = if !is_down && !is_call && k > h {
        // Up put K > H: digital decomposition
        let im_h = compute_intermediates(s, h, h, r, sigma, t, false)?;
        let (_, b_h, _, d_h) = all_blocks(&im_h)?;
        // b_h, d_h are Haug blocks at HP ∈ [-~1e20, ~1e20] (price × N(·)); d_h ≤ b_h by construction; no underflow
        let p_uo_h_hp = b_h - d_h;

        // Digital: disc · [N(σ√T - x₁) - (H/S)^α · N(σ√T - y₁)]
        // sigma_sqrt_t_hp - x1_hp: both ∈ [-~1e15, ~1e15]; difference ≤ ~2e15, fits i128
        // sigma_sqrt_t_hp - y1_hp: same reasoning
        // N(·) ∈ [0, SCALE_HP]; fp_mul_hp_i result ∈ [0, SCALE_HP]; difference ≤ SCALE_HP, fits i128
        let digital_hp = fp_mul_hp_i(
            im.discount_hp,
            norm_cdf_poly_hp(im.sigma_sqrt_t_hp - im.x1_hp)?
                - fp_mul_hp_i(
                    im.pow_2lm2_hp,
                    norm_cdf_poly_hp(im.sigma_sqrt_t_hp - im.y1_hp)?,
                )?,
        )?;

        // p_uo_h_hp ∈ [-~1e20, ~1e20]; fp_mul_hp_i of (k-h) upscaled × digital ∈ [-~1e20, ~1e20]; sum ≤ ~2e20, fits i128
        p_uo_h_hp + fp_mul_hp_i(upscale_std_to_hp(k - h)?, digital_hp)?
    } else if is_down && is_call && k < h {
        // Down call K < H: out = A - C
        // a, c: Haug blocks at HP ∈ [-~1e20, ~1e20]; a ≥ c by formula construction; difference ∈ [-~1e20, ~1e20], fits i128
        a - c
    } else if (is_down && !is_call && k > h) || (!is_down && is_call) {
        // Straddling: down put K>H or up call K<H: out = B - A + C - D
        // a, b, c, d all Haug blocks at HP ∈ [-~1e20, ~1e20]; cumulative sum of four terms ≤ ~4e20, fits i128
        b - a + c - d
    } else {
        // Non-straddling: down call K≥H or up put K≤H: out = B - D
        // b, d: Haug blocks at HP ∈ [-~1e20, ~1e20]; d ≤ b by formula construction; no underflow; fits i128
        b - d
    };

    let vanilla = downscale_hp_to_std(vanilla_hp);
    let out_price = core::cmp::min(downscale_hp_to_std(out_hp), vanilla);
    let price = if is_out {
        out_price
    } else {
        vanilla - out_price
    };

    Ok(BarrierResult { price, vanilla })
}

/// Path-state-aware barrier pricing.
///
/// Set `barrier_was_breached` from persisted contract/oracle state. Once
/// breached, a knock-out is worth zero and a knock-in is worth the vanilla
/// option regardless of the current spot.
pub fn barrier_option_with_state(
    s: u128,
    k: u128,
    h: u128,
    r: u128,
    sigma: u128,
    t: u128,
    is_call: bool,
    barrier_type: BarrierType,
    barrier_was_breached: bool,
) -> Result<BarrierResult, SolMathError> {
    if !barrier_was_breached {
        return barrier_option(s, k, h, r, sigma, t, is_call, barrier_type);
    }
    if s == 0 || k == 0 || h == 0 || sigma == 0 || t == 0 {
        return Err(SolMathError::DomainError);
    }
    let (call, put) = black_scholes_price_hp(s, k, r, sigma, t)?;
    let vanilla = if is_call { call } else { put };
    let knocked_out = matches!(
        barrier_type,
        BarrierType::DownAndOut | BarrierType::UpAndOut
    );
    Ok(BarrierResult {
        price: if knocked_out { 0 } else { vanilla },
        vanilla,
    })
}

#[cfg(test)]
mod path_state_tests {
    use super::*;

    #[test]
    fn historical_breach_overrides_current_safe_spot() {
        let out = barrier_option_with_state(
            100 * SCALE,
            100 * SCALE,
            90 * SCALE,
            50_000_000_000,
            200_000_000_000,
            SCALE,
            true,
            BarrierType::DownAndOut,
            true,
        )
        .unwrap();
        let knocked_in = barrier_option_with_state(
            100 * SCALE,
            100 * SCALE,
            90 * SCALE,
            50_000_000_000,
            200_000_000_000,
            SCALE,
            true,
            BarrierType::DownAndIn,
            true,
        )
        .unwrap();
        assert_eq!(out.price, 0);
        assert_eq!(knocked_in.price, knocked_in.vanilla);
    }
}