spec_math 0.1.6

/*
* Cephes Math Library Release 2.8:  June, 2000
* Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*/

#![allow(clippy::excessive_precision)]

const C1: f64 = 0.35502805388781723926;
const C2: f64 = 0.258819403792806798405;
const SQRT3: f64 = 1.732050807568877293527;
const SQPII: f64 = 5.64189583547756286948E-1;

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

// #ifdef UNK
// #define MAXAIRY 25.77
// #endif

const MAXAIRY: f64 = 103.892;


const AN: [f64; 8] = [
    3.46538101525629032477E-1,
    1.20075952739645805542E1,
    7.62796053615234516538E1,
    1.68089224934630576269E2,
    1.59756391350164413639E2,
    7.05360906840444183113E1,
    1.40264691163389668864E1,
    9.99999999999999995305E-1,
];

const AD: [f64; 8] = [
    5.67594532638770212846E-1,
    1.47562562584847203173E1,
    8.45138970141474626562E1,
    1.77318088145400459522E2,
    1.64234692871529701831E2,
    7.14778400825575695274E1,
    1.40959135607834029598E1,
    1.00000000000000000470E0,
];

const APN: [f64; 8] = [
    6.13759184814035759225E-1,
    1.47454670787755323881E1,
    8.20584123476060982430E1,
    1.71184781360976385540E2,
    1.59317847137141783523E2,
    6.99778599330103016170E1,
    1.39470856980481566958E1,
    1.00000000000000000550E0,
];

const APD: [f64; 8] = [
    3.34203677749736953049E-1,
    1.11810297306158156705E1,
    7.11727352147859965283E1,
    1.58778084372838313640E2,
    1.53206427475809220834E2,
    6.86752304592780337944E1,
    1.38498634758259442477E1,
    9.99999999999999994502E-1,
];

const BN16: [f64; 5] = [
    -2.53240795869364152689E-1,
    5.75285167332467384228E-1,
    -3.29907036873225371650E-1,
    6.44404068948199951727E-2,
    -3.82519546641336734394E-3,
];

const BD16: [f64; 5] = [
    /* 1.00000000000000000000E0, */
    -7.15685095054035237902E0,
    1.06039580715664694291E1,
    -5.23246636471251500874E0,
    9.57395864378383833152E-1,
    -5.50828147163549611107E-2,
];

const BPPN: [f64; 5] = [
    4.65461162774651610328E-1,
    -1.08992173800493920734E0,
    6.38800117371827987759E-1,
    -1.26844349553102907034E-1,
    7.62487844342109852105E-3,
];

const BPPD: [f64; 5] = [
    /* 1.00000000000000000000E0, */
    -8.70622787633159124240E0,
    1.38993162704553213172E1,
    -7.14116144616431159572E0,
    1.34008595960680518666E0,
    -7.84273211323341930448E-2,
];

const AFN: [f64; 9] = [
    -1.31696323418331795333E-1,
    -6.26456544431912369773E-1,
    -6.93158036036933542233E-1,
    -2.79779981545119124951E-1,
    -4.91900132609500318020E-2,
    -4.06265923594885404393E-3,
    -1.59276496239262096340E-4,
    -2.77649108155232920844E-6,
    -1.67787698489114633780E-8,
];

const AFD: [f64; 9] = [
    /* 1.00000000000000000000E0, */
    1.33560420706553243746E1,
    3.26825032795224613948E1,
    2.67367040941499554804E1,
    9.18707402907259625840E0,
    1.47529146771666414581E0,
    1.15687173795188044134E-1,
    4.40291641615211203805E-3,
    7.54720348287414296618E-5,
    4.51850092970580378464E-7,
];

const AGN: [f64; 11] = [
    1.97339932091685679179E-2,
    3.91103029615688277255E-1,
    1.06579897599595591108E0,
    9.39169229816650230044E-1,
    3.51465656105547619242E-1,
    6.33888919628925490927E-2,
    5.85804113048388458567E-3,
    2.82851600836737019778E-4,
    6.98793669997260967291E-6,
    8.11789239554389293311E-8,
    3.41551784765923618484E-10,
];

const AGD: [f64; 10] = [
    /*  1.00000000000000000000E0, */
    9.30892908077441974853E0,
    1.98352928718312140417E1,
    1.55646628932864612953E1,
    5.47686069422975497931E0,
    9.54293611618961883998E-1,
    8.64580826352392193095E-2,
    4.12656523824222607191E-3,
    1.01259085116509135510E-4,
    1.17166733214413521882E-6,
    4.91834570062930015649E-9,
];

const APFN: [f64; 9] = [
    1.85365624022535566142E-1,
    8.86712188052584095637E-1,
    9.87391981747398547272E-1,
    4.01241082318003734092E-1,
    7.10304926289631174579E-2,
    5.90618657995661810071E-3,
    2.33051409401776799569E-4,
    4.08718778289035454598E-6,
    2.48379932900442457853E-8,
];

const APFD: [f64; 9] = [
    /*  1.00000000000000000000E0, */
    1.47345854687502542552E1,
    3.75423933435489594466E1,
    3.14657751203046424330E1,
    1.09969125207298778536E1,
    1.78885054766999417817E0,
    1.41733275753662636873E-1,
    5.44066067017226003627E-3,
    9.39421290654511171663E-5,
    5.65978713036027009243E-7,
];

const APGN: [f64; 11] = [
    -3.55615429033082288335E-2,
    -6.37311518129435504426E-1,
    -1.70856738884312371053E0,
    -1.50221872117316635393E0,
    -5.63606665822102676611E-1,
    -1.02101031120216891789E-1,
    -9.48396695961445269093E-3,
    -4.60325307486780994357E-4,
    -1.14300836484517375919E-5,
    -1.33415518685547420648E-7,
    -5.63803833958893494476E-10,
];

const APGD: [f64; 10] = [
    /*  1.00000000000000000000E0, */
    9.85865801696130355144E0,
    2.16401867356585941885E1,
    1.73130776389749389525E1,
    6.17872175280828766327E0,
    1.08848694396321495475E0,
    9.95005543440888479402E-2,
    4.78468199683886610842E-3,
    1.18159633322838625562E-4,
    1.37480673554219441465E-6,
    5.79912514929147598821E-9,
];

pub fn airy(x: f64) -> (f64, f64, f64, f64) // ai, aip, bi, bip
{
    //! Airy Functions
    //!
    //! ## DESCRIPTION:
    //!
    //! Solution of the differential equation
    //! 
    //! y"(x) = xy
    //! 
    //! The function returns the two independent solutions Ai, Bi
    //! and their first derivatives Ai'(x), Bi'(x).
    //!
    //! Evaluation is by power series summation for small x,
    //! by rational minimax approximations for large x.
    //!
    //! ## ACCURACY:
    //!
    //! Error criterion is absolute when function <= 1, relative
    //! when function > 1, except * denotes relative error criterion.
    //! For large negative x, the absolute error increases as x^1.5.
    //! For large positive x, the relative error increases as x^1.5.
    //!
    //!<table>
    //! <tr>
    //!     <th>Arithmetic</th>
    //!     <th>Domain</th>
    //!     <th>Function</th>
    //!     <th># Trials</th>
    //!     <th>Peak</th>
    //!     <th>RMS</th>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>-10, 0</td>
    //!     <td>Ai</td>
    //!     <td>10000</td>
    //!     <td>1.6e-15</td>
    //!     <td>2.7e-16</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>0, 10</td>
    //!     <td>Ai</td>
    //!     <td>10000</td>
    //!     <td>2.3e-14*</td>
    //!     <td>1.8e-15*</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>-10, 0</td>
    //!     <td>Ai'</td>
    //!     <td>10000</td>
    //!     <td>4.6e-15</td>
    //!     <td>7.6e-16</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>0, 10</td>
    //!     <td>Ai'</td>
    //!     <td>10000</td>
    //!     <td>1.8e-14*</td>
    //!     <td>1.5e-15*</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>-10, 10</td>
    //!     <td>Bi</td>
    //!     <td>30000</td>
    //!     <td>4.2e-15</td>
    //!     <td>5.3e-16*</td>
    //! </tr>
    //! <tr>
    //!     <td>IEEE</td>
    //!     <td>-10, 10</td>
    //!     <td>Bi'</td>
    //!     <td>30000</td>
    //!     <td>4.9e-15</td>
    //!     <td>7.3e-16*</td>
    //! </tr>
    //!</table>
    if x > MAXAIRY {
        return (0.0, 0.0, f64::INFINITY, f64::INFINITY);
    }

    if x < -2.09 {
        let t = (-x).sqrt();
        let zeta = -2.0 * x * t / 3.0;
        let t = t.sqrt();
        let k = SQPII / t;
        let z = 1.0 / zeta;
        let zz = z * z;
        let uf = 1.0 + zz * polevl(zz, &AFN, 8) / p1evl(zz, &AFD, 9);
        let ug = z * polevl(zz, &AGN, 10) / p1evl(zz, &AGD, 10);
        let theta = zeta + 0.25 * M_PI;
        let f = theta.sin();
        let g = theta.cos();
        let ai = k * (f * uf - g * ug);
        let bi = k * (g * uf + f * ug);
        let uf = 1.0 + zz * polevl(zz, &APFN, 8) / p1evl(zz, &APFD, 9);
        let ug = z * polevl(zz, &APGN, 10) / p1evl(zz, &APGD, 10);
        let k = SQPII * t;
        let aip = -k * (g * uf + f * ug);
        let bip = k * (f * uf - g * ug);
        return (ai, aip, bi, bip);
    }
    let mut ai = 0.0; // Value doesn't matter, just declaring variable
    let mut aip = 0.0; // Value doesn't matter, just declaring vairable
    let ai_needed = if x >= 2.09 {		/* cbrt(9) */
        let t = x.sqrt();
        let zeta = 2.0 * x * t / 3.0;
        let g = zeta.exp();
        let t = t.sqrt();
        let k = 2.0 * t * g;
        let z = 1.0 / zeta;
        let f = polevl(z, &AN, 7) / polevl(z, &AD, 7);
        ai = SQPII * f / k;
        let k = -0.5 * SQPII * t / g;
        let f = polevl(z, &APN, 7) / polevl(z, &APD, 7);
        aip = f * k;

        if x > 8.3203353 {	/* zeta > 16 */
            let f = z * polevl(z, &BN16, 4) / p1evl(z, &BD16, 5);
            let k = SQPII * g;
            let bi = k * (1.0 + f) / t;
            let f = z * polevl(z, &BPPN, 4) / p1evl(z, &BPPD, 5);
            let bip = k * t * (1.0 + f);
            return (ai, aip, bi, bip);
        }
        false
    } else {
        true
    };

    let mut f = 1.0;
    let mut g = x;
    let mut t = 1.0;
    let mut uf = 1.0;
    let mut ug = x;
    let mut k = 1.0;
    let z = x * x * x;
    while t > MACHEP {
        uf *= z;
        k += 1.0;
        uf /= k;
        ug *= z;
        k += 1.0;
        ug /= k;
        uf /= k;
        f += uf;
        k += 1.0;
        ug /= k;
        g += ug;
        t = (uf / f).abs();
    }
    uf = C1 * f;
    ug = C2 * g;
    if ai_needed {
        ai = uf - ug;
    }
    let bi = SQRT3 * (uf + ug);

    /* the deriviative of ai */
    k = 4.0;
    uf = x * x / 2.0;
    ug = z / 3.0;
    f = uf;
    g = 1.0 + ug;
    uf /= 3.0;
    t = 1.0;

    while t > MACHEP {
        uf *= z;
        ug /= k;
        k += 1.0;
        ug *= z;
        uf /= k;
        f += uf;
        k += 1.0;
        ug /= k;
        uf /= k;
        g += ug;
        k += 1.0;
        t = (ug / g).abs();
    }

    uf = C1 * f;
    ug = C2 * g;
    if ai_needed {
        aip = uf - ug;
    }
    let bip = SQRT3 * (uf + ug);
    
    (ai, aip, bi, bip)
}

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

    fn tuple_is_nan(t: (f64, f64, f64, f64)) -> bool {
        t.0.is_nan() && t.1.is_nan() && t.2.is_nan() && t.3.is_nan()
    }

    #[test]
    fn airy_trivials() {
        assert_eq!(airy(103.892 + 1e-10), (0.0, 0.0, f64::INFINITY, f64::INFINITY));
        assert_eq!(airy(f64::INFINITY), (0.0, 0.0, f64::INFINITY, f64::INFINITY));
        assert_eq!(tuple_is_nan(airy(f64::NAN)), true);
        // TODO: Boundary conditions at -infinity
        //assert_eq!(airy(-f64::INFINITY), ???);
    }

    #[test]
    fn airy_mdeium_x() {
        assert_eq!(airy(-2.09), (0.17005055173203007, 0.654846906732163, -0.43393986253871786, 0.20079740491720918));
        assert_eq!(airy(-2.0), (0.22740742820168564, 0.618259020741691, -0.4123025879563985, 0.27879516692116973));
        assert_eq!(airy(-1.0), (0.5355608832923522, -0.010160567116645175, 0.10399738949694468, 0.5923756264227923));
        assert_eq!(airy(0.0), (0.3550280538878172, -0.2588194037928068, 0.6149266274460007, 0.4482883573538264));
        assert_eq!(airy(1.0), (0.13529241631288147, -0.15914744129679328, 1.2074235949528715, 0.9324359333927756));
        assert_eq!(airy(2.0), (0.03492413042327436, -0.05309038443365388, 3.2980949999782143, 4.10068204993289));
        assert_eq!(airy(2.09 - 1e-10), (0.03042031836790704, -0.0470883854429367, 3.6953287957997896, 4.743632783969338));
    }

    #[test]
    fn airy_small_x() {
        // TODO: Accuracy diverges with large negative x
        //assert_eq!(airy(-1e+20), (0.0, 0.0, 0.0, 0.0));
        //assert_eq!(airy(-1e+10), 23.141782445598864);
        assert_eq!(airy(-100.0), (0.1767533932395512, -0.24229703166070582, 0.02427388768017233, 1.7675948932340444));
        assert_eq!(airy(-10.0), (0.040241238486441955, 0.9962650441327905, -0.3146798296438388, 0.11941411339990535));
        assert_eq!(airy(-3.0), (-0.37881429367765806, 0.314583769216599, -0.19828962637492664, -0.6756112226852584));
        assert_eq!(airy(-2.1), (0.1634845129992927, 0.6583406928143434, -0.4359023482307267, 0.19168563232842983));
        assert_eq!(airy(-2.09 - 1e-20), (0.17005055173203007, 0.654846906732163, -0.43393986253871786, 0.20079740491720918));
    }

    #[test]
    fn airy_large_x() {
        assert_eq!(airy(2.09), (0.03042031836319837, -0.04708838543657822, 3.6953287962741523, 4.743632784741663));
        assert_eq!(airy(2.1), (0.02995260211586652, -0.046455994032674586, 3.7431535649561756, 4.821549925717496));
        assert_eq!(airy(3.0), (0.006591139357460717, -0.011912976705951313, 14.037328963730229, 22.92221496638217));
        assert_eq!(airy(10.0), (1.1047532552898654e-10, -3.520633676738912e-10, 455641153.54822654, 1429236134.48287));
        assert_eq!(airy(100.0), (2.6344821520882847e-291, -2.63514036160451e-290, 6.041223996669972e+288, 6.039712745310374e+289));
        assert_eq!(airy(103.892), (2.240778287387011e-308, -2.2845065026831452e-307, 6.968354736681834e305, 7.100986741522436e306));
    }
}