libitofin 0.5.0

//! Sine and cosine integrals `Si` and `Ci`.
//!
//! Port of `ql/math/integrals/exponentialintegrals.*`. The real-valued sine
//! integral [`si`] and cosine integral [`ci`] use rational (minimax)
//! approximations for `|x| <= 4` and an asymptotic form built from the
//! auxiliary functions `f` and `g` beyond that. The complex-valued family
//! ([`ei`], [`e1`]) follows Pegoraro & Slusallek, "On the Evaluation of the
//! Complex-Valued Exponential Integral", combining a power series, an
//! asymptotic series and a continued fraction by region.

#![allow(clippy::excessive_precision)]

use crate::errors::QlResult;
use crate::math::comparison::sign;
use crate::types::{Complex, Real, Size};
use crate::{fail, require};

use std::f64::consts::PI;

/// The Euler-Mascheroni constant (QuantLib's `M_EULER_MASCHERONI`).
const EULER_MASCHERONI: Real = 0.5772156649015328606065120900824024;

/// Evaluates `coeffs[0] + coeffs[1] w + coeffs[2] w^2 + ...` by Horner's rule.
fn horner(w: Real, coeffs: &[Real]) -> Real {
    coeffs.iter().rev().fold(0.0, |acc, &c| acc * w + c)
}

// Auxiliary-function `f` numerator/denominator polynomials in w = 1/x^2.
const F_NUM: [Real; 11] = [
    1.0,
    7.44437068161936700618e2,
    1.96396372895146869801e5,
    2.37750310125431834034e7,
    1.43073403821274636888e9,
    4.33736238870432522765e10,
    6.40533830574022022911e11,
    4.20968180571076940208e12,
    1.00795182980368574617e13,
    4.94816688199951963482e12,
    -4.94701168645415959931e11,
];
const F_DEN: [Real; 10] = [
    1.0,
    7.46437068161927678031e2,
    1.97865247031583951450e5,
    2.41535670165126845144e7,
    1.47478952192985464958e9,
    4.58595115847765779830e10,
    7.08501308149515401563e11,
    5.06084464593475076774e12,
    1.43468549171581016479e13,
    1.11535493509914254097e13,
];

/// Auxiliary function `f` for the large-argument asymptotics of `Si`/`Ci`.
fn f_aux(x: Real) -> Real {
    let w = 1.0 / (x * x);
    horner(w, &F_NUM) / (x * horner(w, &F_DEN))
}

// Auxiliary-function `g` numerator (times w) / denominator polynomials.
const G_NUM: [Real; 11] = [
    1.0,
    8.1359520115168615e2,
    2.35239181626478200e5,
    3.12557570795778731e7,
    2.06297595146763354e9,
    6.83052205423625007e10,
    1.09049528450362786e12,
    7.57664583257834349e12,
    1.81004487464664575e13,
    6.43291613143049485e12,
    -1.36517137670871689e12,
];
const G_DEN: [Real; 10] = [
    1.0,
    8.19595201151451564e2,
    2.40036752835578777e5,
    3.26026661647090822e7,
    2.23355543278099360e9,
    7.87465017341829930e10,
    1.39866710696414565e12,
    1.17164723371736605e13,
    4.01839087307656620e13,
    3.99653257887490811e13,
];

/// Auxiliary function `g` for the large-argument asymptotics of `Si`/`Ci`.
fn g_aux(x: Real) -> Real {
    let w = 1.0 / (x * x);
    w * horner(w, &G_NUM) / horner(w, &G_DEN)
}

// Rational approximation of Si(x) for x <= 4, in w = x^2.
const SI_NUM: [Real; 8] = [
    1.0,
    -4.54393409816329991e-2,
    1.15457225751016682e-3,
    -1.41018536821330254e-5,
    9.43280809438713025e-8,
    -3.53201978997168357e-10,
    7.08240282274875911e-13,
    -6.05338212010422477e-16,
];
const SI_DEN: [Real; 7] = [
    1.0,
    1.01162145739225565e-2,
    4.99175116169755106e-5,
    1.55654986308745614e-7,
    3.28067571055789734e-10,
    4.5049097575386581e-13,
    3.21107051193712168e-16,
];

// Rational approximation of Ci(x) for x <= 4, in w = x^2. The numerator is
// multiplied by w (there is no constant term).
const CI_NUM: [Real; 7] = [
    -0.25,
    7.51851524438898291e-3,
    -1.27528342240267686e-4,
    1.05297363846239184e-6,
    -4.68889508144848019e-9,
    1.06480802891189243e-11,
    -9.93728488857585407e-15,
];
const CI_DEN: [Real; 8] = [
    1.0,
    1.1592605689110735e-2,
    6.72126800814254432e-5,
    2.55533277086129636e-7,
    6.97071295760958946e-10,
    1.38536352772778619e-12,
    1.89106054713059759e-15,
    1.39759616731376855e-18,
];

/// The sine integral `Si(x) = \int_0^x sin(t)/t dt`, defined for all real `x`
/// (odd: `Si(-x) = -Si(x)`).
pub fn si(x: Real) -> Real {
    if x < 0.0 {
        return -si(-x);
    }
    if x <= 4.0 {
        let w = x * x;
        return x * horner(w, &SI_NUM) / horner(w, &SI_DEN);
    }
    std::f64::consts::FRAC_PI_2 - f_aux(x) * x.cos() - g_aux(x) * x.sin()
}

/// The cosine integral `Ci(x) = \gamma + ln(x) + \int_0^x (cos(t)-1)/t dt`,
/// defined for `x > 0`.
///
/// For `x < 0` the cosine integral is complex (`Ci(x) = Ci(-x) + i pi`), so this
/// real-valued function returns `NaN` there, following the `f64` convention for
/// out-of-domain arguments; `Ci(0)` is `-inf`.
pub fn ci(x: Real) -> Real {
    if x < 0.0 {
        return Real::NAN;
    }
    if x <= 4.0 {
        let w = x * x;
        return EULER_MASCHERONI + x.ln() + w * horner(w, &CI_NUM) / horner(w, &CI_DEN);
    }
    f_aux(x) * x.sin() - g_aux(x) * x.cos()
}

/// `Ei(z) + acc`, where `acc` carries the branch adjustment used by [`e1`].
fn ei_with_acc(z: Complex, acc: Complex) -> QlResult<Complex> {
    if z.re == 0.0 && z.im == 0.0 {
        return Ok(Complex::new(Real::NEG_INFINITY, 0.0));
    }

    const DIST: Real = 4.5;
    const MAX_ERROR: Real = 5.0 * Real::EPSILON;

    let z_inf = (0.01 * Real::MAX).ln() + 100.0_f64.ln();
    let below_overflow_threshold = z.re < z_inf;
    require!(below_overflow_threshold, "argument error {z}");

    let z_asym = 2.0 - 1.035 * MAX_ERROR.ln();
    let abs_z = z.norm();

    let matches = |z1: Complex, z2: Complex| -> bool {
        let d = z1 - z2;
        d.re.abs() <= MAX_ERROR * z1.re.abs() && d.im.abs() <= MAX_ERROR * z1.im.abs()
    };

    if abs_z > 1.1 * z_asym {
        let mut ei = acc + Complex::new(0.0, sign(z.im) * PI);
        let mut s = z.exp() / z;
        let last = abs_z.floor() + 1.0;
        let mut i: Real = 1.0;
        while i <= last {
            if matches(ei + s, ei) {
                return Ok(ei + s);
            }
            ei += s;
            s *= i / z;
            i += 1.0;
        }
        fail!("series conversion issue for Ei({z})");
    }

    if abs_z > DIST && (z.re < 0.0 || z.im.abs() > DIST) {
        let mut cf = Complex::new(0.0, 0.0);
        for k in (1..=47u32).rev() {
            cf = -((k * k) as Real / (2.0 * k as Real + 1.0 - z + cf));
        }
        return Ok((acc + Complex::new(0.0, sign(z.im) * PI)) - z.exp() / (1.0 - z + cf));
    }

    let mut s = Complex::new(0.0, 0.0);
    let mut sn = z;
    let mut nn: Real = 1.0;
    let mut n: Size = 2;
    while n < 1000 && s + sn * nn != s {
        s += sn * nn;
        if (n & 1) != 0 {
            nn += 1.0 / (2.0 * (n / 2) as Real + 1.0);
        }
        sn *= -z / (2.0 * n as Real);
        n += 1;
    }
    require!(n < 1000, "series conversion issue for Ei({z})");

    let r = (EULER_MASCHERONI + acc) + z.ln() + (0.5 * z).exp() * s;
    if z.im != 0.0 {
        Ok(r)
    } else {
        Ok(Complex::new(r.re, acc.im))
    }
}

/// The complex exponential integral `Ei(z)` (principal branch).
///
/// `Ei(0)` is `-inf`; arguments with a real part beyond the overflow threshold
/// of `exp` are out of domain and return `Err`.
pub fn ei(z: Complex) -> QlResult<Complex> {
    ei_with_acc(z, Complex::new(0.0, 0.0))
}

/// The complex exponential integral `E1(z) = -Ei(-z)` on the principal branch.
pub fn e1(z: Complex) -> QlResult<Complex> {
    let acc = if z.im < 0.0 {
        Complex::new(0.0, -PI)
    } else if z.im > 0.0 || z.re < 0.0 {
        Complex::new(0.0, PI)
    } else {
        Complex::new(0.0, 0.0)
    };
    Ok(-(ei_with_acc(-z, acc)?))
}

/// The complex sine integral `Si(z)`, via a Taylor series for `|z| <= 0.2` and
/// [`e1`] beyond that.
pub fn si_complex(z: Complex) -> QlResult<Complex> {
    if z.norm() <= 0.2 {
        let mut s = Complex::new(0.0, 0.0);
        let mut nn = z;
        let mut k: Size = 2;
        while k < 100 && s != s + nn {
            s += nn;
            let kr = k as Real;
            nn *= -z * z / ((2.0 * kr - 2.0) * (2.0 * kr - 1.0) * (2.0 * kr - 1.0))
                * (2.0 * kr - 3.0);
            k += 1;
        }
        require!(k < 100, "series conversion issue for Si({z})");
        Ok(s)
    } else {
        let i = Complex::new(0.0, 1.0);
        let branch = if (z.re >= 0.0 && z.im >= 0.0) || (z.re > 0.0 && z.im < 0.0) {
            PI
        } else {
            -PI
        };
        Ok(0.5 * i * (e1(-i * z)? - e1(i * z)? - Complex::new(0.0, branch)))
    }
}

/// The complex cosine integral `Ci(z)` (principal branch), via [`e1`].
pub fn ci_complex(z: Complex) -> QlResult<Complex> {
    let i = Complex::new(0.0, 1.0);
    let mut acc = Complex::new(0.0, 0.0);
    if z.re < 0.0 && z.im >= 0.0 {
        acc.im = PI;
    } else if z.re <= 0.0 && z.im <= 0.0 {
        acc.im = -PI;
    }
    Ok(-0.5 * (e1(-i * z)? + e1(i * z)?) + acc)
}

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

    // (x, Si(x), Ci(x)) reference values from QuantLib's testRealSiCiIntegrals
    // (computed with Mathematica / mpmath).
    const DATA: [(Real, Real, Real); 17] = [
        (1e-12, 1e-12, -27.0538054510270153677),
        (0.1, 0.09994446110827695570, -1.7278683866572965838),
        (1.0, 0.9460830703671830149, 0.3374039229009681347),
        (1.9999, 1.6053675097543679041, 0.4230016343635392),
        (3.9999, 1.758222058430840841, -0.140965355646150101),
        (4.0001, 1.758184218306157867, -0.140998037827177150),
        (5.0, 1.5499312449446741373, -0.19002974965664387862),
        (7.0, 1.4545966142480935906, 0.076695278482184518383),
        (10.0, 1.6583475942188740493, -0.045456433004455372635),
        (15.0, 1.6181944437083687391, 0.046278677674360439604),
        (20.0, 1.5482417010434398402, 0.04441982084535331654),
        (24.9, 1.532210740207620024, -0.010788215638781789846),
        (25.1, 1.5311526281483412938, -0.0028719014454227088097),
        (30.0, 1.566756540030351111, -0.033032417282071143779),
        (40.0, 1.5869851193547845068, 0.019020007896208766962),
        (400.0, 1.5721148692738117518, -0.00212398883084634893),
        (4000.0, 1.5709788562309441985, -0.00017083030544201591130),
    ];

    #[test]
    fn matches_reference_si_and_ci() {
        let tol = 1e-12;
        for (x, si_ref, ci_ref) in DATA {
            assert!(
                (si(x) - si_ref).abs() < tol,
                "Si({x}) = {} vs {si_ref}",
                si(x)
            );
            assert!(
                (ci(x) - ci_ref).abs() < tol,
                "Ci({x}) = {} vs {ci_ref}",
                ci(x)
            );
            // Si is odd: Si(-x) = -Si(x).
            assert!(
                (si(-x) + si_ref).abs() < tol,
                "Si(-{x}) = {} vs {}",
                si(-x),
                -si_ref
            );
        }
    }

    #[test]
    fn ci_is_nan_below_zero_and_neg_infinity_at_zero() {
        assert!(ci(-1.0).is_nan());
        assert_eq!(ci(0.0), Real::NEG_INFINITY);
    }

    #[test]
    fn si_at_zero_is_zero() {
        assert_eq!(si(0.0), 0.0);
    }

    // Port of testExponentialIntegralLimits: QL_CHECK_CLOSE(a, b, t) is
    // BOOST_CHECK_CLOSE, whose tolerance is a percentage, hence the /100.
    fn assert_close(label: &str, value: Real, reference: Real, rel_tol_percent: Real) {
        let rel = rel_tol_percent / 100.0;
        assert!(
            value == reference
                || ((value - reference).abs() <= rel * value.abs()
                    && (value - reference).abs() <= rel * reference.abs()),
            "{label}: {value} vs {reference} (rel tol {rel})"
        );
    }

    #[test]
    fn ei_limit_for_large_arguments() {
        let tol = 1000.0 * Real::EPSILON;
        let large_value = 0.75 * (0.1 * Real::MAX).ln();
        let expected_real = large_value.exp() / large_value;

        let pos_imag = ei(Complex::new(large_value, Real::MIN_POSITIVE)).unwrap();
        assert_close("Ei imag (+0 side)", pos_imag.im, PI, tol);
        assert_close(
            "Ei real (+0 side)",
            pos_imag.re,
            expected_real,
            1e3 / large_value,
        );

        let neg_imag = ei(Complex::new(large_value, -Real::MIN_POSITIVE)).unwrap();
        assert_close("Ei imag (-0 side)", neg_imag.im, -PI, tol);
        assert_close(
            "Ei real (-0 side)",
            neg_imag.re,
            expected_real,
            1e3 / large_value,
        );

        let zero_imag = ei(Complex::new(large_value, 0.0)).unwrap();
        assert_eq!(zero_imag.im, 0.0);
    }

    #[test]
    fn ei_huge_imaginary_argument_converges() {
        let value = ei(Complex::new(0.0, 1e20)).unwrap();
        assert_close("Ei huge-imag imag", value.im, PI, 1000.0 * Real::EPSILON);
        assert!(value.re.abs() < 1e-15, "Ei huge-imag real: {}", value.re);
    }

    #[test]
    fn ei_at_zero_is_negative_infinity() {
        let value = ei(Complex::new(0.0, 0.0)).unwrap();
        assert_eq!(value, Complex::new(Real::NEG_INFINITY, 0.0));
    }

    #[test]
    fn ei_is_out_of_domain_beyond_exp_overflow() {
        assert!(ei(Complex::new(710.0, 0.0)).is_err());
    }

    #[test]
    fn ei_small_circle_limit_is_euler_plus_log() {
        let tol = 1000.0 * Real::EPSILON;
        let small_r = Real::EPSILON * Real::EPSILON;
        for x in -100..100 {
            let phi = x as Real / 100.0 * PI;
            let z = Complex::from_polar(small_r, phi);
            let value = ei(z).unwrap();
            let limit = EULER_MASCHERONI + z.ln();
            assert_close("Ei small-circle real", value.re, limit.re, tol);
            assert_close("Ei small-circle imag", value.im, limit.im, tol);
        }
    }

    #[test]
    fn ei_large_circle_limit_is_signed_pi() {
        let tol = 1000.0 * Real::EPSILON;
        let large_r = 0.75 * (0.1 * Real::MAX).ln();
        for x in -10..10 {
            let phi = x as Real / 10.0 * PI;
            if phi.abs() > 0.5 * PI {
                let z = Complex::from_polar(large_r, phi);
                let value = ei(z).unwrap();
                let limit_imag = sign(z.im) * PI;
                assert!(
                    crate::math::comparison::close_enough(value.re, 0.0),
                    "Ei large-circle real: {} not close enough to 0",
                    value.re
                );
                assert_close("Ei large-circle imag", value.im, limit_imag, tol);
            }
        }
    }

    // Port of testExponentialIntegral: rows are (z.re, z.im, Si.re, Si.im,
    // Ci.re, Ci.im, Ei.re, Ei.im, E1.re, E1.im), reference values calculated
    // with Mathematica or Python/mpmath.
    #[allow(clippy::approx_constant)]
    #[rustfmt::skip]
    const COMPLEX_DATA: [[Real; 10]; 99] = [
        [1e-10, 0.0, 1.0e-10, 0.0, -22.4486352650389, 0.0, -22.4486352649389, 0.0, 22.4486352651389, 0.0],
        [7.0710678118655e-11, 7.0710678118655e-11, 7.0710678118655e-11, 7.0710678118655e-11, -22.4486352650389, 0.785398163397448, -22.4486352649682, 0.785398163468159, 22.4486352651096, -0.785398163326738],
        [3.0901699437495e-11, 9.5105651629515e-11, 3.0901699437495e-11, 9.5105651629515e-11, -22.4486352650389, 1.25663706143591, -22.448635265008, 1.25663706153102, 22.4486352650698, -1.25663706134081],
        [0.0, 1e-10, 0.0, 1.0e-10, -22.4486352650389, 1.5707963267949, -22.4486352650389, 1.5707963268949, 22.4486352650389, -1.5707963266949],
        [0.0, 1e-10, 0.0, 1.0e-10, -22.4486352650389, 1.5707963267949, -22.4486352650389, 1.5707963268949, 22.4486352650389, -1.5707963266949],
        [-8.0901699437495e-11, 5.8778525229247e-11, -8.0901699437495e-11, 5.8778525229247e-11, -22.4486352650389, 2.51327412287184, -22.4486352651198, 2.51327412293062, 22.448635264958, -2.51327412281306],
        [-1e-10, 0.0, -1.0e-10, 0.0, -22.4486352650389, 3.14159265358979, -22.4486352651389, 0.0, 22.4486352649389, -3.14159265358979],
        [-8.0901699437495e-11, -5.8778525229247e-11, -8.0901699437495e-11, -5.8778525229247e-11, -22.4486352650389, -2.51327412287184, -22.4486352651198, -2.51327412293062, 22.448635264958, 2.51327412281306],
        [0.0, -1e-10, 0.0, -1.0e-10, -22.4486352650389, -1.5707963267949, -22.4486352650389, -1.5707963268949, 22.4486352650389, 1.5707963266949],
        [3.0901699437495e-11, -9.5105651629515e-11, 3.0901699437495e-11, -9.5105651629515e-11, -22.4486352650389, -1.25663706143591, -22.448635265008, -1.25663706153102, 22.4486352650698, 1.25663706134081],
        [9.8768834059514e-11, -1.5643446504023002e-11, 9.8768834059514e-11, -1.5643446504023e-11, -22.4486352650389, -0.157079632679488, -22.4486352649402, -0.157079632695132, 22.4486352651377, 0.157079632663845],
        [0.15, 0.0, 0.149812626514082, 0.0, -1.32552404918277, 0.0, -1.16408641729839, 0.0, 1.46446167052028, 0.0],
        [0.1060660171779825, 0.1060660171779825, 0.106198510172016, 0.105933345197561, -1.31990959342105, 0.779773166034167, -1.21397624822349, 0.897221670932746, 1.42584293861861, -0.684824650588713],
        [0.0463525491562425, 0.14265847744427249, 0.0465043664443717, 0.1427686871506, -1.31535197062462, 1.25332575154654, -1.27825242518864, 1.40248660838809, 1.37065439517488, -1.11739007291224],
        [0.0, 0.15, 0.0, 0.150187626610941, -1.31427404390933, 1.5707963267949, -1.32552404918277, 1.72060895330898, 1.32552404918277, -1.42098370028081],
        [0.0, 0.15, 0.0, 0.150187626610941, -1.31427404390933, 1.5707963267949, -1.32552404918277, 1.72060895330898, 1.32552404918277, -1.42098370028081],
        [-0.1213525491562425, 0.0881677878438705, -0.121410363295163, 0.0879894647931175, -1.32164680474487, 2.51862071457814, -1.43946484971679, 2.59626744276408, 1.19687588593211, -2.41957522097486],
        [-0.15, 0.0, -0.149812626514082, 0.0, -1.32552404918277, 3.14159265358979, -1.46446167052028, 0.0, 1.16408641729839, -3.14159265358979],
        [-0.1213525491562425, -0.0881677878438705, -0.121410363295163, -0.0879894647931175, -1.32164680474487, -2.51862071457814, -1.43946484971679, -2.59626744276408, 1.19687588593211, 2.41957522097486],
        [0.0, -0.15, 0.0, -0.150187626610941, -1.31427404390933, -1.5707963267949, -1.32552404918277, -1.72060895330898, 1.32552404918277, 1.42098370028081],
        [0.0463525491562425, -0.14265847744427249, 0.0465043664443717, -0.1427686871506, -1.31535197062462, -1.25332575154654, -1.27825242518864, -1.40248660838809, 1.37065439517488, 1.11739007291224],
        [0.148153251089271, -0.0234651697560345, 0.147986276837203, -0.0233801359873959, -1.32524974813753, -0.155344509602526, -1.16622995490181, -0.182371337566645, 1.46287076355731, 0.135270572544445],
        [0.25, 0.0, 0.249133570319757, 0.0, -0.824663062580946, 0.0, -0.542543264661914, 0.0, 1.04428263444374, 0.0],
        [0.1767766952966375, 0.1767766952966375, 0.177389351153991, 0.17616173766105, -0.809119386275216, 0.769773219911456, -0.632957648614166, 0.978412458037432, 0.985281123936265, -0.623633755729451],
        [0.0772542485937375, 0.2377641290737875, 0.0779581492943877, 0.238274358309521, -0.796425249249655, 1.24741416450428, -0.745153392294084, 1.50303646097033, 0.898260598498369, -1.02852866129867],
        [0.0, 0.25, 0.0, 0.250869684890912, -0.793412949552826, 1.5707963267949, -0.824663062580946, 1.81992989711465, 0.824663062580946, -1.32166275647514],
        [0.0, 0.25, 0.0, 0.250869684890912, -0.793412949552826, 1.5707963267949, -0.824663062580946, 1.81992989711465, 0.824663062580946, -1.32166275647514],
        [-0.2022542485937375, 0.1469463130731175, -0.20252086544385, 0.146120744825161, -0.813939960005834, 2.52811043072268, -1.00626764691037, 2.64616186234439, 0.60229889383601, -2.35061809970499],
        [-0.25, 0.0, -0.249133570319757, 0.0, -0.824663062580946, 3.14159265358979, -1.04428263444374, 0.0, 0.542543264661914, -3.14159265358979],
        [-0.2022542485937375, -0.1469463130731175, -0.20252086544385, -0.146120744825161, -0.813939960005834, -2.52811043072268, -1.00626764691037, -2.64616186234439, 0.60229889383601, 2.35061809970499],
        [0.0, -0.25, 0.0, -0.250869684890912, -0.793412949552826, -1.5707963267949, -0.824663062580946, -1.81992989711465, 0.824663062580946, 1.32166275647514],
        [0.0772542485937375, -0.2377641290737875, 0.0779581492943877, -0.238274358309521, -0.796425249249655, -1.24741416450428, -0.745153392294084, -1.50303646097033, 0.898260598498369, 1.02852866129867],
        [0.246922085148785, -0.0391086162600575, 0.24614979209014, -0.0387156766342252, -0.823906068503191, -0.152275113509673, -0.546488805945054, -0.201435843693654, 1.04188216592042, 0.122428128357486],
        [1.0, 0.0, 0.946083070367183, 0.0, 0.337403922900968, 0.0, 1.89511781635594, 0.0, 0.21938393439552, 0.0],
        [0.70710678118655, 0.70710678118655, 0.745192155353662, 0.666664817419508, 0.566802098259312, 0.535629617322428, 1.23346691567882, 1.78035886482613, 0.0998627191601961, -0.289974554118806],
        [0.30901699437495, 0.95105651629515, 0.355652074843551, 0.983694298574337, 0.782614772996823, 1.09956193553216, 0.643964830804846, 2.31231301720838, -0.112533957890793, -0.475476714030747],
        [0.0, 1.0, 0.0, 1.05725087537573, 0.837866940980208, 1.5707963267949, 0.337403922900968, 2.51687939716208, -0.337403922900968, -0.624713256427714],
        [0.0, 1.0, 0.0, 1.05725087537573, 0.837866940980208, 1.5707963267949, 0.337403922900968, 2.51687939716208, -0.337403922900968, -0.624713256427714],
        [-0.80901699437495, 0.58778525229247, -0.824526943360603, 0.5349755552469, 0.491722358913221, 2.74478237579885, -0.14431784116889, 2.91012082986304, -1.43603057378731, -1.62893165104155],
        [-1.0, 0.0, -0.946083070367183, 0.0, 0.337403922900968, 3.14159265358979, -0.21938393439552, 0.0, -1.89511781635594, -3.14159265358979],
        [-0.80901699437495, -0.58778525229247, -0.824526943360603, -0.5349755552469, 0.491722358913221, -2.74478237579885, -0.14431784116889, -2.91012082986304, -1.43603057378731, 1.62893165104155],
        [0.0, -1.0, 0.0, -1.05725087537573, 0.837866940980208, -1.5707963267949, 0.337403922900968, -2.51687939716208, -0.337403922900968, 0.624713256427714],
        [0.30901699437495, -0.95105651629515, 0.355652074843551, -0.983694298574337, 0.782614772996823, -1.09956193553216, 0.643964830804846, -2.31231301720838, -0.112533957890793, 0.475476714030747],
        [0.98768834059514, -0.15643446504023, 0.939353669480516, -0.132366326809511, 0.347743692745538, -0.0857637957494435, 1.86192420379474, -0.4235071237, 0.214836056406461, 0.0577866622153682],
        [5.0, 0.0, 1.54993124494467, 0.0, -0.190029749656644, 0.0, 40.1852753558032, 0.0, 0.00114829559127533, 0.0],
        [3.53553390593275, 3.53553390593275, 3.68715086115432, -3.15718137390906, -3.15476810467167, -2.11185029092794, -6.31194947858072, 7.36979747887716, -0.00241326923739065, 0.00450424343148012],
        [1.5450849718747501, 4.75528258147575, 14.299679516973, 6.85221185491562, 6.85257226323722, -12.7303117750282, -0.931350039879264, 2.99045284011251, 0.0356665739529384, 0.0160488285537158],
        [0.0, 5.0, 0.0, 20.0932118256972, 20.0920635301059, 1.5707963267949, -0.190029749656644, 3.12072757173957, 0.190029749656644, -0.0208650818502225],
        [0.0, 5.0, 0.0, 20.0932118256972, 20.0920635301059, 1.5707963267949, -0.190029749656644, 3.12072757173957, 0.190029749656644, -0.0208650818502225],
        [-4.04508497187475, 2.93892626146235, -2.0577013528011, -1.96223940975232, -1.9637046590567, 3.61921566552724, 0.00286020292932927, 3.14261835694337, 6.84905720502975, 11.1883945116728],
        [-5.0, 0.0, -1.54993124494467, 0.0, -0.190029749656644, 3.14159265358979, -0.00114829559127533, 0.0, -40.1852753558032, -3.14159265358979],
        [-4.04508497187475, -2.93892626146235, -2.0577013528011, 1.96223940975232, -1.9637046590567, -3.61921566552724, 0.00286020292932927, -3.14261835694337, 6.84905720502975, -11.1883945116728],
        [0.0, -5.0, 0.0, -20.0932118256972, 20.0920635301059, -1.5707963267949, -0.190029749656644, -3.12072757173957, 0.190029749656644, 0.0208650818502225],
        [1.5450849718747501, -4.75528258147575, 14.299679516973, -6.85221185491562, 6.85257226323722, 12.7303117750282, -0.931350039879264, -2.99045284011251, 0.0356665739529384, -0.0160488285537158],
        [4.9384417029757, -0.7821723252011501, 1.53351371140353, 0.167535111630988, -0.252671967618136, -0.0455545136665558, 31.7637646606649, -20.6127722347705, 0.000742118122850436, 0.000971589948194675],
        [10.0, 0.0, 1.65834759421887, 0.0, -0.0454564330044554, 0.0, 2492.22897624188, 0.0, 4.15696892968532e-6, 0.0],
        [7.0710678118655, 7.0710678118655, -3.77451753034182, 62.6425755592338, 62.6425711229056, 5.34523470197841, 125.285146682139, -7.54895590552534, 4.43632828562146e-6, -7.91551583068017e-5],
        [3.0901699437495003, 9.5105651629515, 303.07292777526, -690.037761260879, -690.037754650298, -301.502129842997, -0.659900725018632, 5.27667742385125, -0.00134856502993308, 0.00415958644984393],
        [0.0, 10.0, 0.0, 1246.11449019942, 1246.11448604245, 1.5707963267949, -0.0454564330044554, 3.22914392101377, 0.0454564330044554, 0.0875512674239774],
        [0.0, 10.0, 0.0, 1246.11449019942, 1246.11448604245, 1.5707963267949, -0.0454564330044554, 3.22914392101377, 0.0454564330044554, 0.0875512674239774],
        [-8.0901699437495, 5.8778525229247, -14.6236949578037, 13.4643508624518, 13.4645870261785, 16.1946084513107, -2.79815608075126e-5, 3.14158769865141, -157.085481478947, -317.2439811058],
        [-10.0, 0.0, -1.65834759421887, 0.0, -0.0454564330044554, 3.14159265358979, -4.15696892968532e-6, 0.0, -2492.22897624188, -3.14159265358979],
        [-8.0901699437495, -5.8778525229247, -14.6236949578037, -13.4643508624518, 13.4645870261785, -16.1946084513107, -2.79815608075126e-5, -3.14158769865141, -157.085481478947, 317.2439811058],
        [0.0, -10.0, 0.0, -1246.11449019942, 1246.11448604245, -1.5707963267949, -0.0454564330044554, -3.22914392101377, 0.0454564330044554, -0.0875512674239774],
        [3.0901699437495003, -9.5105651629515, 303.07292777526, 690.037761260879, -690.037754650298, 301.502129842997, -0.659900725018632, -5.27667742385125, -0.00134856502993308, -0.00415958644984393],
        [9.8768834059514, -1.5643446504023002, 1.78956084261706, 0.114701769782499, -0.118816490702582, 0.198823504802007, 411.904076239608, -2157.22483235914, -6.48699583272709e-7, 4.66032253043785e-6],
        [25.0, 0.0, 1.53148255099996, 0.0, -0.00684859717970259, 0.0, 3005950906.52555, 0.0, 5.34889975534022e-13, 0.0],
        [17.67766952966375, 17.67766952966375, -894423.548678786, -396595.979622699, -396595.9796227, 894425.119475113, -793191.959245399, -1788847.09735757, 7.48981460647877e-10, 3.27816276287981e-10],
        [7.72542485937375, 23.77641290737875, 395787595.545024, 194501516.12134, 194501516.12134, -395787593.974227, -80.7948153607822, -39.8888851700048, 1.72503667797818e-5, 2.36415887840135e-6],
        [0.0, 25.0, 0.0, 1502975453.26277, 1502975453.26277, 1.5707963267949, -0.00684859717970259, 3.10227887779486, 0.00684859717970259, -0.0393137757949353],
        [0.0, 25.0, 0.0, 1502975453.26277, 1502975453.26277, 1.5707963267949, -0.00684859717970259, 3.10227887779486, 0.00684859717970259, -0.0393137757949353],
        [-20.22542485937375, 14.69463130731175, -19129.3494470458, 45406.0213041107, 45406.0213041213, 19130.9202433848, 5.85665949258649e-11, 3.14159265356458, -2432061.38760638, 25010638.0968068],
        [-25.0, 0.0, -1.53148255099996, 0.0, -0.00684859717970259, 3.14159265358979, -5.34889975534022e-13, 0.0, -3005950906.52555, -3.14159265358979],
        [-20.22542485937375, -14.69463130731175, -19129.3494470458, -45406.0213041107, 45406.0213041213, -19130.9202433848, 5.85665949258649e-11, -3.14159265356458, -2432061.38760638, -25010638.0968068],
        [0.0, -25.0, 0.0, -1502975453.26277, 1502975453.26277, -1.5707963267949, -0.00684859717970259, -3.10227887779486, 0.00684859717970259, 0.0393137757949353],
        [7.72542485937375, -23.77641290737875, 395787595.545024, -194501516.12134, 194501516.12134, 395787593.974227, -80.7948153607822, 39.8888851700048, 1.72503667797818e-5, -2.36415887840135e-6],
        [24.6922085148785, -3.91086162600575, 0.61973692887531, 0.318459426938049, -0.318931296543192, -0.950420524151913, -1816162045.63054, 1255955799.5082, -4.40593065675657e-13, -5.79490191675286e-13],
        [50.0, 0.0, 1.55161707248594, 0.0, -0.00562838632411631, 0.0, 1.05856368971317e+20, 0.0, 3.78326402955046e-24, 0.0],
        [35.3553390593275, 35.3553390593275, 53807668130.5995, -22948660925283.2, -22948660925283.2, -53807668129.0287, -45897321850566.4, 107615336261.199, -9.9766761181828e-21, 8.71502630154959e-18],
        [15.4508497187475, 47.5528258147575, 2.49903843573354e+18, -3.83240358282137e+18, -3.83240358282137e+18, -2.49903843573354e+18, -68343.3715391731, 77339.6040605891, 3.71621275609622e-10, 3.85406628982992e-9],
        [0.0, 50.0, 0.0, 5.29281844856585e+19, 5.29281844856585e+19, 1.5707963267949, -0.00562838632411631, 3.12241339928083, 0.00562838632411631, -0.0191792543089607],
        [0.0, 50.0, 0.0, 5.29281844856585e+19, 5.29281844856585e+19, 1.5707963267949, -0.00562838632411631, 3.12241339928083, 0.00562838632411631, -0.0191792543089607],
        [-40.4508497187475, 29.3892626146235, -57258797567.9644, -12906669326.6389, -12906669326.6389, 57258797569.5352, -8.55226617604501e-21, 3.14159265358979, 6.68228261723918e+15, -3.43017053184612e+15],
        [-50.0, 0.0, -1.55161707248594, 0.0, -0.00562838632411631, 3.14159265358979, -3.78326402955046e-24, 0.0, -1.05856368971317e+20, -3.14159265358979],
        [-40.4508497187475, -29.3892626146235, -57258797567.9644, 12906669326.6389, -12906669326.6389, -57258797569.5352, -8.55226617604501e-21, -3.14159265358979, 6.68228261723918e+15, 3.43017053184612e+15],
        [0.0, -50.0, 0.0, -5.29281844856585e+19, 5.29281844856585e+19, -1.5707963267949, -0.00562838632411631, -3.12241339928083, 0.00562838632411631, 0.0191792543089607],
        [15.4508497187475, -47.5528258147575, 2.49903843573354e+18, 3.83240358282137e+18, -3.83240358282137e+18, 2.49903843573354e+18, -68343.3715391731, -77339.6040605891, 3.71621275609622e-10, -3.85406628982992e-9],
        [49.384417029757, -7.8217232520115, -16.8292457944994, 16.9326906903424, -16.9326976506474, -18.4000381995002, 1.09489979806082e+19, -5.61228684199658e+19, -8.51344869310291e-25, 6.95142343223447e-24],
        [700.0, 0.0, 1.57199393223749, 0.0, 0.000778810012739756, 0.0, 1.45097873605256e+301, 0.0, 1.40651876623403e-307, 0.0],
        [494.974746830585, 494.974746830585, -5.39480977313549e+211, -3.7907051625115e+211, -3.7907051625115e+211, 5.39480977313549e+211, -7.58141032502299e+211, -1.0789619546271e+212, 1.26627531288803e-218, 8.89746644202181e-219],
        [216.311896062465, 665.7395614066049, 6.68861022474796e+285, -6.86204916856497e+285, -6.86204916856497e+285, -6.68861022474796e+285, 4.35129688126332e+89, -1.25283433405018e+91, 9.10599247691995e-98, -1.3494793845188e-97],
        [0.0, 700.0, 0.0, 7.2548936802628e+300, 7.2548936802628e+300, 1.5707963267949, 0.000778810012739756, 3.14279025903239, -0.000778810012739756, 0.00119760544259495],
        [0.0, 700.0, 0.0, 7.2548936802628e+300, 7.2548936802628e+300, 1.5707963267949, 0.000778810012739756, 3.14279025903239, -0.000778810012739756, 0.00119760544259495],
        [-566.311896062465, 411.449676604729, 4.13964135191794e+174, 3.47943069430311e+175, 3.47943069430311e+175, -4.13964135191794e+174, 1.39494929258574e-249, 3.14159265358979, 9.43022777090499e+242, 8.40743888884655e+242],
        [-700.0, 0.0, -1.57199393223749, 0.0, 0.000778810012739756, 3.14159265358979, -1.40651876623403e-307, 0.0, -1.45097873605256e+301, -3.14159265358979],
        [-566.311896062465, -411.449676604729, 4.13964135191794e+174, -3.47943069430311e+175, 3.47943069430311e+175, 4.13964135191794e+174, 1.39494929258574e-249, -3.14159265358979, 9.43022777090499e+242, -8.40743888884655e+242],
        [0.0, -700.0, 0.0, -7.2548936802628e+300, 7.2548936802628e+300, -1.5707963267949, 0.000778810012739756, -3.14279025903239, -0.000778810012739756, -0.00119760544259495],
        [216.311896062465, -665.7395614066049, 6.68861022474796e+285, 6.86204916856497e+285, -6.86204916856497e+285, 6.68861022474796e+285, 4.35129688126332e+89, 1.25283433405018e+91, 9.10599247691995e-98, 1.3494793845188e-97],
        [691.381838416598, -109.50412552816101, -2.38570018769502e+44, -9.72638025849046e+43, 9.72638025849046e+43, -2.38570018769502e+44, -2.15172979114587e+297, -1.50043260461905e+297, -7.44435180959991e-304, 2.26013762375079e-304],
    ];

    fn assert_complex_close(
        name: &str,
        z: Complex,
        value: Complex,
        reference: Complex,
        diff: Real,
        diff_tol: Real,
        tol: Real,
    ) {
        assert!(
            diff <= diff_tol
                && !diff.is_nan()
                && (reference.re.abs() >= tol || value.re.abs() <= tol)
                && (reference.im.abs() >= tol || value.im.abs() <= tol),
            "{name}({z}) = {value} vs {reference} (diff {diff}, tol {diff_tol})"
        );
    }

    #[test]
    fn matches_reference_complex_exponential_integrals() {
        let tol = 100.0 * Real::EPSILON;
        for row in &COMPLEX_DATA {
            let x = row[0];
            let y = if row[1].abs() < 1e-12 { 0.0 } else { row[1] };
            let z = Complex::new(x, y);

            let si = si_complex(z).unwrap();
            let mut reference = Complex::new(row[2], row[3]);
            let mut diff = (si - reference).norm() / reference.norm();
            assert_complex_close("Si", z, si, reference, diff, tol, tol);

            let ci = ci_complex(z).unwrap();
            reference = Complex::new(row[4], row[5]);
            diff = (ci - reference)
                .norm()
                .min((ci - reference).norm() / reference.norm());
            assert_complex_close("Ci", z, ci, reference, diff, tol, tol);

            let ei_value = ei(z).unwrap();
            reference = Complex::new(row[6], row[7]);
            diff = (ei_value - reference).norm() / reference.norm();
            assert_complex_close("Ei", z, ei_value, reference, diff, tol, tol);

            let e1_value = e1(z).unwrap();
            reference = Complex::new(row[8], row[9]);
            diff = (e1_value - reference).norm() / reference.norm();
            assert_complex_close("E1", z, e1_value, reference, diff, 10.0 * tol, tol);
        }
    }
}