libitofin 0.52.1

//! Stateless, seeded Monte Carlo building blocks shared across bindings.
//!
//! Draws use the core Mersenne Twister / inverse-normal policy, consumed in
//! path, time, asset order. Nonzero seeds reproduce the native sequence; zero
//! is rejected because the native RNG interprets it as a nondeterministic seed.

use crate::errors::{QlError, QlResult};
use crate::handle::Handle;
use crate::interestrate::Compounding;
use crate::math::array::Array;
use crate::math::matrix::Matrix;
use crate::math::matrixutilities::choleskydecomposition::cholesky_decomposition;
use crate::math::randomnumbers::rngtraits::{McRngTraits, PseudoRandom, SequenceGenerator};
use crate::processes::{HestonDiscretization as Discretization, HestonProcess};
use crate::quotes::make_quote_handle;
use crate::require;
use crate::shared::{Shared, shared};
use crate::stochasticprocess::StochasticProcess;
use crate::termstructures::yields::FlatForward;
use crate::termstructures::yieldtermstructure::YieldTermStructure;
use crate::time::date::{Date, Month};
use crate::time::daycounters::actual365fixed::Actual365Fixed;
use crate::time::frequency::Frequency;
use crate::types::{Natural, Rate, Real, Size, Time, Volatility};

fn error(message: &str) -> QlError {
    QlError::new(message, file!(), line!())
}

fn zeros(count: Size) -> QlResult<Vec<Real>> {
    let mut result = Vec::new();
    result
        .try_reserve_exact(count)
        .map_err(|_| error("allocation failed"))?;
    result.resize(count, 0.0);
    Ok(result)
}

fn product(a: Size, b: Size) -> QlResult<Size> {
    a.checked_mul(b)
        .ok_or_else(|| error("simulation dimensions overflow"))
}

/// Independent standard-normal draws in native sequence order.
///
/// An empty request returns an empty vector, but still requires a nonzero seed.
pub fn gaussian_draws(count: Size, seed: Natural) -> QlResult<Vec<Real>> {
    require!(seed != 0, "seed must be nonzero");
    let mut result = zeros(count)?;
    // A scalar sequence keeps native internal allocations constant-sized.
    let mut rng = PseudoRandom::make_sequence_generator(1, seed)?;
    for value in &mut result {
        *value = rng.next_sequence().value[0];
    }
    Ok(result)
}

/// Parameters for exact-discretization geometric Brownian motion.
///
/// `drift` is the annualized arithmetic drift (risk-neutral callers supply
/// their own rate minus yield); volatility and horizon use matching time units.
/// Correlation is a dense row-major, symmetric, unit-diagonal **positive
/// definite** matrix. Singular positive-semidefinite matrices are rejected;
/// no identity fallback, jitter, or eigenvalue repair is performed.
pub struct GbmRequest<'a> {
    pub initial: &'a [Real],
    pub drift: &'a [Rate],
    pub volatility: &'a [Volatility],
    pub correlation: &'a [Real],
    pub horizon: Time,
    pub steps: Size,
    pub paths: Size,
    pub seed: Natural,
    /// Return only `[path, asset]`, with the same terminal values as full mode.
    pub terminal_only: bool,
}

/// Number of result elements, with overflow checked before allocation.
pub fn output_len(request: &GbmRequest<'_>) -> QlResult<Size> {
    let times = if request.terminal_only {
        1
    } else {
        request
            .steps
            .checked_add(1)
            .ok_or_else(|| error("step count overflows"))?
    };
    product(product(request.paths, times)?, request.initial.len())
}

fn correlation_factor(values: &[Real], assets: Size) -> QlResult<Matrix> {
    require!(
        values.len() == product(assets, assets)?,
        "correlation shape mismatch"
    );
    // Strict input checks prevent the native decomposition's flexible mode
    // from silently interpreting an asymmetric matrix using its upper half.
    for i in 0..assets {
        require!(
            values[i * assets + i] == 1.0,
            "correlation diagonal must equal one"
        );
        for j in 0..assets {
            let value = values[i * assets + j];
            require!(
                value.is_finite() && value.abs() <= 1.0,
                "invalid correlation entry"
            );
            require!(
                value == values[j * assets + i],
                "correlation must be symmetric"
            );
        }
    }
    let mut matrix = Matrix::with_size(assets, assets);
    for i in 0..assets {
        matrix
            .row_mut(i)
            .copy_from_slice(&values[i * assets..(i + 1) * assets]);
    }
    // Flexible mode avoids the core strict-mode panic. Require every pivot
    // positive and verify reconstruction so this remains an SPD-only API.
    let factor = cholesky_decomposition(&matrix, true);
    let tolerance = 64.0 * Real::EPSILON * assets as Real;
    for i in 0..assets {
        require!(
            matches!(
                factor[(i, i)].partial_cmp(&0.0),
                Some(std::cmp::Ordering::Greater)
            ),
            "correlation must be positive definite"
        );
        for j in 0..=i {
            let rebuilt: Real = (0..=j).map(|k| factor[(i, k)] * factor[(j, k)]).sum();
            require!(
                rebuilt.is_finite() && (rebuilt - matrix[(i, j)]).abs() <= tolerance,
                "correlation factorization failed"
            );
        }
    }
    Ok(factor)
}

/// Generate GBM paths, laid out `[path, time, asset]`, including time zero.
///
/// Uses `dt = horizon / steps` and `S *= exp((mu - sigma²/2)dt +
/// sigma*sqrt(dt)*(L*z))`. All paths consume the same number of draws,
/// including zero-volatility assets. Zero horizon returns initial values
/// exactly. Zero volatility follows deterministic `initial*exp(mu*t)`.
///
/// # Errors
/// Rejects nonfinite parameters, nonpositive spots, negative volatilities or
/// horizon, zero steps/paths/seed, inconsistent shapes, non-SPD correlation,
/// and nonfinite or underflowed prices. Supports at most 1,024 assets to bound
/// the core Cholesky routine's infallibly allocated matrix workspace. Output
/// and simulation-owned scratch allocations use `try_reserve_exact`.
pub fn gbm_paths(request: &GbmRequest<'_>) -> QlResult<Vec<Real>> {
    let assets = request.initial.len();
    require!(
        (1..=1024).contains(&assets),
        "asset count must be between 1 and 1024"
    );
    require!(
        request.steps > 0 && request.paths > 0,
        "steps and paths must be positive"
    );
    require!(request.seed != 0, "seed must be nonzero");
    require!(
        request.horizon.is_finite() && request.horizon >= 0.0,
        "invalid horizon"
    );
    require!(
        request.drift.len() == assets && request.volatility.len() == assets,
        "asset shape mismatch"
    );
    for i in 0..assets {
        require!(
            request.initial[i].is_finite() && request.initial[i] > 0.0,
            "invalid initial value"
        );
        require!(request.drift[i].is_finite(), "invalid drift");
        require!(
            request.volatility[i].is_finite() && request.volatility[i] >= 0.0,
            "invalid volatility"
        );
    }
    // Validate dimension products even when terminal mode does not store time.
    product(product(request.paths, request.steps)?, assets)?;
    let count = output_len(request)?;
    let factor = correlation_factor(request.correlation, assets)?;
    let mut output = zeros(count)?;
    let mut current = zeros(assets)?;
    let mut normals = zeros(assets)?;
    let mut rng = PseudoRandom::make_sequence_generator(1, request.seed)?;
    let dt = request.horizon / request.steps as Real;
    let sqrt_dt = dt.sqrt();
    let mut cursor = 0;
    for _ in 0..request.paths {
        current.copy_from_slice(request.initial);
        if !request.terminal_only {
            output[cursor..cursor + assets].copy_from_slice(&current);
            cursor += assets;
        }
        for step in 0..request.steps {
            for z in &mut normals {
                *z = rng.next_sequence().value[0];
            }
            for i in 0..assets {
                let sigma = request.volatility[i];
                let mu = request.drift[i];
                if request.horizon == 0.0 {
                    current[i] = request.initial[i];
                } else if sigma == 0.0 {
                    let time = if step + 1 == request.steps {
                        request.horizon
                    } else {
                        (step + 1) as Real * dt
                    };
                    current[i] = request.initial[i] * (mu * time).exp();
                } else {
                    let z: Real = (0..=i).map(|j| factor[(i, j)] * normals[j]).sum();
                    current[i] *= ((mu - 0.5 * sigma * sigma) * dt + sigma * sqrt_dt * z).exp();
                }
                require!(
                    current[i].is_finite() && current[i] > 0.0,
                    "simulated price overflow or underflow"
                );
            }
            if !request.terminal_only {
                output[cursor..cursor + assets].copy_from_slice(&current);
                cursor += assets;
            }
        }
        if request.terminal_only {
            output[cursor..cursor + assets].copy_from_slice(&current);
            cursor += assets;
        }
    }
    Ok(output)
}

/// Parameters for stateless Heston paths under flat risk-free and dividend rates.
///
/// The two independent normals at each step are consumed in path, time,
/// spot-factor, variance-factor order. Correlation is applied by the process.
pub struct HestonRequest {
    pub spot: Real,
    pub variance: Real,
    pub risk_free_rate: Rate,
    pub dividend_yield: Rate,
    pub kappa: Real,
    pub theta: Real,
    pub sigma: Real,
    pub rho: Real,
    pub horizon: Time,
    pub steps: Size,
    pub paths: Size,
    pub seed: Natural,
    pub discretization: Discretization,
    pub terminal_only: bool,
}

/// Number of spot/variance result elements, checked before allocation.
pub fn heston_output_len(request: &HestonRequest) -> QlResult<Size> {
    let times = if request.terminal_only {
        1
    } else {
        request
            .steps
            .checked_add(1)
            .ok_or_else(|| error("step count overflows"))?
    };
    product(product(request.paths, times)?, 2)
}

fn flat_curve(rate: Rate) -> Handle<dyn YieldTermStructure> {
    Handle::new(shared(FlatForward::with_rate(
        Date::new(1, Month::January, 2026),
        rate,
        Actual365Fixed::new(),
        Compounding::Continuous,
        Frequency::Annual,
    )) as Shared<dyn YieldTermStructure>)
}

/// Generate Heston spot/variance paths in `[path, time, component]` order.
///
/// Component zero is spot and component one is variance. Full output includes
/// time zero; terminal output is `[path, component]` and equals the last full
/// time row for the same seed. Zero horizon repeats the initial state exactly.
///
/// # Errors
/// Rejects invalid dimensions, zero seed, nonfinite values, spot <= 0,
/// variance < 0, kappa/theta/sigma <= 0, rho outside [-1,1], negative
/// horizon, unsupported schemes, allocation failures, and nonfinite states.
pub fn heston_paths(request: &HestonRequest) -> QlResult<Vec<Real>> {
    require!(
        request.steps > 0 && request.paths > 0,
        "steps and paths must be positive"
    );
    require!(request.seed != 0, "seed must be nonzero");
    require!(
        request.spot.is_finite() && request.spot > 0.0,
        "invalid initial spot"
    );
    require!(
        request.variance.is_finite() && request.variance >= 0.0,
        "invalid initial variance"
    );
    require!(
        request.risk_free_rate.is_finite() && request.dividend_yield.is_finite(),
        "invalid rate"
    );
    require!(
        request.kappa.is_finite() && request.kappa > 0.0,
        "invalid kappa"
    );
    require!(
        request.theta.is_finite() && request.theta > 0.0,
        "invalid theta"
    );
    require!(
        request.sigma.is_finite() && request.sigma > 0.0,
        "invalid sigma"
    );
    require!(
        request.rho.is_finite() && request.rho.abs() <= 1.0,
        "invalid rho"
    );
    require!(
        request.horizon.is_finite() && request.horizon >= 0.0,
        "invalid horizon"
    );
    require!(
        matches!(
            request.discretization,
            Discretization::QuadraticExponential | Discretization::QuadraticExponentialMartingale
        ),
        "unsupported Heston discretization"
    );
    product(product(request.paths, request.steps)?, 2)?;
    let count = heston_output_len(request)?;
    let dt = request.horizon / request.steps as Real;
    require!(request.horizon == 0.0 || dt > 0.0, "time step underflow");
    let process = HestonProcess::with_discretization(
        flat_curve(request.risk_free_rate),
        flat_curve(request.dividend_yield),
        make_quote_handle(request.spot).handle(),
        request.variance,
        request.kappa,
        request.theta,
        request.sigma,
        request.rho,
        request.discretization,
    );
    let mut output = zeros(count)?;
    let mut rng = PseudoRandom::make_sequence_generator(1, request.seed)?;
    let initial = [request.spot, request.variance];
    let mut cursor = 0;
    for _ in 0..request.paths {
        let mut state = Array::from(initial);
        if !request.terminal_only {
            output[cursor..cursor + 2].copy_from_slice(&initial);
            cursor += 2;
        }
        for step in 0..request.steps {
            let z_spot = rng.next_sequence().value[0];
            let z_variance = rng.next_sequence().value[0];
            if request.horizon > 0.0 {
                state = process.evolve(
                    step as Real * dt,
                    &state,
                    dt,
                    &Array::from([z_spot, z_variance]),
                )?;
                require!(
                    state[0].is_finite() && state[0] > 0.0,
                    "simulated spot overflow or underflow"
                );
                require!(
                    state[1].is_finite() && state[1] >= 0.0,
                    "invalid simulated variance"
                );
            }
            if !request.terminal_only {
                output[cursor..cursor + 2].copy_from_slice(&[state[0], state[1]]);
                cursor += 2;
            }
        }
        if request.terminal_only {
            output[cursor..cursor + 2].copy_from_slice(&[state[0], state[1]]);
            cursor += 2;
        }
    }
    Ok(output)
}

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

    fn heston_request() -> HestonRequest {
        HestonRequest {
            spot: 100.0,
            variance: 0.04,
            risk_free_rate: 0.05,
            dividend_yield: 0.02,
            kappa: 1.2,
            theta: 0.06,
            sigma: 0.3,
            rho: -0.5,
            horizon: 1.0,
            steps: 4,
            paths: 3,
            seed: 42,
            discretization: Discretization::QuadraticExponentialMartingale,
            terminal_only: false,
        }
    }

    #[test]
    fn heston_draw_order_layout_and_terminal_match_process() {
        let request = heston_request();
        let full = heston_paths(&request).unwrap();
        assert_eq!(full.len(), 30);
        assert_eq!(&full[..2], &[100.0, 0.04]);
        let process = HestonProcess::with_discretization(
            flat_curve(request.risk_free_rate),
            flat_curve(request.dividend_yield),
            make_quote_handle(request.spot).handle(),
            request.variance,
            request.kappa,
            request.theta,
            request.sigma,
            request.rho,
            request.discretization,
        );
        let mut draws = PseudoRandom::make_sequence_generator(1, request.seed).unwrap();
        let z0 = draws.next_sequence().value[0];
        let z1 = draws.next_sequence().value[0];
        let first = process
            .evolve(
                0.0,
                &Array::from([100.0, 0.04]),
                0.25,
                &Array::from([z0, z1]),
            )
            .unwrap();
        assert_eq!(&full[2..4], &[first[0], first[1]]);
        assert!((full[2] - 93.218_736_641_315_03).abs() < 1e-11);
        assert!((full[3] - 0.065_675_158_526_020_28).abs() < 1e-13);
        let terminal = heston_paths(&HestonRequest {
            terminal_only: true,
            ..request
        })
        .unwrap();
        for (path, row) in terminal.chunks_exact(2).enumerate() {
            assert_eq!(row, &full[path * 10 + 8..path * 10 + 10]);
        }
        assert_eq!(full, heston_paths(&heston_request()).unwrap());
        assert_ne!(
            full,
            heston_paths(&HestonRequest {
                seed: 43,
                ..heston_request()
            })
            .unwrap()
        );
    }

    #[test]
    fn heston_qe_high_psi_fixed_fixture() {
        let request = HestonRequest {
            variance: 0.01,
            kappa: 0.5,
            theta: 0.01,
            sigma: 0.2,
            horizon: 1.0,
            steps: 1,
            paths: 1,
            discretization: Discretization::QuadraticExponential,
            terminal_only: true,
            ..heston_request()
        };
        let result = heston_paths(&request).unwrap();
        assert!((result[0] - 96.553_174_001_572_44).abs() < 1e-11);
        assert!((result[1] - 0.018_076_059_597_846_472).abs() < 1e-13);
    }

    #[test]
    fn heston_zero_horizon_and_invalid_inputs() {
        let zero = HestonRequest {
            horizon: 0.0,
            ..heston_request()
        };
        for state in heston_paths(&zero).unwrap().chunks_exact(2) {
            assert_eq!(state, &[100.0, 0.04]);
        }
        let cases = [
            HestonRequest {
                steps: 0,
                ..heston_request()
            },
            HestonRequest {
                paths: 0,
                ..heston_request()
            },
            HestonRequest {
                seed: 0,
                ..heston_request()
            },
            HestonRequest {
                spot: 0.0,
                ..heston_request()
            },
            HestonRequest {
                variance: -0.1,
                ..heston_request()
            },
            HestonRequest {
                risk_free_rate: Real::NAN,
                ..heston_request()
            },
            HestonRequest {
                dividend_yield: Real::INFINITY,
                ..heston_request()
            },
            HestonRequest {
                kappa: 0.0,
                ..heston_request()
            },
            HestonRequest {
                theta: 0.0,
                ..heston_request()
            },
            HestonRequest {
                sigma: 0.0,
                ..heston_request()
            },
            HestonRequest {
                rho: 1.1,
                ..heston_request()
            },
            HestonRequest {
                horizon: -1.0,
                ..heston_request()
            },
            HestonRequest {
                steps: Size::MAX,
                ..heston_request()
            },
            HestonRequest {
                paths: Size::MAX,
                ..heston_request()
            },
            HestonRequest {
                discretization: Discretization::FullTruncation,
                ..heston_request()
            },
        ];
        for request in cases {
            assert!(heston_paths(&request).is_err());
        }
    }

    #[test]
    fn heston_variance_moment_and_discounted_spot() {
        let request = HestonRequest {
            variance: 0.09,
            theta: 0.09,
            kappa: 1.0,
            sigma: 0.3,
            rho: -0.8,
            horizon: 0.1,
            steps: 1,
            paths: 30_000,
            terminal_only: true,
            ..heston_request()
        };
        let result = heston_paths(&request).unwrap();
        let n = request.paths as Real;
        let mean_v = result.chunks_exact(2).map(|v| v[1]).sum::<Real>() / n;
        let var_v = result
            .chunks_exact(2)
            .map(|v| (v[1] - mean_v).powi(2))
            .sum::<Real>()
            / n;
        let ex = (-request.kappa * request.horizon).exp();
        let target_v = request.theta + (request.variance - request.theta) * ex;
        let target_var = request.variance * request.sigma.powi(2) * ex / request.kappa * (1.0 - ex)
            + request.theta * request.sigma.powi(2) / (2.0 * request.kappa) * (1.0 - ex).powi(2);
        assert!((mean_v - target_v).abs() < 6.0 * (target_var / n).sqrt());
        assert!((var_v - target_var).abs() < 0.05 * target_var);
        let mean_spot = result.chunks_exact(2).map(|v| v[0]).sum::<Real>() / n;
        let spot_var = result
            .chunks_exact(2)
            .map(|v| (v[0] - mean_spot).powi(2))
            .sum::<Real>()
            / n;
        let target_spot = request.spot
            * ((request.risk_free_rate - request.dividend_yield) * request.horizon).exp();
        assert!((mean_spot - target_spot).abs() < 6.0 * (spot_var / n).sqrt());
    }

    fn request() -> GbmRequest<'static> {
        GbmRequest {
            initial: &[100.0, 80.0],
            drift: &[0.05, -0.02],
            volatility: &[0.2, 0.3],
            correlation: &[1.0, 0.4, 0.4, 1.0],
            horizon: 1.0,
            steps: 4,
            paths: 3,
            seed: 42,
            terminal_only: false,
        }
    }

    #[test]
    fn gaussian_matches_native_sequence_across_batch_boundaries() {
        let draws = gaussian_draws(35, 42).unwrap();
        let mut native = PseudoRandom::make_sequence_generator(7, 42).unwrap();
        for chunk in draws.chunks_exact(7) {
            assert_eq!(chunk, native.next_sequence().value);
        }
        assert_eq!(draws, gaussian_draws(35, 42).unwrap());
        assert_ne!(draws, gaussian_draws(35, 43).unwrap());
        assert!(gaussian_draws(0, 42).unwrap().is_empty());
        assert!(gaussian_draws(1, 0).is_err());
        assert!(gaussian_draws(Size::MAX, 42).is_err());
    }

    #[test]
    fn correlated_paths_match_native_rng_and_analytic_two_asset_factor() {
        let r = request();
        let paths = gbm_paths(&r).unwrap();
        let mut native = PseudoRandom::make_sequence_generator(8, 42).unwrap();
        for path in paths.chunks_exact(10) {
            assert_eq!(&path[..2], r.initial);
            let z = &native.next_sequence().value;
            let mut expected = [100.0, 80.0];
            for step in 0..4 {
                let correlated = [
                    z[2 * step],
                    0.4 * z[2 * step] + (1.0_f64 - 0.16).sqrt() * z[2 * step + 1],
                ];
                for i in 0..2 {
                    let sigma = r.volatility[i];
                    expected[i] *= ((r.drift[i] - 0.5 * sigma * sigma) * 0.25
                        + sigma * 0.5 * correlated[i])
                        .exp();
                    assert_eq!(path[(step + 1) * 2 + i], expected[i]);
                }
            }
        }
        let terminal = gbm_paths(&GbmRequest {
            terminal_only: true,
            ..r
        })
        .unwrap();
        for (i, chunk) in terminal.chunks_exact(2).enumerate() {
            assert_eq!(chunk, &paths[i * 10 + 8..i * 10 + 10]);
        }
    }

    #[test]
    fn zero_horizon_and_volatility_are_exact() {
        let r = GbmRequest {
            horizon: 0.0,
            ..request()
        };
        for pair in gbm_paths(&r).unwrap().chunks_exact(2) {
            assert_eq!(pair, r.initial);
        }
        let r = GbmRequest {
            volatility: &[0.0, 0.0],
            ..request()
        };
        let result = gbm_paths(&r).unwrap();
        for path in result.chunks_exact(10) {
            for time in 0..=4 {
                for asset in 0..2 {
                    assert_eq!(
                        path[2 * time + asset],
                        r.initial[asset] * (r.drift[asset] * time as Real / 4.0).exp()
                    );
                }
            }
        }
        let changed_seed = gbm_paths(&GbmRequest { seed: 91, ..r }).unwrap();
        assert_eq!(result, changed_seed);
    }

    #[test]
    fn scalar_terminal_log_returns_have_theoretical_moments() {
        let r = GbmRequest {
            initial: &[100.0],
            drift: &[0.07],
            volatility: &[0.2],
            correlation: &[1.0],
            horizon: 2.0,
            steps: 4,
            paths: 50_000,
            seed: 1234,
            terminal_only: true,
        };
        let values = gbm_paths(&r).unwrap();
        let n = values.len() as Real;
        let mean = values.iter().map(|s| (s / 100.0).ln()).sum::<Real>() / n;
        let variance = values
            .iter()
            .map(|s| ((s / 100.0).ln() - mean).powi(2))
            .sum::<Real>()
            / n;
        // Six standard errors for the mean and variance of normal log returns.
        assert!((mean - 0.1).abs() < 6.0 * (0.08 / n).sqrt());
        assert!((variance - 0.08).abs() < 6.0 * 0.08 * (2.0 / n).sqrt());
    }

    #[test]
    fn correlated_log_returns_reproduce_requested_correlation() {
        let r = GbmRequest {
            paths: 50_000,
            steps: 1,
            terminal_only: true,
            ..request()
        };
        let values = gbm_paths(&r).unwrap();
        let n = r.paths as Real;
        let mut sum = [0.0; 2];
        let mut square = [0.0; 2];
        let mut cross = 0.0;
        for pair in values.chunks_exact(2) {
            let x = (pair[0] / 100.0).ln();
            let y = (pair[1] / 80.0).ln();
            sum[0] += x;
            sum[1] += y;
            square[0] += x * x;
            square[1] += y * y;
            cross += x * y;
        }
        let covariance = cross / n - sum[0] * sum[1] / (n * n);
        let variance = [
            square[0] / n - (sum[0] / n).powi(2),
            square[1] / n - (sum[1] / n).powi(2),
        ];
        let correlation = covariance / (variance[0] * variance[1]).sqrt();
        assert!((correlation - 0.4).abs() < 6.0 * (1.0 - 0.16) / n.sqrt());
    }

    #[test]
    fn rejects_invalid_parameters_and_dimensions() {
        let cases = [
            GbmRequest {
                steps: 0,
                ..request()
            },
            GbmRequest {
                paths: 0,
                ..request()
            },
            GbmRequest {
                seed: 0,
                ..request()
            },
            GbmRequest {
                horizon: -1.0,
                ..request()
            },
            GbmRequest {
                horizon: Real::NAN,
                ..request()
            },
            GbmRequest {
                initial: &[],
                ..request()
            },
            GbmRequest {
                initial: &[0.0, 1.0],
                ..request()
            },
            GbmRequest {
                initial: &[Real::INFINITY, 1.0],
                ..request()
            },
            GbmRequest {
                drift: &[0.0],
                ..request()
            },
            GbmRequest {
                drift: &[Real::NAN, 0.0],
                ..request()
            },
            GbmRequest {
                volatility: &[-0.1, 0.2],
                ..request()
            },
            GbmRequest {
                volatility: &[Real::INFINITY, 0.2],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0, 0.1, 0.2, 1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0, Real::NAN, Real::NAN, 1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[0.9, 0.0, 0.0, 1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0, 1.0, 1.0, 1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0, -1.0, -1.0, 1.0],
                ..request()
            },
            GbmRequest {
                correlation: &[1.0, 1.1, 1.1, 1.0],
                ..request()
            },
            GbmRequest {
                paths: Size::MAX,
                ..request()
            },
            GbmRequest {
                steps: Size::MAX,
                ..request()
            },
            GbmRequest {
                volatility: &[0.0, 0.0],
                drift: &[1e308, 0.0],
                ..request()
            },
        ];
        for r in cases {
            assert!(gbm_paths(&r).is_err());
        }
        // All entries are valid correlations, but this matrix is indefinite.
        assert!(correlation_factor(&[1.0, 0.9, 0.9, 0.9, 1.0, -0.9, 0.9, -0.9, 1.0], 3).is_err());
    }
}