spec_math 0.1.6

/*
* Cephes Math Library Release 2.1:  January, 1989
* Copyright 1984, 1987, 1989 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#![allow(clippy::excessive_precision)]

const SN: [f64; 6] = [
    -8.39167827910303881427E-11,
    4.62591714427012837309E-8,
    -9.75759303843632795789E-6,
    9.76945438170435310816E-4,
    -4.13470316229406538752E-2,
    1.00000000000000000302E0,
];

const SD: [f64; 6] = [
    2.03269266195951942049E-12,
    1.27997891179943299903E-9,
    4.41827842801218905784E-7,
    9.96412122043875552487E-5,
    1.42085239326149893930E-2,
    9.99999999999999996984E-1,
];

const CN: [f64; 6] = [
    2.02524002389102268789E-11,
    -1.35249504915790756375E-8,
    3.59325051419993077021E-6,
    -4.74007206873407909465E-4,
    2.89159652607555242092E-2,
    -1.00000000000000000080E0,
];

const CD: [f64; 6] = [
    4.07746040061880559506E-12,
    3.06780997581887812692E-9,
    1.23210355685883423679E-6,
    3.17442024775032769882E-4,
    5.10028056236446052392E-2,
    4.00000000000000000080E0,
];

const FN4: [f64; 7] = [
    4.23612862892216586994E0,
    5.45937717161812843388E0,
    1.62083287701538329132E0,
    1.67006611831323023771E-1,
    6.81020132472518137426E-3,
    1.08936580650328664411E-4,
    5.48900223421373614008E-7,
];

const FD4: [f64; 7] = [
    /*  1.00000000000000000000E0, */
    8.16496634205391016773E0,
    7.30828822505564552187E0,
    1.86792257950184183883E0,
    1.78792052963149907262E-1,
    7.01710668322789753610E-3,
    1.10034357153915731354E-4,
    5.48900252756255700982E-7,
];

const FN8: [f64; 9] = [
    4.55880873470465315206E-1,
    7.13715274100146711374E-1,
    1.60300158222319456320E-1,
    1.16064229408124407915E-2,
    3.49556442447859055605E-4,
    4.86215430826454749482E-6,
    3.20092790091004902806E-8,
    9.41779576128512936592E-11,
    9.70507110881952024631E-14,
];

const FD8: [f64; 8] = [
    /*  1.00000000000000000000E0, */
    9.17463611873684053703E-1,
    1.78685545332074536321E-1,
    1.22253594771971293032E-2,
    3.58696481881851580297E-4,
    4.92435064317881464393E-6,
    3.21956939101046018377E-8,
    9.43720590350276732376E-11,
    9.70507110881952025725E-14,
];

const GN4: [f64; 8] = [
    8.71001698973114191777E-2,
    6.11379109952219284151E-1,
    3.97180296392337498885E-1,
    7.48527737628469092119E-2,
    5.38868681462177273157E-3,
    1.61999794598934024525E-4,
    1.97963874140963632189E-6,
    7.82579040744090311069E-9,
];

const GD4: [f64; 7] = [
    /*  1.00000000000000000000E0, */
    1.64402202413355338886E0,
    6.66296701268987968381E-1,
    9.88771761277688796203E-2,
    6.22396345441768420760E-3,
    1.73221081474177119497E-4,
    2.02659182086343991969E-6,
    7.82579218933534490868E-9,
];

const GN8: [f64; 9] = [
    6.97359953443276214934E-1,
    3.30410979305632063225E-1,
    3.84878767649974295920E-2,
    1.71718239052347903558E-3,
    3.48941165502279436777E-5,
    3.47131167084116673800E-7,
    1.70404452782044526189E-9,
    3.85945925430276600453E-12,
    3.14040098946363334640E-15,
];

const GD8: [f64; 9] = [
    /*  1.00000000000000000000E0, */
    1.68548898811011640017E0,
    4.87852258695304967486E-1,
    4.67913194259625806320E-2,
    1.90284426674399523638E-3,
    3.68475504442561108162E-5,
    3.57043223443740838771E-7,
    1.72693748966316146736E-9,
    3.87830166023954706752E-12,
    3.14040098946363335242E-15,
];

use crate::cephes64::consts::{M_PI_2, EULER};
use crate::cephes64::polevl::{polevl, p1evl};

// $$\mathrm{Ci}(x) = \mathrm{EUL} + \ln(x) + \int_0^x{\frac{\cos(t)-1}{t}\,dt}$$
// $$\mathrm{Si}(x) = \int_0^x{\frac{\sin(t)}{t}\,dt}$$

pub fn sici(x: f64) -> (f64, f64) // si, ci
{
    //! Sine and cosine integrals
    //!
    //! ## DESCRIPTION:
    //!
    //! Evaluates the integrals
    //!
    #![doc=include_str!("ci.svg")]
    //!
    #![doc=include_str!("si.svg")]
    //!
    //! where eul = 0.57721566490153286061 is Euler's constant.
    //! The integrals are approximated by rational functions.
    //! For x > 8 auxiliary functions f(x) and g(x) are employed
    //! such that
    //!
    //! `Ci(x) = f(x) sin(x) - g(x) cos(x)`
    //!
    //! `Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)`
    //!
    //! ## ACCURACY:
    //!
    //! Absolute error, except relative when > 1:
    //!<table>
    //! <tr>
    //!     <th>Arithmetic</th>
    //!     <th>Function</th>
    //!     <th>Domain</th>
    //!     <th># Trials</th>
    //!     <th>Peak</th>
    //!     <th>RMS</th>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>Si(x)</td>
    //!     <td>0, 50</td>
    //!     <td>30000</td>
    //!     <td>4.4e-16</td>
    //!     <td>7.3e-17</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>Ci(x)</td>
    //!     <td>0, 50</td>
    //!     <td>30000</td>
    //!     <td>6.9e-16</td>
    //!     <td>5.1e-17</td>
    //! </tr>
    //!</table>

    let mut x = x;

    let sign = if x < 0.0 {
        x = -x;
        true
    } else {
        false
    };


    if x == 0.0 {
        return (0.0, -f64::INFINITY);
    }


    let (si, ci) = if x > 1.0e9 { 
        // This case is imprecise since roundoff error in doubles of this size
        // becomes significant
        if x.is_infinite() {
            if sign {
                return (-M_PI_2, f64::NAN);
            } else {
                return (M_PI_2, 0.0);
            }
        }
        (M_PI_2 - x.cos() / x, x.sin() / x)

    } else if x > 4.0 {
        /* The auxiliary functions are:
        *
        *
        * *si = *si - M_PI_2;
        * c = cos(x);
        * s = sin(x);
        *
        * t = *ci * s - *si * c;
        * a = *ci * c + *si * s;
        *
        * *si = t;
        * *ci = -a;
        */
        let s = x.sin();
        let c = x.cos();
        let z = 1.0 / (x * x);
        let (f, g) = if x < 8.0 {
            (polevl(z, &FN4, 6) / (x * p1evl(z, &FD4, 7)),
            z * polevl(z, &GN4, 7) / p1evl(z, &GD4, 7))
        } else {
            (polevl(z, &FN8, 8) / (x * p1evl(z, &FD8, 8)), 
            z * polevl(z, &GN8, 8) / p1evl(z, &GD8, 9))
        };
        (M_PI_2 - f * c - g * s, f * s - g * c)
    } else {
        let z = x * x;
        let s = x * polevl(z, &SN, 5) / polevl(z, &SD, 5);
        let c = z * polevl(z, &CN, 5) / polevl(z, &CD, 5);

        (s, EULER + x.ln() + c)	/* real part if x < 0 */
    };

    if sign {
        (-si, ci)
    } else {
        (si, ci)
    }

}

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

    #[test]
    fn sici_trivials() {
        assert_eq!(sici(f64::NAN).0.is_nan(), true);
        assert_eq!(sici(f64::NAN).1.is_nan(), true);
        assert_eq!(sici(0.0), (0.0, -f64::INFINITY));
        assert_eq!(sici(f64::INFINITY), (M_PI_2, 0.0));
        assert_eq!(sici(-f64::INFINITY).0, -M_PI_2);
        assert_eq!(sici(-f64::INFINITY).1.is_nan(), true);
    }

    #[test]
    fn sici_large() {
        assert_eq!(sici(4.0 + 1e-5), (1.7582012469370094, -0.14098333198447985));
        assert_eq!(sici(4.5), (1.654140414379244, -0.1934911221017388));
        assert_eq!(sici(10.0), (1.658347594218874, -0.04545643300445537));
        assert_eq!(sici(100.0), (1.5622254668890563, -0.005148825142610493));
        assert_eq!(sici(2e5), (1.5707913395764488, -3.57284412145905e-07));
        assert_eq!(sici(1e9), (1.5707963259570095, 5.458434486108123e-10));

        assert_eq!(sici(-4.0 - 1e-5), (-1.7582012469370094, -0.14098333198447985));
        assert_eq!(sici(-4.5), (-1.654140414379244, -0.1934911221017388));
        assert_eq!(sici(-10.0), (-1.658347594218874, -0.04545643300445537));
        assert_eq!(sici(-100.0), (-1.5622254668890563, -0.005148825142610493));
        assert_eq!(sici(-2e5), (-1.5707913395764488, -3.57284412145905e-07));
        assert_eq!(sici(-1e9), (-1.5707963259570095, 5.458434486108123e-10));
    }

    #[test]
    fn sici_small() {
        assert_eq!(sici(1e-20), (1e-20, -45.474486194979384));
        assert_eq!(sici(1e-10), (1e-10, -22.448635265038924));
        assert_eq!(sici(1e-5), (9.999999999944446e-06, -10.935709800093695));
        assert_eq!(sici(0.1), (0.09994446110827694, -1.7278683866572966));
        assert_eq!(sici(1.0), (0.9460830703671831, 0.33740392290096816));
        assert_eq!(sici(2.0), (1.605412976802695, 0.422980828774865));
        assert_eq!(sici(3.0), (1.848652527999468, 0.11962978600800067));
        assert_eq!(sici(4.0), (1.758203138949053, -0.1409816978869305));

        assert_eq!(sici(-1e-20), (-1e-20, -45.474486194979384));
        assert_eq!(sici(-1e-10), (-1e-10, -22.448635265038924));
        assert_eq!(sici(-1e-5), (-9.999999999944446e-06, -10.935709800093695));
        assert_eq!(sici(-0.1), (-0.09994446110827694, -1.7278683866572966));
        assert_eq!(sici(-1.0), (-0.9460830703671831, 0.33740392290096816));
        assert_eq!(sici(-2.0), (-1.605412976802695, 0.422980828774865));
        assert_eq!(sici(-3.0), (-1.848652527999468, 0.11962978600800067));
        assert_eq!(sici(-4.0), (-1.758203138949053, -0.1409816978869305));
    }

    #[test]
    fn sici_asy() {
        // Note: The Ci values are imprecise since input precision is reduced at this scale
        assert_eq!(sici(1e9 + 1e-5), (1.5707963259570148, 5.458518396719191e-10));
        assert_eq!(sici(1.1e9), (1.5707963275339918, 5.293245707467659e-10));
        assert_eq!(sici(1e10), (1.5707963267075846, -4.875060250875107e-11));
        assert_eq!(sici(1e15), (1.570796326794897, 8.582727931702359e-16));

        assert_eq!(sici(-1e9 - 1e-5), (-1.5707963259570148, 5.458518396719191e-10));
        assert_eq!(sici(-1.1e9), (-1.5707963275339918, 5.293245707467659e-10));
        assert_eq!(sici(-1e10), (-1.5707963267075846, -4.875060250875107e-11));
        assert_eq!(sici(-1e15), (-1.570796326794897, 8.582727931702359e-16));
    }
}