Time doesn’t flow uniformly.
Neither should your pricing.

// Black-76 call — Heritage C++ (15 years production)
constexpr inline double call_price(
    double F, double K, double D,
    double sigma, double T) noexcept {
    const auto total_vol = sigma * std::sqrt(T);
    const auto d1 = (std::log(F / K) + 0.5 * sigma * sigma * T) / total_vol;
    const auto d2 = d1 - total_vol;
    return D * (F * N(d1) - K * N(d2));
}

// Black-76 call — C# at parity with Rust
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static OneOf<double, Black76Error> PriceOfCallNd(
    INormDist nd,
    double forward, double strike,
    double timeToExpiry, double volatility) {
    if (forward <= 0) return Black76Error.NegativeForward;
    var tVol = volatility * Math.Sqrt(timeToExpiry);
    var d1 = Math.Log(forward / strike) / tVol + 0.5 * tVol;
    var d2 = d1 - tVol;
    return Math.Max(forward * nd.Cdf(d1) - strike * nd.Cdf(d2), 0);
}

// nano-math.wit — package xll-ai:nano-math@6.4.12
interface special-math {
    phi:  func(x: f64) -> result<f64, math-error>;
    erfc: func(x: f64) -> result<f64, math-error>;
}

interface decimal-math {
    record decimal28 { lo: u32, mid: u32, hi: u32, flags: u32 }
    phi: func(x: f64) -> result<decimal28, math-error>;
}

// CUDA kernel — actual code from the repo
__host__ __device__ inline float norm_cdf(float x) {
    return 0.5f * erfcf(-x * 0.70710678118654752440f);
}

__global__ void black76_kernel(
    int n, const float* F, const float* K,
    const float* T, const float* V,
    const float* DF, float* C, float* P) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) black76_one(F[i], K[i], T[i], V[i], DF[i], C[i], P[i]);
}

# Generated from MathIR — same test vectors as Rust
def call_price(F: float, K: float, D: float,
               sigma: float, T: float) -> MathResult:
    if math.isnan(F): return (False, float('nan'))
    total_vol = sigma * math.sqrt(T)
    d1 = (math.log(F / K) + 0.5 * sigma**2 * T) / total_vol
    d2 = d1 - total_vol
    return (True, D * (F * norm_cdf(d1) - K * norm_cdf(d2)))

BLACK76.CALL = LAMBDA(F, K, D, sigma, T,
  LET(
    tv,  sigma * SQRT(T),
    d1,  (LN(F/K) + 0.5*sigma^2*T) / tv,
    d2,  d1 - tv,
    D * (F * NORM.S.DIST(d1,1) - K * NORM.S.DIST(d2,1))
  )
)

// Svensson zero rate — Heritage C++
constexpr inline double zero_rate(
    const SvenssonParams& p, double tau) noexcept {
    const auto e1 = std::exp(-tau / p.tau1);
    const auto e2 = std::exp(-tau / p.tau2);
    const auto t1 = (1.0 - e1) * p.tau1 / tau;
    const auto t2 = (1.0 - e2) * p.tau2 / tau;
    return p.beta0 + p.beta1*t1 + p.beta2*(t1-e1) + p.beta3*(t2-e2);
}

// Svensson zero rate — C# at parity with Rust
public static OneOf<double, CurveError> ZeroRate(
    SvenssonParams p, double tau) {
    if (tau <= 0) return CurveError.NonPositiveTenor;
    var e1 = Math.Exp(-tau / p.Tau1);
    var e2 = Math.Exp(-tau / p.Tau2);
    var t1 = (1.0 - e1) * p.Tau1 / tau;
    var t2 = (1.0 - e2) * p.Tau2 / tau;
    return p.Beta0 + p.Beta1*t1 + p.Beta2*(t1-e1) + p.Beta3*(t2-e2);
}

// Svensson — generated from MathIR
export function zeroRate(p, tau) {
    if (tau <= 0) return { ok: false, error: 'non_positive_tenor' };
    const e1 = Math.exp(-tau / p.tau1);
    const e2 = Math.exp(-tau / p.tau2);
    const t1 = (1 - e1) * p.tau1 / tau;
    const t2 = (1 - e2) * p.tau2 / tau;
    return { ok: true, value: p.beta0+p.beta1*t1+p.beta2*(t1-e1)+p.beta3*(t2-e2) };
}

# Svensson — generated from MathIR
def zero_rate(p: SvenssonParams, tau: float) -> Result:
    if tau <= 0: return (False, 0.0)
    e1 = math.exp(-tau / p.tau1)
    e2 = math.exp(-tau / p.tau2)
    t1 = (1 - e1) * p.tau1 / tau
    t2 = (1 - e2) * p.tau2 / tau
    return (True, p.beta0+p.beta1*t1+p.beta2*(t1-e1)+p.beta3*(t2-e2))

// SABR implied vol — C# at parity with Rust
public static OneOf<double, SabrError> ImpliedVol(
    INormDist nd, SabrParams p,
    double forward, double strike, double expiry) {
    var fk = Math.Pow(forward * strike, (1.0 - p.Beta) / 2.0);
    var z = (p.Nu / p.Alpha) * fk * Math.Log(forward / strike);
    var x = Math.Log((Math.Sqrt(1 - 2*p.Rho*z + z*z) + z - p.Rho)
            / (1 - p.Rho));
    // Hagan expansion + correction
    return p.Alpha * z / (fk * x) * correction;
}

// SABR — Heritage C++, production since 2012
inline double implied_vol(
    const SabrParams& p,
    double F, double K, double T) noexcept {
    const auto fk = std::pow(F * K, (1.0 - p.beta) / 2.0);
    const auto z = (p.nu / p.alpha) * fk * std::log(F / K);
    // ... Hagan expansion
    return p.alpha * z / (fk * x) * correction;
}

// reflexlibs-options.wit — SABR interface
interface sabr-model {
    record sabr-params {
        alpha: f64, beta: f64,
        rho: f64, nu: f64,
    }
    implied-vol: func(
        params: sabr-params,
        forward: f64, strike: f64, expiry: f64,
    ) -> result<f64, pricing-error>;
}

// Batch option pricing — one thread per option
__global__ void black76_kernel(
    int n, const float* F, const float* K,
    const float* T, const float* V,
    const float* DF, float* C, float* P) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) black76_one(F[i], K[i], T[i], V[i], DF[i], C[i], P[i]);
}

# SABR implied vol — Python research bridge
def implied_vol(p: SabrParams,
                F: float, K: float, T: float) -> Result:
    fk = (F * K) ** ((1 - p.beta) / 2)
    z = (p.nu / p.alpha) * fk * math.log(F / K)
    x = math.log((math.sqrt(1 - 2*p.rho*z + z**2) + z - p.rho) / (1 - p.rho))
    # Hagan expansion
    return (True, p.alpha * z / (fk * x) * correction)

/// Heston Normal — kernel signatures, full bodies in heston_normal.rs
pub fn price_of_call(
    forward: f64, strike: f64, tau: f64,
    kappa: f64, theta: f64, omega: f64,
    rho: f64, v0: f64,
) -> Result<f64, HestonError> { /* … Gauss-Legendre over φ(iε−ξ) */ }

pub fn delta_of_call(/* same signature */) -> Result<f64, HestonError> { /* … */ }
pub fn gamma        (/* same signature */) -> Result<f64, HestonError> { /* … */ }
pub fn theta_forward(/* same signature */) -> Result<f64, HestonError> { /* … */ }

/// Characteristic function φ — paper Φ(τ,u,ν₀) = exp(C + Dν₀ − iux₀)
fn phi(u: Complex<f64>, /* params … */) -> Complex<f64> { /* … */ }

// HestonNormal.h — 6 templated entry points, one file per pricing family
template<class T>
T price_of_call(T F, T K, T tau,
                  T kappa, T theta, T omega, T rho, T v0) noexcept;

template<class T> T price_of_put (/* same signature */) noexcept;
template<class T> T delta_of_call(/* … */)                    noexcept;
template<class T> T delta_of_put (/* … */)                    noexcept;
template<class T> T gamma        (/* … */)                    noexcept;
template<class T> T theta_forward(/* … */)                    noexcept;

// Implementation: bodies in HestonNormal.cpp — Gauss-Legendre quadrature,
// characteristic function in HestonNormalCharacteristic.cpp

// HestonNormalFromPaper.cs — Path A (clean-room from SSRN 2445328)
public static class HestonNormalFromPaper {
    public static double PriceOfCall(
        double forward, double strike, double timeToExpiry,
        double kappa, double theta, double omega,
        double rho, double v0,
        int nodeCount = 1024) { /* … integrate Re[φ·(−e^…K)/(ξ−iε)²] / π */ }

    public static double PriceOfPut  (/* same signature */) { /* … put-call parity */ }
    public static double DeltaOfCall (/* … */) { /* … multiplier −iu */ }
    public static double Gamma       (/* … */) { /* … multiplier −u² */ }
    public static double ThetaForward(/* … */) { /* … ∂φ/∂τ via ODE 9a/9b */ }

    // Φ(τ,u,ν₀) = exp(C(τ) + D(τ)ν₀ − iu·x₀)  — paper §3, negative-i convention
    private static (Complex phi, Complex phiDot) CharacteristicFunction(/* … */) { /* … */ }
}

// HestonNormalFromHeritage.cs — Path B (under permission grant)
public static class HestonNormalFromHeritage { /* same 6 entry points, +i convention */ }

/// Bates Lognormal — kernel signatures (Heston params + jump params)
pub fn price_of_call(
    forward: f64, strike: f64, tau: f64,
    kappa: f64, theta: f64, omega: f64, rho: f64, v0: f64,
    lambda_j: f64, mu_j: f64, sigma_j: f64,
) -> Result<f64, BatesError> { /* … Heston φ · jump-CF … */ }

pub fn delta_of_call(/* same signature */) -> Result<f64, BatesError> { /* … */ }
pub fn gamma        (/* same signature */) -> Result<f64, BatesError> { /* … */ }
pub fn theta_forward(/* same signature */) -> Result<f64, BatesError> { /* … */ }

/// φ_Bates(u) = φ_Heston(u) · exp(τ·λ(e^{iμ_j u − ½σ_j² u²} − 1 − iμu))
fn jump_characteristic(u: Complex<f64>, /* … */) -> Complex<f64> { /* … */ }

// BatesLognormal.h — 6 entry points per variant (Normal/Lognormal/LognormalShifted)
template<class T>
T price_of_call(T F, T K, T tau,
                  T kappa, T theta, T omega, T rho, T v0,
                  T lambda_j, T mu_j, T sigma_j) noexcept;

template<class T> T price_of_put (/* same signature */) noexcept;
template<class T> T delta_of_call(/* … */)                    noexcept;
template<class T> T delta_of_put (/* … */)                    noexcept;
template<class T> T gamma        (/* … */)                    noexcept;
template<class T> T theta_forward(/* … */)                    noexcept;

// Variants: BatesLognormal, BatesNormal, BatesLognormalShifted
// Shared CF plumbing with HestonLognormal, +Merton jump component

// BatesLognormalFromPaper.cs — Path A (queued, Bates 1996 RFS)
public static class BatesLognormalFromPaper {
    public static double PriceOfCall(
        double forward, double strike, double timeToExpiry,
        double kappa, double theta, double omega, double rho, double v0,
        double lambdaJ, double muJ, double sigmaJ,
        int nodeCount = 1024) { /* … Heston φ × jump-CF … */ }

    public static double PriceOfPut  (/* same signature */) { /* … */ }
    public static double DeltaOfCall (/* … */)                         { /* … */ }
    public static double Gamma       (/* … */)                         { /* … */ }
    public static double ThetaForward(/* … */)                         { /* … */ }
}

// BatesLognormalFromHeritage.cs — Path B under permission grant
public static class BatesLognormalFromHeritage { /* same 6 entry points */ }

// IFixedRateSwap — C++17 Delta One
template <class C>
class IFixedRateSwap : public ISwap {
public:
    virtual const IFixedLeg<C>& fixed_leg() const = 0;
    virtual const IFloatLeg<C>& float_leg() const = 0;
    virtual const C& paying_currency() const = 0;
};

// AccrualPeriod — no chrono dependency at contract level
struct AccrualPeriod {
    Ymd start;
    Ymd end;
    std::optional<StubData> stub;
};

// Discriminants match Rust #[repr(u8)] and C# : byte exactly
enum class Basis : std::uint8_t {
    ActualActualIsda = 4,
    Actual360Fixed   = 5,
    Actual365Fixed   = 6,
    // ... 14 conventions total
};

// IFixedRateSwap — Rustified C# from Heritage
public interface IFixedRateSwap : ISwap {
    IFixedLeg FixedLeg { get; }
    IFloatLeg FloatLeg { get; }
    ICurrency PayingCurrency { get; }
}

// AccrualPeriod — day count, compounding, spread
public readonly record struct AccrualPeriod(
    DateOnly Start, DateOnly End,
    DayCountBasis Basis, double Notional);

// ISwaption — Heritage C#
public interface ISwaption : IOption {
    IFixedRateSwap UnderlyingSwap { get; }
    SettlementType Settlement { get; }
    CallOrPut CallOrPut { get; }
    DateOnly Expiry { get; }
}

public enum SettlementType { Cash, Physical }
public enum CallOrPut { Call, Put }

// Curve sensitivity — C# bump-and-reprice
public static double[] ZcSensitivities(
    IDiscountCurve baseCurve,
    ReadOnlySpan<double> tenors,
    double bumpSize = 0.0001) {
    var result = new double[tenors.Length];
    for (var i = 0; i < tenors.Length; i++) {
        var bumped = baseCurve.Bump(tenors[i], bumpSize);
        result[i] = Reprice(bumped) - Reprice(baseCurve);
    }
    return result;
}

// EMA — C# with Span<T> streaming
public sealed class Ema {
    private readonly double _alpha;
    private double _value;
    private bool _init;

    public Ema(int period) => _alpha = 2.0 / (period + 1);

    public double Update(double price) {
        if (!_init) { _value = price; _init = true; }
        else _value = _alpha * price + (1 - _alpha) * _value;
        return _value;
    }
}