r2rs-nmath 0.1.1

// Translation of nmath's bessel_y
//
// Copyright (C) 2020 neek-sss
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.

use num_traits::Float;
use strafe_type::{FloatConstraint, Positive64, Real64};

use crate::{
    func::bessel::{bessel_j, bessel_j_ex},
    traits::TrigPI,
};

pub fn bessel_y<P: Into<Positive64>, R: Into<Real64>>(x: P, alpha: R) -> Real64 {
    let mut x = x.into().unwrap();
    let mut alpha = alpha.into().unwrap();

    let mut nb = 0;
    let mut ncalc = 0;
    let mut na = 0.0;
    // let mut by: *mut f64 = 0 as *mut f64;
    // let mut vmax: *const libc::c_void = 0 as *const libc::c_void;
    /* NaNs propagated correctly */

    na = alpha.floor();
    if alpha < 0.0 {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
        return ((if alpha - na == 0.5 {
            0.0
        } else {
            bessel_y(x, -alpha).unwrap() * alpha.cos_pi()
        }) - (if alpha == na {
            0.0
        } else {
            bessel_j(x, -alpha).unwrap() * alpha.sin_pi()
        }))
        .into();
    } else if alpha > 1e7 {
        warn!(
            "besselY(x, nu): nu={} too large for bessel_y() algorithm",
            alpha
        );
        return f64::nan().into();
    }
    nb = 1 + na as usize;
    alpha -= nb as f64 - 1.0;
    // vmax = vmaxget();
    let mut by = vec![0.0; nb];
    // by = R_alloc(
    // nb as u64,
    // ::std::mem::size_of::<f64>() as u64,
    // ) as *mut f64;
    Y_bessel(&mut x, &mut alpha, &mut nb, &mut by, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc == -(1) {
            // vmaxset(vmax);
            return f64::infinity().into();
        } else if ncalc < -1 {
            warn!(
                "bessel_y({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            /* ncalc >= 0 */
            warn!(
                "bessel_y({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = by[nb - 1];
    // vmaxset(vmax);
    x.into()
}

/// Modified version of `bessel_y`, accepting a work array
/// instead of allocating one.
pub fn bessel_y_ex(mut x: f64, mut alpha: f64, by: &mut [f64]) -> f64 {
    let mut nb = 0;
    let mut ncalc = 0;
    let mut na = 0.0;
    /* NaNs propagated correctly */
    if x.is_nan() || alpha.is_nan() {
        return x + alpha;
    }
    if x < 0.0 {
        warn!("value out of range in bessel_y");
    }
    na = alpha.floor();
    if alpha < 0.0 {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
        return (if alpha - na == 0.5 {
            0.0
        } else {
            (bessel_y_ex(x, -alpha, by)) * alpha.cos_pi()
        }) - (if alpha == na {
            0.0
        } else {
            (bessel_j_ex(x, -alpha, by)) * alpha.sin_pi()
        });
    } else if alpha > 1e7 {
        warn!(
            "besselY(x, nu): nu={} too large for bessel_y() algorithm",
            alpha
        );
        return f64::nan();
    }
    nb = 1 + na as usize;
    alpha -= nb as f64 - 1.0;
    Y_bessel(&mut x, &mut alpha, &mut nb, by, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc == -(1) {
            return f64::infinity();
        } else if ncalc < -(1) {
            warn!(
                "bessel_y({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            /* ncalc >= 0 */
            warn!(
                "bessel_y({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = by[nb - 1];
    x
}

/// From <http://www.netlib.org/specfun/rybesl> Fortran translated by f2c,...
/// * ------------------------------=#---- Martin Maechler, ETH Zurich
///
/// ----------------------------------------------------------------------
///
/// This routine calculates Bessel functions Y_(N+ALPHA) (X)
///v for non-negative argument X, and non-negative order N+ALPHA.
///
///
/// Explanation of variables in the calling sequence
///
/// X     - Non-negative argument for which
/// Y's are to be calculated.
/// ALPHA - Fractional part of order for which
/// Y's are to be calculated.  0 <= ALPHA < 1.0.
/// NB    - Number of functions to be calculated, NB > 0.
/// The first function calculated is of order ALPHA, and the
/// last is of order (NB - 1 + ALPHA).
/// BY    - Output vector of length NB. If the
/// routine terminates normally (NCALC=NB), the vector BY
/// contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
/// If (0 < NCALC < NB), BY(I) contains correct function
/// values for I <= NCALC, and contains the ratios
/// Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
/// NCALC - Output variable indicating possible errors.
/// Before using the vector BY, the user should check that
/// NCALC=NB, i.e., all orders have been calculated to
/// the desired accuracy. See error returns below.
///
///
/// *******************************************************************
///
/// Error returns
///
/// In case of an error, NCALC != NB, and not all Y's are
/// calculated to the desired accuracy.
///
/// NCALC < -1:  An argument is out of range. For example,
/// NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >=
/// XMAX.  In this case, BY\[0\] = 0.0, the remainder of the
/// BY-vector is not calculated, and NCALC is set to
/// MIN0(NB,0)-2  so that NCALC != NB.
/// NCALC = -1:  Y(ALPHA,X) >= XINF.  The requested function
/// values are set to 0.0.
/// 1 < NCALC < NB: Not all requested function values could
/// be calculated accurately.  BY(I) contains correct function
/// values for I <= NCALC, and and the remaining NB-NCALC
/// array elements contain 0.0.
fn Y_bessel(x: &mut f64, alpha: &mut f64, nb: &mut usize, by: &mut [f64], ncalc: &mut i32) {
    /* ----------------------------------------------------------------------
     Mathematical constants
       FIVPI = 5*PI
       PIM5 = 5*PI - 15
    ----------------------------------------------------------------------*/
    static fivpi: f64 = 15.707963267948966192;
    static pim5: f64 = 0.70796326794896619231;
    /*----------------------------------------------------------------------
    Coefficients for Chebyshev polynomial expansion of
    1/gamma(1-x), abs(x) <= .5
    ----------------------------------------------------------------------*/
    static ch: [f64; 21] = [
        -6.7735241822398840964e-24,
        -6.1455180116049879894e-23,
        2.9017595056104745456e-21,
        1.3639417919073099464e-19,
        2.3826220476859635824e-18,
        -9.0642907957550702534e-18,
        -1.4943667065169001769e-15,
        -3.3919078305362211264e-14,
        -1.7023776642512729175e-13,
        9.1609750938768647911e-12,
        2.4230957900482704055e-10,
        1.7451364971382984243e-9,
        -3.3126119768180852711e-8,
        -8.6592079961391259661e-7,
        -4.9717367041957398581e-6,
        7.6309597585908126618e-5,
        0.0012719271366545622927,
        0.0017063050710955562222,
        -0.07685284084478667369,
        -0.28387654227602353814,
        0.92187029365045265648,
    ];
    /* Local variables */
    let mut i = 0; /* -Wall */
    let mut k = 0;
    let mut na = 0;
    let mut alfa = 0.0;
    let mut div = 0.0;
    let mut ddiv = 0.0;
    let mut even = 0.0;
    let mut gamma = 0.0;
    let mut term = 0.0;
    let mut cosmu = 0.0;
    let mut sinmu = 0.0;
    let mut b = 0.0;
    let mut c = 0.0;
    let mut d = 0.0;
    let mut e = 0.0;
    let mut f = 0.0;
    let mut g = 0.0;
    let mut h = 0.0;
    let mut p = 0.0;
    let mut q = 0.0;
    let mut r = 0.0;
    let mut s = 0.0;
    let mut d1 = 0.0;
    let mut d2 = 0.0;
    let mut q0 = 0.0;
    let mut pa = 0.0;
    let mut pa1 = 0.0;
    let mut qa = 0.0;
    let mut qa1 = 0.0;
    let mut en = 0.0;
    let mut en1 = 0.0;
    let mut nu = 0.0;
    let mut ex = 0.0;
    let mut ya = 0.0;
    let mut ya1 = 0.0;
    let mut twobyx = 0.0;
    let mut den = 0.0;
    let mut odd = 0.0;
    let mut aye = 0.0;
    let mut dmu = 0.0;
    let mut x2 = 0.0;
    let mut xna = 0.0;
    ya1 = 0.0;
    ya = ya1;
    en1 = ya;
    ex = *x;
    nu = *alpha;
    if *nb > 0 && (0.0..1.0).contains(&nu) {
        if !(2.2250738585072014e-308..=1e8).contains(&ex) {
            /* Warning is not really appropriate, give
             * proper limit:
             * ML_ERROR(ME_RANGE, "Y_bessel"); */
            *ncalc = *nb as i32; /*was ML_POSINF */
            if ex > 1e8 {
                by[0] = 0.0
            } else if ex < 2.2250738585072014e-308 {
                by[0] = f64::NEG_INFINITY
            }
            i = 0;
            while i < *nb {
                by[i] = by[0];
                i += 1
            }
            return;
        }
        xna = (nu + 0.5).trunc();
        na = xna as usize;
        if na == 1 {
            /* was 0 */
            /* <==>  .5 <= *alpha < 1  <==>  -5. <= nu < 0 */
            nu -= xna
        }
        if nu == -0.5 {
            p = strafe_consts::SQRT_2DPI / ex.sqrt();
            ya = p * ex.sin();
            ya1 = -p * ex.cos()
        } else if ex < 3.0 {
            /* -------------------------------------------------------------
            Use Temme's scheme for small X
            ------------------------------------------------------------- */
            b = ex * 0.5;
            d = -b.ln();
            f = nu * d;
            e = b.powf(-nu);
            if nu.abs() < 2.149e-8 {
                c = std::f64::consts::FRAC_1_PI
            } else {
                c = nu / nu.sin_pi()
            }
            /* ------------------------------------------------------------
            Computation of sinh(f)/f
            ------------------------------------------------------------ */
            if f.abs() < 1.0 {
                x2 = f * f;
                en = 19.0;
                s = 1.0;
                i = 1;
                while i <= 9 {
                    s = s * x2 / en / (en - 1.0) + 1.0;
                    en -= 2.0;
                    i += 1
                }
            } else {
                s = (e - 1.0 / e) * 0.5 / f
            }
            /* --------------------------------------------------------
            Computation of 1/gamma(1-a) using Chebyshev polynomials */
            x2 = nu * nu * 8.0;
            aye = ch[0];
            even = 0.0;
            alfa = ch[1];
            odd = 0.0;
            i = 3;
            while i <= 19 {
                even = -(aye + aye + even);
                aye = -even * x2 - aye + ch[i - 1];
                odd = -(alfa + alfa + odd);
                alfa = -odd * x2 - alfa + ch[i];
                i += 2
            }
            even = (even * 0.5 + aye) * x2 - aye + ch[20];
            odd = (odd + alfa) * 2.0;
            gamma = odd * nu + even;
            /* End of computation of 1/gamma(1-a)
            ----------------------------------------------------------- */
            g = e * gamma; /* x > thresh_BESS_Y */
            e = (e + 1.0 / e) * 0.5;
            f = 2.0 * c * (odd * e + even * s * d);
            e = nu * nu;
            p = g * c;
            q = std::f64::consts::FRAC_1_PI / g;
            c = nu * std::f64::consts::FRAC_PI_2;
            if c.abs() < 2.149e-8 {
                r = 1.0
            } else {
                r = (nu / 2.0).sin_pi() / c
            }
            r = std::f64::consts::PI * c * r * r;
            c = 1.0;
            d = -b * b;
            h = 0.0;
            ya = f + r * q;
            ya1 = p;
            en = 1.0;
            while (g / (1.0 + ya.abs())).abs() + (h / (1.0 + ya1.abs())).abs()
                > 2.2204460492503131e-16
            {
                f = (f * en + p + q) / (en * en - e);
                c *= d / en;
                p /= en - nu;
                q /= en + nu;
                g = c * (f + r * q);
                h = c * p - en * g;
                ya += g;
                ya1 += h;
                en += 1.0
            }
            ya = -ya;
            ya1 = -ya1 / b
        } else if ex < 16.0 {
            /* --------------------------------------------------------------
            Use Temme's scheme for moderate X :  3 <= x < 16
            -------------------------------------------------------------- */
            c = (0.5 - nu) * (0.5 + nu);
            b = ex + ex;
            e = ex * std::f64::consts::FRAC_1_PI * nu.cos_pi() / 2.2204460492503131e-16;
            e *= e;
            p = 1.0;
            q = -ex;
            r = 1.0 + ex * ex;
            s = r;
            en = 2.0;
            while r * en * en < e {
                en1 = en + 1.0;
                d = (en - 1.0 + c / en) / s;
                p = (en + en - p * d) / en1;
                q = (-b + q * d) / en1;
                s = p * p + q * q;
                r *= s;
                en = en1
            }
            f = p / s;
            p = f;
            g = -q / s;
            q = g;
            loop {
                en -= 1.0;
                if !(en > 0.0) {
                    break;
                }
                r = en1 * (2.0 - p) - 2.0;
                s = b + en1 * q;
                d = (en - 1.0 + c / en) / (r * r + s * s);
                p = d * r;
                q = d * s;
                e = f + 1.0;
                f = p * e - g * q;
                g = q * e + p * g;
                en1 = en
            }
            f += 1.0;
            d = f * f + g * g;
            pa = f / d;
            qa = -g / d;
            d = nu + 0.5 - p;
            q += ex;
            pa1 = (pa * q - qa * d) / ex;
            qa1 = (qa * q + pa * d) / ex;
            b = ex - std::f64::consts::FRAC_PI_2 * (nu + 0.5);
            c = b.cos();
            s = b.sin();
            d = strafe_consts::SQRT_2DPI / ex.sqrt();
            ya = d * (pa * s + qa * c);
            ya1 = d * (qa1 * s - pa1 * c)
        } else {
            /* ----------------------------------------------------------
            Use Campbell's asymptotic scheme.
            ---------------------------------------------------------- */
            na = 0;
            d1 = (ex / fivpi).trunc();
            i = d1 as usize;
            dmu = ex - 15.0 * d1 - d1 * pim5 - (*alpha + 0.5) * std::f64::consts::FRAC_PI_2;
            if i - ((i / 2) << 1) == 0 {
                cosmu = dmu.cos();
                sinmu = dmu.sin()
            } else {
                cosmu = -dmu.cos();
                sinmu = -dmu.sin()
            }
            ddiv = 8.0 * ex;
            dmu = *alpha;
            den = ex.sqrt();
            k = 1;
            while k <= 2 {
                p = cosmu;
                cosmu = sinmu;
                sinmu = -p;
                d1 = (2.0 * dmu - 1.0) * (2.0 * dmu + 1.0);
                d2 = 0.0;
                div = ddiv;
                p = 0.0;
                q = 0.0;
                q0 = d1 / div;
                term = q0;
                i = 2;
                while i <= 20 {
                    d2 += 8.0;
                    d1 -= d2;
                    div += ddiv;
                    term = -term * d1 / div;
                    p += term;
                    d2 += 8.0;
                    d1 -= d2;
                    div += ddiv;
                    term *= d1 / div;
                    q += term;
                    if term.abs() <= 2.2204460492503131e-16 {
                        break;
                    }
                    i += 1
                }
                p += 1.0;
                q += q0;
                if k == 1 {
                    ya = strafe_consts::SQRT_2DPI * (p * cosmu - q * sinmu) / den
                } else {
                    ya1 = strafe_consts::SQRT_2DPI * (p * cosmu - q * sinmu) / den
                }
                dmu += 1.0;
                k += 1
            }
        }
        if na == 1 {
            h = 2.0 * (nu + 1.0) / ex;
            if h > 1.0 && ya1.abs() > f64::MAX / h {
                h = 0.0;
                ya = 0.0
            }
            h = h * ya1 - ya;
            ya = ya1;
            ya1 = h
        }
        /* ---------------------------------------------------------------
        Now have first one or two Y's
        --------------------------------------------------------------- */
        by[0] = ya;
        *ncalc = 1;
        if *nb > 1 {
            by[1] = ya1;
            if ya1 != 0.0 {
                aye = 1.0 + *alpha;
                twobyx = 2.0 / ex;
                *ncalc = 2;
                i = 2;
                while i < *nb {
                    if twobyx < 1.0 {
                        if (by[i - 1]).abs() * twobyx >= f64::MAX / aye {
                            break;
                        }
                    } else if (by[i - 1]).abs() >= f64::MAX / aye / twobyx {
                        break;
                    }
                    by[i] = twobyx * aye * by[i - 1] - by[i - 2];
                    aye += 1.0;
                    *ncalc += 1;
                    i += 1
                }
            }
        }
        i = *ncalc as usize;
        while i < *nb {
            by[i] = f64::NEG_INFINITY;
            i += 1
        }
    } else {
        by[0] = 0.0;
        *ncalc = (if *nb == 0 { *nb } else { 0 }) as i32 - 1
    };
}