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)]


/* S(x) for small x */
const SN: [f64; 6] = [
    -2.99181919401019853726E3,
    7.08840045257738576863E5,
    -6.29741486205862506537E7,
    2.54890880573376359104E9,
    -4.42979518059697779103E10,
    3.18016297876567817986E11,
];

const SD: [f64; 6] = [
    /* 1.00000000000000000000E0, */
    2.81376268889994315696E2,
    4.55847810806532581675E4,
    5.17343888770096400730E6,
    4.19320245898111231129E8,
    2.24411795645340920940E10,
    6.07366389490084639049E11,
];

/* C(x) for small x */
const CN: [f64; 6] = [
    -4.98843114573573548651E-8,
    9.50428062829859605134E-6,
    -6.45191435683965050962E-4,
    1.88843319396703850064E-2,
    -2.05525900955013891793E-1,
    9.99999999999999998822E-1,
];

const CD: [f64; 7] = [
    3.99982968972495980367E-12,
    9.15439215774657478799E-10,
    1.25001862479598821474E-7,
    1.22262789024179030997E-5,
    8.68029542941784300606E-4,
    4.12142090722199792936E-2,
    1.00000000000000000118E0,
];

/* Auxiliary function f(x) */
const FN: [f64; 10] = [
    4.21543555043677546506E-1,
    1.43407919780758885261E-1,
    1.15220955073585758835E-2,
    3.45017939782574027900E-4,
    4.63613749287867322088E-6,
    3.05568983790257605827E-8,
    1.02304514164907233465E-10,
    1.72010743268161828879E-13,
    1.34283276233062758925E-16,
    3.76329711269987889006E-20,
];

const FD: [f64; 10] = [
    /*  1.00000000000000000000E0, */
    7.51586398353378947175E-1,
    1.16888925859191382142E-1,
    6.44051526508858611005E-3,
    1.55934409164153020873E-4,
    1.84627567348930545870E-6,
    1.12699224763999035261E-8,
    3.60140029589371370404E-11,
    5.88754533621578410010E-14,
    4.52001434074129701496E-17,
    1.25443237090011264384E-20,
];

/* Auxiliary function g(x) */
const GN: [f64; 11] = [
    5.04442073643383265887E-1,
    1.97102833525523411709E-1,
    1.87648584092575249293E-2,
    6.84079380915393090172E-4,
    1.15138826111884280931E-5,
    9.82852443688422223854E-8,
    4.45344415861750144738E-10,
    1.08268041139020870318E-12,
    1.37555460633261799868E-15,
    8.36354435630677421531E-19,
    1.86958710162783235106E-22,
];

const GD: [f64; 11] = [
    /*  1.00000000000000000000E0, */
    1.47495759925128324529E0,
    3.37748989120019970451E-1,
    2.53603741420338795122E-2,
    8.14679107184306179049E-4,
    1.27545075667729118702E-5,
    1.04314589657571990585E-7,
    4.60680728146520428211E-10,
    1.10273215066240270757E-12,
    1.38796531259578871258E-15,
    8.39158816283118707363E-19,
    1.86958710162783236342E-22,
];

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

// $$\mathrm{C}(x) = \int_{0}^{x}{\cos\left(\frac{\pi}{2}\,t^2\right)\,dt}$$
// $$\mathrm{S}(x) = \int_{0}^{x}{\sin\left(\frac{\pi}{2}\,t^2\right)\,dt}$$

pub fn fresnl(xxa: f64) -> (f64, f64) // ssa, cca
{
    //! Fresnel integral
    //!
    //! ## DESCRIPTION:
    //!
    //! Evaluates the Fresnel integrals
    #![doc=include_str!("fresnl_c.svg")]
    //!
    #![doc=include_str!("fresnl_s.svg")]
    //!
    //! The integrals are evaluated by a power series for x < 1.
    //! For x >= 1 auxiliary functions f(x) and g(x) are employed
    //! such that
    //!
    //! `C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )`
    //! `S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )`
    //!
    //! ## ACCURACY:
    //!
    //! Relative error.
    //!
    //!<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>S(x)</td>
    //!     <td>0, 10</td>
    //!     <td>10000</td>
    //!     <td>2.0e-15</td>
    //!     <td>3.2e-16</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>C(x)</td>
    //!     <td>0, 10</td>
    //!     <td>10000</td>
    //!     <td>1.8e-15</td>
    //!     <td>3.3e-16</td>
    //! </tr>
    //!</table>

    let (cc, ss) = if xxa.is_infinite() {
        (0.5, 0.5)
    } else {
        let x = xxa.abs();
        let x2 = x * x;
        if x2 < 2.5625 {
            let t = x2 * x2;
            (x * polevl(t, &CN, 5) / polevl(t, &CD, 6),
            x * x2 * polevl(t, &SN, 5) / p1evl(t, &SD, 6))
        } else if x > 36974.0 {
            /*
            * http://functions.wolfram.com/GammaBetaErf/FresnelC/06/02/
            * http://functions.wolfram.com/GammaBetaErf/FresnelS/06/02/
            */
            (0.5 + 1.0 / ( M_PI * x) * (M_PI * x2 / 2.0).sin(),
            0.5 - 1.0 / ( M_PI * x) * (M_PI * x2 / 2.0).cos())
        } else {


            /*             Asymptotic power series auxiliary functions
            *             for large argument
            */
            let mut t = M_PI * x2;
            let u = 1.0 / (t * t);
            t = 1.0 / t;
            let f = 1.0 - u * polevl(u, &FN, 9) / p1evl(u, &FD, 10);
            let g = t * polevl(u, &GN, 10) / p1evl(u, &GD, 11);

            t = M_PI_2 * x2;
            let c = t.cos();
            let s = t.sin();
            t = M_PI * x;

            (0.5 + (f * s - g * c) / t,
            0.5 - (f * c + g * s) / t)
        }
    };

    //done:
    if xxa < 0.0 {
        (-ss, -cc)
    } else {
        (ss, cc)
    }
}

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

    #[test]
    fn fresnl_trivials() {
        assert_eq!(fresnl(f64::INFINITY), (0.5, 0.5));
        assert_eq!(fresnl(-f64::INFINITY), (-0.5, -0.5));
    }

    #[test]
    fn fresnl_small_x() { // x < 1.6007810593582121
        assert_eq!(fresnl(0.0), (0.0, 0.0));

        assert_eq!(fresnl(0.1), (0.0005235895476122108, 0.09999753262708506));
        assert_eq!(fresnl(0.5), (0.06473243285999929, 0.4923442258714464));
        assert_eq!(fresnl(1.0), (0.4382591473903547, 0.779893400376823));
        assert_eq!(fresnl(1.5), (0.697504960082093, 0.44526117603982157));
        assert_eq!(fresnl(1.60078105935821), (0.6382848910003114, 0.36496500031406365));

        assert_eq!(fresnl(-0.1), (-0.0005235895476122108, -0.09999753262708506));
        assert_eq!(fresnl(-0.5), (-0.06473243285999929, -0.4923442258714464));
        assert_eq!(fresnl(-1.0), (-0.4382591473903547, -0.779893400376823));
        assert_eq!(fresnl(-1.5), (-0.697504960082093, -0.44526117603982157));
        assert_eq!(fresnl(-1.60078105935821), (-0.6382848910003114, -0.36496500031406365));
    }

    #[test]
    fn fresnl_large_x() { // x > 36974.0
        assert_eq!(fresnl(36974.0 + 1e-5), (0.49999657449629364, 0.5000078981732503));
        assert_eq!(fresnl(36975.0), (0.4999999999973447, 0.5000086087866445));
        assert_eq!(fresnl(40000.0), (0.4999920422528454, 0.4999999999975119));
        // TODO: Accuracy diverges for large x
        assert_eq!(fresnl(1e10), (0.5000000000311804, 0.5000000000064027));
        assert_eq!(fresnl(1e15), (0.49999999999999983, 0.5000000000000002));

        assert_eq!(fresnl(-36974.0 - 1e-5), (-0.49999657449629364, -0.5000078981732503));
        assert_eq!(fresnl(-36975.0), (-0.4999999999973447, -0.5000086087866445));
        assert_eq!(fresnl(-40000.0), (-0.4999920422528454, -0.4999999999975119));
        // TODO: Accuracy diverges for large x
        assert_eq!(fresnl(-1e10), (-0.5000000000311804, -0.5000000000064027));
        assert_eq!(fresnl(-1e15), (-0.49999999999999983, -0.5000000000000002));
    }

    #[test]
    fn fresnl_nominal() { // x > 1.6007810593582121 AND x <= 36974.0
        assert_eq!(fresnl(1.60078105935822), (0.6382848910003036, 0.3649650003140573));
        assert_eq!(fresnl(2.0), (0.34341567836369824, 0.48825340607534073));
        assert_eq!(fresnl(10.0), (0.46816997858488224, 0.49989869420551575));
        assert_eq!(fresnl(200.0), (0.49840845056938343, 0.4999999873348505));
        assert_eq!(fresnl(1500.0), (0.4997877934092108, 0.4999999999699133));
        assert_eq!(fresnl(5000.0), (0.4999363380227632, 0.4999999999988784));
        assert_eq!(fresnl(10000.0), (0.49996816901138164, 0.49999999999927663));
        assert_eq!(fresnl(20000.0), (0.4999840845056908, 0.4999999999987433));
        assert_eq!(fresnl(30000.0), (0.4999893896704605, 0.5000000000003436));
        assert_eq!(fresnl(36974.0), (0.49999139098052187, 0.4999999999989329));

        assert_eq!(fresnl(-1.60078105935822), (-0.6382848910003036, -0.3649650003140573));
        assert_eq!(fresnl(-2.0), (-0.34341567836369824, -0.48825340607534073));
        assert_eq!(fresnl(-10.0), (-0.46816997858488224, -0.49989869420551575));
        assert_eq!(fresnl(-200.0), (-0.49840845056938343, -0.4999999873348505));
        assert_eq!(fresnl(-1500.0), (-0.4997877934092108, -0.4999999999699133));
        assert_eq!(fresnl(-5000.0), (-0.4999363380227632, -0.4999999999988784));
        assert_eq!(fresnl(-10000.0), (-0.49996816901138164, -0.49999999999927663));
        assert_eq!(fresnl(-20000.0), (-0.4999840845056908, -0.4999999999987433));
        assert_eq!(fresnl(-30000.0), (-0.4999893896704605, -0.5000000000003436));
        assert_eq!(fresnl(-36974.0), (-0.49999139098052187, -0.4999999999989329));
    }

}