r2rs-nmath 0.1.1

// Translation of nmath's bessel_j
//
// 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_y, bessel_y_ex},
        betagam::gamma_cody,
    },
    traits::TrigPI,
};

pub fn bessel_j<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 bj: *mut f64 = 0 as *mut f64;
    // let mut vmax: *const libc::c_void = 0 as *const libc::c_void;

    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_j(x, -alpha).unwrap() * alpha.cos_pi()
        }) + (if alpha == na {
            0.0
        } else {
            bessel_y(x, -alpha).unwrap() * alpha.sin_pi()
        }))
        .into();
    } else if alpha > 1e7 {
        warn!(
            "besselJ(x, nu): nu={} too large for bessel_j() algorithm",
            alpha
        ); /* nb-1 <= alpha < nb */
        return f64::nan().into();
    } // ==> alpha' in [0, 1)
    nb = 1 + na as usize;
    alpha -= nb as f64 - 1.0;
    let mut bj = vec![0.0; nb];
    J_bessel(&mut x, &mut alpha, &mut nb, &mut bj, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc < 0 {
            warn!(
                "bessel_j({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            warn!(
                "bessel_j({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = bj[nb - 1];
    // vmaxset(vmax);
    x.into()
}

/// Modified version of `bessel_j`, accepting a work array
/// instead of allocating one.
pub fn bessel_j_ex(mut x: f64, mut alpha: f64, bj: &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_j");
    }
    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_j_ex(x, -alpha, bj)) * alpha.cos_pi()
        }) + (if alpha == na {
            0.0
        } else {
            (bessel_y_ex(x, -alpha, bj)) * alpha.sin_pi()
        });
    } else if alpha > 1e7 {
        warn!(
            "besselJ(x, nu): nu={} too large for bessel_j() algorithm",
            alpha
        ); /* nb-1 <= alpha < nb */
        return f64::nan();
    } // ==> alpha' in [0, 1)
    nb = 1 + na as usize;
    alpha -= nb as f64 - 1.0;
    J_bessel(&mut x, &mut alpha, &mut nb, bj, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc < 0 {
            warn!(
                "bessel_j({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            warn!(
                "bessel_j({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = bj[nb - 1];
    x
}

/// From <http://www.netlib.org/specfun/rjbesl> Fortran translated by f2c,...
/// ------------------------------=#---- Martin Maechler, ETH Zurich
/// Additional code for nu == alpha < 0  MM
///
/// Calculates Bessel functions J_{n+alpha} (x)
/// for non-negative argument x, and non-negative order n+alpha, n = 0,1,..,nb-1.
///
/// Explanation of variables in the calling sequence.
///
/// X     - Non-negative argument for which J's are to be calculated.
/// ALPHA - Fractional part of order for which
/// J's are to be calculated.  0 <= ALPHA < 1.
/// NB    - Number of functions to be calculated, NB >= 1.
/// The first function calculated is of order ALPHA, and the
/// last is of order (NB - 1 + ALPHA).
/// B     - Output vector of length NB.  If RJBESL
/// terminates normally (NCALC=NB), the vector B contains the
/// functions J/ALPHA/(X) through J/NB-1+ALPHA/(X).
/// NCALC - Output variable indicating possible errors.
/// Before using the vector B, the user should check that
/// NCALC=NB, i.e., all orders have been calculated to
/// the desired accuracy. See the following
///
/// ****************************************************************
///
/// Error return codes
///
/// In case of an error,  NCALC != NB, and not all J's are
/// calculated to the desired accuracy.
///
/// NCALC < 0: An argument is out of range. For example,
/// NBES <= 0, ALPHA < 0 or > 1, or X is too large.
/// In this case, b\[1\] is set to zero, the remainder of the
/// B-vector is not calculated, and NCALC is set to
/// MIN(NB,0)-1 so that NCALC != NB.
///
/// NB > NCALC > 0: Not all requested function values could
/// be calculated accurately.  This usually occurs because NB is
/// much larger than ABS(X).  In this case, b\[N\] is calculated
/// to the desired accuracy for N <= NCALC, but precision
/// is lost for NCALC < N <= NB.  If b\[N\] does not vanish
/// for N > NCALC (because it is too small to be represented),
/// and b\[N\]/b\[NCALC\] = 10^(-K), then only the first NSIG - K
/// significant figures of b\[N\] can be trusted.
fn J_bessel(x: &mut f64, alpha: &mut f64, nb: &mut usize, b: &mut [f64], ncalc: &mut i32) {
    let mut current_block: u64;
    /* ---------------------------------------------------------------------
     Mathematical constants

      _2TPI   = 2 / PI
      TWOPI1 = first few significant digits of 2 * PI
      TWOPI2 = (2*PI - TWOPI1) to working precision, i.e.,
           TWOPI1 + TWOPI2 = 2 * PI to extra precision.
    --------------------------------------------------------------------- */
    static pi2: f64 = std::f64::consts::FRAC_2_PI;
    static twopi1: f64 = 6.28125;
    static twopi2: f64 = 0.001935307179586476925286767;
    /*---------------------------------------------------------------------
     *  Factorial(N)
     *--------------------------------------------------------------------- */
    static fact: [f64; 25] = [
        1.0,
        1.0,
        2.0,
        6.0,
        24.0,
        120.0,
        720.0,
        5040.0,
        40320.0,
        362880.0,
        3628800.0,
        39916800.0,
        479001600.0,
        6227020800.0,
        87178291200.0,
        1.307674368e12,
        2.0922789888e13,
        3.55687428096e14,
        6.402373705728e15,
        1.21645100408832e17,
        2.43290200817664e18,
        5.109094217170944e19,
        1.12400072777760768e21,
        2.585201673888497664e22,
        6.2044840173323943936e23,
    ];
    /* Local variables */
    let mut nend = 0;
    let mut intx = 0;
    let mut nbmx = 0;
    let mut i = 0;
    let mut j = 0;
    let mut k = 0;
    let mut l = 0;
    let mut m = 0;
    let mut n = 0;
    let mut nstart = 0;
    let mut nu = 0.0;
    let mut twonu = 0.0;
    let mut capp = 0.0;
    let mut capq = 0.0;
    let mut pold = 0.0;
    let mut vcos = 0.0;
    let mut test = 0.0;
    let mut vsin = 0.0;
    let mut p = 0.0;
    let mut s = 0.0;
    let mut t = 0.0;
    let mut z = 0.0;
    let mut alpem = 0.0;
    let mut halfx = 0.0;
    let mut aa = 0.0;
    let mut bb = 0.0;
    let mut cc = 0.0;
    let mut psave = 0.0;
    let mut plast = 0.0;
    let mut tover = 0.0;
    let mut t1 = 0.0;
    let mut alp2em = 0.0;
    let mut em = 0.0;
    let mut en = 0.0;
    let mut xc = 0.0;
    let mut xk = 0.0;
    let mut xm = 0.0;
    let mut psavel = 0.0;
    let mut gnu = 0.0;
    let mut xin = 0.0;
    let mut sum = 0.0;
    /* Parameter adjustment */
    // b = b.offset(-1);
    nu = *alpha;
    twonu = nu + nu;
    /*-------------------------------------------------------------------
    Check for out of range arguments.
    -------------------------------------------------------------------*/
    if *nb > 0 && *x >= 0.0 && (0.0..1.0).contains(&nu) {
        *ncalc = *nb as i32;
        if *x > 1e5 {
            warn!("value out of range in J_bessel");
            /* indeed, the limit is 0,
             * but the cutoff happens too early */
            i = 1; /*was ML_POSINF (really nonsense) */
            while i <= *nb {
                b[i - 1] = 0.0;
                i += 1
            }
            return;
        }
        intx = *x as usize;
        /* Initialize result array to zero. */
        i = 1;
        while i <= *nb {
            b[i - 1] = 0.0;
            i += 1
        }
        /*===================================================================
        Branch into  3 cases :
        1) use 2-term ascending series for small X
        2) use asymptotic form for large X when NB is not too large
        3) use recursion otherwise
        ===================================================================*/
        if *x < 1e-4 {
            /* ---------------------------------------------------------------
            Two-term ascending series for small X.
            --------------------------------------------------------------- */
            alpem = 1.0 + nu;
            halfx = if *x > 8.9e-308 { (0.5) * *x } else { 0.0 };
            aa = if nu != 0.0 {
                (halfx.powf(nu)) / (nu * gamma_cody(nu).unwrap())
            } else {
                1.0
            };
            bb = if *x + 1.0 > 1.0 {
                (-halfx) * halfx
            } else {
                0.0
            };
            b[0] = aa + aa * bb / alpem;
            if *x != 0.0 && b[0] == 0.0 {
                *ncalc = 0
            }
            if *nb != 1 {
                if *x <= 0.0 {
                    n = 2;
                    while n <= *nb {
                        b[n - 1] = 0.0;
                        n += 1
                    }
                } else {
                    /* ----------------------------------------------
                    Calculate higher order functions.
                    ---------------------------------------------- */
                    if bb == 0.0 {
                        tover = (8.9e-308 + 8.9e-308) / *x
                    } else {
                        tover = 8.9e-308 / bb
                    }
                    cc = halfx;
                    n = 2;
                    while n <= *nb {
                        aa /= alpem;
                        alpem += 1.0;
                        aa *= cc;
                        if aa <= tover * alpem {
                            aa = 0.0
                        }
                        b[n - 1] = aa + aa * bb / alpem;
                        if b[n - 1] == 0.0 && *ncalc > n as i32 {
                            *ncalc = n as i32 - 1
                        }
                        n += 1
                    }
                }
            }
        } else if *x > 25.0 && *nb <= intx + 1 {
            /* ------------------------------------------------------------
            Asymptotic series for X > 25 (and not too large nb)
            ------------------------------------------------------------ */
            xc = (pi2 / *x).sqrt();
            xin = 1.0 / (64.0 * *x * *x);
            if *x >= 130.0 {
                m = 4
            } else if *x >= 35.0 {
                m = 8
            } else {
                m = 11
            }
            xm = 4.0 * m as f64;
            /* ------------------------------------------------
            Argument reduction for SIN and COS routines.
            ------------------------------------------------ */
            t = (*x / (twopi1 + twopi2) + 0.5).trunc();
            z = *x - t * twopi1 - t * twopi2 - (nu + 0.5) / pi2;
            vsin = z.sin();
            vcos = z.cos();
            gnu = twonu;
            i = 1;
            while i <= 2 {
                s = (xm - 1.0 - gnu) * (xm - 1.0 + gnu) * xin * 0.5;
                t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
                t1 = (gnu - (xm + 1.0)) * (gnu + (xm + 1.0));
                k = m + m;
                capp = s * t / fact[k];
                capq = s * t1 / fact[k + 1];
                xk = xm;
                while k >= 4 {
                    /* k + 2(j-2) == 2m */
                    xk -= 4.0;
                    s = (xk - 1.0 - gnu) * (xk - 1.0 + gnu);
                    t1 = t;
                    t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
                    capp = (capp + 1.0 / fact[k - 2]) * s * t * xin;
                    capq = (capq + 1.0 / fact[k - 1]) * s * t1 * xin;
                    k -= 2
                }
                capp += 1.0;
                capq = (capq + 1.0) * (gnu * gnu - 1.0) * (0.125 / *x);
                b[i - 1] = xc * (capp * vcos - capq * vsin);
                if *nb == 1 {
                    return;
                }
                /* vsin <--> vcos */
                t = vsin;
                vsin = -vcos;
                vcos = t;
                gnu += 2.0;
                i += 1
            }
            /* -----------------------------------------------
            If  NB > 2, compute J(X,ORDER+I) for I = 2, NB-1
            ----------------------------------------------- */
            if *nb > 2 {
                gnu = twonu + 2.0;
                j = 3;
                while j <= *nb {
                    b[j - 1] = gnu * b[(j - 1) - 1] / *x - b[(j - 2) - 1];
                    j += 1;
                    gnu += 2.0
                }
            }
        } else {
            /* rtnsig_BESS <= x && ( x <= 25 || intx+1 < *nb ) :
            --------------------------------------------------------
            Use recurrence to generate results.
            First initialize the calculation of P*S.
            -------------------------------------------------------- */
            nbmx = *nb as i32 - intx as i32;
            n = intx + 1;
            en = (n + n) as f64 + twonu;
            plast = 1.0;
            p = en / *x;
            /* ---------------------------------------------------
            Calculate general significance test.
            --------------------------------------------------- */
            test = 1e16 + 1e16;
            if nbmx >= 3 {
                /* ------------------------------------------------------------
                Calculate P*S until N = NB-1.0  Check for possible overflow.
                ---------------------------------------------------------- */
                tover = 1e308 / 1e16;
                nstart = intx + 2;
                nend = *nb - 1;
                en = (nstart + nstart) as f64 - 2.0 + twonu;
                k = nstart;
                's_555: loop {
                    if k > nend {
                        current_block = 1915186496383530739;
                        break;
                    }
                    n = k;
                    en += 2.0;
                    pold = plast;
                    plast = p;
                    p = en * plast / *x - pold;
                    if p > tover {
                        /* -------------------------------------------
                        To avoid overflow, divide P*S by TOVER.
                        Calculate P*S until ABS(P) > 1.
                        -------------------------------------------*/
                        tover = 1e308;
                        p /= tover;
                        plast /= tover;
                        psave = p;
                        psavel = plast;
                        nstart = n + 1;
                        loop {
                            n += 1;
                            en += 2.0;
                            pold = plast;
                            plast = p;
                            p = en * plast / *x - pold;
                            if !(p <= 1.0) {
                                break;
                            }
                        }
                        bb = en / *x;
                        /* -----------------------------------------------
                        Calculate backward test and find NCALC,
                        the highest N such that the test is passed.
                        ----------------------------------------------- */
                        test = pold * plast * (0.5 - 0.5 / (bb * bb));
                        test /= 1e16;
                        p = plast * tover;
                        n -= 1;
                        en -= 2.0;
                        nend = if *nb <= n { *nb } else { n };
                        l = nstart;
                        while l <= nend {
                            pold = psavel;
                            psavel = psave;
                            psave = en * psavel / *x - pold;
                            if psave * psavel > test {
                                *ncalc = l as i32 - 1;
                                current_block = 15945884708764381230;
                                break 's_555;
                            } else {
                                l += 1
                            }
                        }
                        *ncalc = nend as i32;
                        current_block = 15945884708764381230;
                        break;
                    } else {
                        k += 1
                    }
                }
                match current_block {
                    15945884708764381230 => {}
                    _ => {
                        n = nend;
                        en = (n + n) as f64 + twonu;
                        /* -----------------------------------------------------
                        Calculate special significance test for NBMX > 2.
                        -----------------------------------------------------*/
                        test = test.max((plast * 1e16).sqrt() * (p + p).sqrt());
                        current_block = 10601179871800211547;
                    }
                }
            } else {
                current_block = 10601179871800211547;
            }
            if current_block == 10601179871800211547 {
                loop
                /* ------------------------------------------------
                Calculate P*S until significance test passes. */
                {
                    n += 1;
                    en += 2.0;
                    pold = plast;
                    plast = p;
                    p = en * plast / *x - pold;
                    if !(p < test) {
                        break;
                    }
                }
            }
            /*---------------------------------------------------------------
            Initialize the backward recursion and the normalization sum.
            --------------------------------------------------------------- */
            n += 1; /* = 2 n - 4 (n/2)
                    = 0 for even, 2 for odd n */
            en += 2.0;
            bb = 0.0;
            aa = 1.0 / p;
            m = n / 2;
            em = m as f64;
            m = (n << 1) - (m << 2);
            if m == 0 {
                sum = 0.0
            } else {
                alpem = em - 1.0 + nu;
                alp2em = em + em + nu;
                sum = aa * alpem * alp2em / em
            }
            nend = n - *nb;
            /* if (nend > 0) */
            /* --------------------------------------------------------
            Recur backward via difference equation, calculating
            (but not storing) b[N], until N = NB.
            -------------------------------------------------------- */
            l = 1; /* m = 2 - m failed on gcc4-20041019 */
            while l <= nend {
                n -= 1;
                en -= 2.0;
                cc = bb;
                bb = aa;
                aa = en * bb / *x - cc;
                m = if m != 0 { 0 } else { 2 };
                if m != 0 {
                    em -= 1.0;
                    alp2em = em + em + nu;
                    if n == 1 {
                        break;
                    }
                    alpem = em - 1.0 + nu;
                    if alpem == 0.0 {
                        alpem = 1.0
                    }
                    sum = (sum + aa * alp2em) * alpem / em
                }
                l += 1
            }
            /*--------------------------------------------------
            Store b[NB].
            --------------------------------------------------*/
            b[n - 1] = aa;
            if *nb <= 1 {
                if nu + 1.0 == 1.0 {
                    alp2em = 1.0
                } else {
                    alp2em = nu
                }
                sum += b[0] * alp2em;
                current_block = 7905447023580051273;
            } else {
                /*-- nb >= 2 : ---------------------------
                Calculate and store b[NB-1].
                ----------------------------------------*/
                n -= 1; /* m = 2 - m failed on gcc4-20041019 */
                en -= 2.0;
                b[n - 1] = en * aa / *x - bb;
                if n == 1 {
                    current_block = 16540037797593851995;
                } else {
                    m = if m != 0 { 0 } else { 2 };
                    if m != 0 {
                        em -= 1.0;
                        alp2em = em + em + nu;
                        alpem = em - 1.0 + nu;
                        if alpem == 0.0 {
                            alpem = 1.0
                        }
                        sum = (sum + b[n - 1] * alp2em) * alpem / em
                    }
                    current_block = 5872168878400681860;
                }
            }
            if current_block == 5872168878400681860 {
                /* if (n - 2 != 0) */
                /* --------------------------------------------------------
                Calculate via difference equation and store b[N],
                until N = 2.
                -------------------------------------------------------- */
                n = n - 1; /* m = 2 - m failed on gcc4-20041019 */
                while n >= 2 {
                    en -= 2.0;
                    b[n - 1] = en * b[(n + 1) - 1] / *x - b[(n + 2) - 1];
                    m = if m != 0 { 0 } else { 2 };
                    if m != 0 {
                        em -= 1.0;
                        alp2em = em + em + nu;
                        alpem = em - 1.0 + nu;
                        if alpem == 0.0 {
                            alpem = 1.0
                        }
                        sum = (sum + b[n - 1] * alp2em) * alpem / em
                    }
                    n -= 1
                }
                /* ---------------------------------------
                Calculate b[1].
                -----------------------------------------*/
                b[0] = 2.0 * (nu + 1.0) * b[1] / *x - b[2];
                current_block = 16540037797593851995;
            }
            if current_block == 16540037797593851995 {
                em -= 1.0;
                alp2em = em + em + nu;
                if alp2em == 0.0 {
                    alp2em = 1.0
                }
                sum += b[0] * alp2em
            }
            /* ---------------------------------------------------
            Normalize.  Divide all b[N] by sum.
            ---------------------------------------------------*/
            /*     if (nu + 1.0 != 1.0) poor test */
            if nu.abs() > 1e-15 {
                sum *= gamma_cody(nu).unwrap() * (0.5 * *x).powf(-nu)
            }
            aa = 8.9e-308;
            if sum > 1.0 {
                aa *= sum
            }
            n = 1;
            while n <= *nb {
                if (b[n - 1]).abs() < aa {
                    b[n - 1] = 0.0
                } else {
                    b[n - 1] /= sum
                }
                n += 1
            }
        }
    } else {
        /* Error return -- X, NB, or ALPHA is out of range : */
        b[0] = 0.0;
        *ncalc = (if *nb == 0 { *nb } else { 0 }) as i32 - 1
    };
}