r2rs-nmath 0.1.1

// Translation of nmath's bessel_k
//
// 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};

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

    let mut nb = 0;
    let mut ncalc = 0;
    let mut ize = 0;
    // let mut bk: *mut f64 = 0 as *mut f64;
    // let mut vmax: *const libc::c_void = 0 as *const libc::c_void;

    ize = if expo { 2 } else { 1 };
    if alpha < 0.0 {
        alpha = -alpha
    }
    nb = 1 + alpha.floor() as usize;
    alpha -= nb as f64 - 1.0;
    // vmax = vmaxget();
    let mut bk = vec![0.0; nb];
    // bk = R_alloc(
    // nb as u64,
    // ::std::mem::size_of::<f64>() as u64,
    // ) as *mut f64;
    K_bessel(&mut x, &mut alpha, &mut nb, &mut ize, &mut bk, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc < 0 {
            warn!(
                "bessel_k({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            warn!(
                "bessel_k({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = bk[nb - 1];
    // vmaxset(vmax);
    return x.into();
}

/// Modified version of `bessel_k` that accepts a work array instead of
/// allocating one.
pub fn bessel_k_ex(mut x: f64, mut alpha: f64, expo: f64, bk: &mut [f64]) -> f64 {
    let mut nb = 0;
    let mut ncalc = 0;
    let mut ize = 0;
    /* NaNs propagated correctly */
    if x.is_nan() || alpha.is_nan() {
        return x + alpha;
    } /* nb-1 <= |alpha| < nb */
    if x < 0.0 {
        warn!("value out of range in bessel_k");
    }
    ize = expo as i32;
    if alpha < 0.0 {
        alpha = -alpha
    }
    nb = 1 + alpha.floor() as usize;
    alpha -= nb as f64 - 1.0;
    K_bessel(&mut x, &mut alpha, &mut nb, &mut ize, bk, &mut ncalc);
    if ncalc != nb as i32 {
        /* error input */
        if ncalc < 0 {
            warn!(
                "bessel_k({}): ncalc (={}) != nb (={}); alpha={}. Arg. out of range?",
                x, ncalc, nb, alpha
            );
        } else {
            warn!(
                "bessel_k({},nu={}): precision lost in result",
                x,
                alpha + nb as f64 - 1.0
            );
        }
    }
    x = bk[nb - 1];
    return x;
}

/// From <http://www.netlib.org/specfun/rkbesl> Fortran translated by f2c,...
/// ------------------------------=#---- Martin Maechler, ETH Zurich
///
///-------------------------------------------------------------------
///
/// This routine calculates modified Bessel functions
/// of the third kind, K_(N+ALPHA) (X), for non-negative
/// argument X, and non-negative order N+ALPHA, with or without
/// exponential scaling.
///
/// Explanation of variables in the calling sequence
///
/// X     - Non-negative argument for which
/// K's or exponentially scaled K's (K*EXP(X))
/// are to be calculated. If K's are to be calculated,
/// X must not be greater than XMAX_BESS_K.
/// ALPHA - Fractional part of order for which
/// K's or exponentially scaled K's (K*EXP(X)) 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).
/// IZE   - Type. IZE = 1 if unscaled K's are to be calculated,
/// = 2 if exponentially scaled K's are to be calculated.
/// BK    - Output vector of length NB. If the
/// routine terminates normally (NCALC=NB), the vector BK
/// contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
/// or the corresponding exponentially scaled functions.
/// If (0 < NCALC < NB), BK(I) contains correct function
/// values for I <= NCALC, and contains the ratios
/// K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
/// NCALC - Output variable indicating possible errors.
/// Before using the vector BK, 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 K'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_BESS_K.
/// In this case, the B-vector is not calculated,
/// and NCALC is set to MIN0(NB,0)-2  so that NCALC != NB.
/// NCALC = -1:  Either  K(ALPHA,X) >= XINF  or
/// K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF.  In this case,
/// the B-vector is not calculated. Note that again
/// NCALC != NB.
///
/// 0 < NCALC < NB: Not all requested function values could
/// be calculated accurately.  BK(I) contains correct function
/// values for I <= NCALC, and contains the ratios
/// K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
fn K_bessel(
    x: &mut f64,
    alpha: &mut f64,
    nb: &mut usize,
    ize: &mut i32,
    bk: &mut [f64],
    ncalc: &mut i32,
) {
    /*---------------------------------------------------------------------
    * Mathematical constants
    * A = LOG(2) - Euler's constant
    * D = SQRT(2/PI)
    ---------------------------------------------------------------------*/
    static a: f64 = 0.11593151565841244881;
    /*---------------------------------------------------------------------
    P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
    Coefficients converted from hex to decimal and modified
    by W. J. Cody, 2/26/82 */
    static p: [f64; 8] = [
        0.805629875690432845,
        20.4045500205365151,
        157.705605106676174,
        536.671116469207504,
        900.382759291288778,
        730.923886650660393,
        229.299301509425145,
        0.822467033424113231,
    ];
    static q: [f64; 7] = [
        29.4601986247850434,
        277.577868510221208,
        1206.70325591027438,
        2762.91444159791519,
        3443.74050506564618,
        2210.63190113378647,
        572.267338359892221,
    ];
    /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */
    static r: [f64; 5] = [
        -0.48672575865218401848,
        13.079485869097804016,
        -101.96490580880537526,
        347.65409106507813131,
        3.495898124521934782e-4,
    ];
    static s: [f64; 4] = [
        -25.579105509976461286,
        212.57260432226544008,
        -610.69018684944109624,
        422.69668805777760407,
    ];
    /* T    - Approximation for SINH(Y)/Y */
    static t: [f64; 6] = [
        1.6125990452916363814e-10,
        2.5051878502858255354e-8,
        2.7557319615147964774e-6,
        1.9841269840928373686e-4,
        0.0083333333333334751799,
        0.16666666666666666446,
    ];
    /*---------------------------------------------------------------------*/
    static estm: [f64; 6] = [52.0583, 5.7607, 2.7782, 14.4303, 185.3004, 9.3715];
    static estf: [f64; 7] = [41.8341, 7.1075, 6.4306, 42.511, 1.35633, 84.5096, 20.0];
    /* Local variables */
    let mut iend = 0; /* -Wall */
    let mut i = 0; /* would only have underflow */
    let mut j = 0;
    let mut k = 0;
    let mut m = 0;
    let mut ii = 0;
    let mut mplus1 = 0;
    let mut x2by4 = 0.0;
    let mut twox = 0.0;
    let mut c = 0.0;
    let mut blpha = 0.0;
    let mut ratio = 0.0;
    let mut wminf = 0.0;
    let mut d1 = 0.0;
    let mut d2 = 0.0;
    let mut d3 = 0.0;
    let mut f0 = 0.0;
    let mut f1 = 0.0;
    let mut f2 = 0.0;
    let mut p0 = 0.0;
    let mut q0 = 0.0;
    let mut t1 = 0.0;
    let mut t2 = 0.0;
    let mut twonu = 0.0;
    let mut dm = 0.0;
    let mut ex = 0.0;
    let mut bk1 = 0.0;
    let mut bk2 = 0.0;
    let mut nu = 0.0;
    ii = 0;
    ex = *x;
    nu = *alpha;
    *ncalc = (if *nb == 0 { *nb } else { 0 }) as i32 - 2;
    if *nb > 0 && (0.0 <= nu && nu < 1.0) && (1 <= *ize && *ize <= 2) {
        let current_block_262: u64;
        if ex <= 0.0 || *ize == 1 && ex > 705.342 {
            if ex <= 0.0 {
                if ex < 0.0 {
                    warn!("value out of range in K_bessel");
                }
                i = 0;
                while i < *nb {
                    bk[i] = f64::infinity();
                    i += 1
                }
            } else {
                i = 0;
                while i < *nb {
                    bk[i] = 0.0;
                    i += 1
                }
            }
            *ncalc = *nb as i32;
            return;
        }
        k = 0;
        if nu < 1.49e-154 {
            nu = 0.0
        } else if nu > 0.5 {
            k = 1;
            nu -= 1.0
        }
        twonu = nu + nu;
        iend = *nb + k - 1;
        c = nu * nu;
        d3 = -c;
        if ex <= 1.0 {
            /* ------------------------------------------------------------
            Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
                   Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
            ------------------------------------------------------------ */
            d1 = 0.0;
            d2 = p[0];
            t1 = 1.0;
            t2 = q[0];
            i = 2;
            while i <= 7 {
                d1 = c * d1 + p[i - 1];
                d2 = c * d2 + p[i];
                t1 = c * t1 + q[i - 1];
                t2 = c * t2 + q[i];
                i += 2
            }
            d1 = nu * d1;
            t1 = nu * t1;
            f1 = ex.ln();
            f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
            q0 = (-nu * (a - nu * (p[7] + nu * (d1 - d2) / (t1 - t2)) - f1)).exp();
            f1 = nu * f0;
            p0 = f1.exp();
            /* -----------------------------------------------------------
            Calculation of F0 =
            ----------------------------------------------------------- */
            d1 = r[4];
            t1 = 1.0;
            i = 0;
            while i < 4 {
                d1 = c * d1 + r[i];
                t1 = c * t1 + s[i];
                i += 1
            }
            /* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0)
             *    = f0 * sinh(f1)/f1 */
            if f1.abs() <= 0.5 {
                f1 *= f1;
                d2 = 0.0;
                i = 0;
                while i < 6 {
                    d2 = f1 * d2 + t[i];
                    i += 1
                }
                d2 = f0 + f0 * f1 * d2
            } else {
                d2 = f1.sinh() / nu
            }
            f0 = d2 - nu * d1 / (t1 * p0);
            if ex <= 1e-10 {
                /* ---------------------------------------------------------
                X <= 1.0E-10
                Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
                --------------------------------------------------------- */
                bk[0] = f0 + ex * f0;
                if *ize == 1 {
                    bk[0] -= ex * bk[0]
                }
                ratio = p0 / f0;
                c = ex * f64::MAX;
                if k != 0 {
                    /* ---------------------------------------------------
                    Calculation of K(ALPHA,X)
                    and  X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2
                    --------------------------------------------------- */
                    *ncalc = -(1);
                    if bk[0] >= c / ratio {
                        return;
                    }
                    bk[0] = ratio * bk[0] / ex;
                    twonu += 2.0;
                    ratio = twonu
                }
                *ncalc = 1;
                if *nb == 1 {
                    return;
                }
                /* -----------------------------------------------------
                Calculate  K(ALPHA+L,X)/K(ALPHA+L-1,X),
                L = 1, 2, ... , NB-1
                ----------------------------------------------------- */
                *ncalc = -(1);
                i = 1;
                while i < *nb {
                    if ratio >= c {
                        return;
                    }
                    bk[i] = ratio / ex;
                    twonu += 2.0;
                    ratio = twonu;
                    i += 1
                }
                *ncalc = 1;
                current_block_262 = 16458853316677622955;
            } else {
                /* ------------------------------------------------------
                10^-10 < X <= 1.0
                ------------------------------------------------------ */
                c = 1.0;
                x2by4 = ex * ex / 4.0;
                p0 = 0.5 * p0;
                q0 = 0.5 * q0;
                d1 = -1.0;
                d2 = 0.0;
                bk1 = 0.0;
                bk2 = 0.0;
                f1 = f0;
                f2 = p0;
                loop {
                    d1 += 2.0;
                    d2 += 1.0;
                    d3 = d1 + d3;
                    c = x2by4 * c / d2;
                    f0 = (d2 * f0 + p0 + q0) / d3;
                    p0 /= d2 - nu;
                    q0 /= d2 + nu;
                    t1 = c * f0;
                    t2 = c * (p0 - d2 * f0);
                    bk1 += t1;
                    bk2 += t2;
                    if !((t1 / (f1 + bk1)).abs() > 2.2204460492503131e-16
                        || (t2 / (f2 + bk2)).abs() > 2.2204460492503131e-16)
                    {
                        break;
                    }
                }
                bk1 = f1 + bk1;
                bk2 = 2.0 * (f2 + bk2) / ex;
                if *ize == 2 {
                    d1 = ex.exp();
                    bk1 *= d1;
                    bk2 *= d1
                }
                wminf = estf[0] * ex + estf[1];
                current_block_262 = 16185292562584120790;
            }
        } else {
            if 2.2204460492503131e-16 * ex > 1.0 {
                /* -------------------------------------------------
                X > 1./EPS
                ------------------------------------------------- */
                *ncalc = *nb as i32;
                bk1 = 1.0 / (strafe_consts::SQRT_2DPI * ex.sqrt());
                i = 0;
                while i < *nb {
                    bk[i] = bk1;
                    i += 1
                }
                return;
            } else {
                /* -------------------------------------------------------
                X > 1.0
                ------------------------------------------------------- */
                twox = ex + ex;
                blpha = 0.0;
                ratio = 0.0;
                if ex <= 4.0 {
                    /* ----------------------------------------------------------
                    Calculation of K(ALPHA+1,X)/K(ALPHA,X),  1.0 <= X <= 4.0
                    ----------------------------------------------------------*/
                    d2 = (estm[0] / ex + estm[1]).trunc();
                    m = d2 as usize;
                    d1 = d2 + d2;
                    d2 -= 0.5;
                    d2 *= d2;
                    i = 2;
                    while i <= m {
                        d1 -= 2.0;
                        d2 -= d1;
                        ratio = (d3 + d2) / (twox + d1 - ratio);
                        i += 1
                    }
                    /* -----------------------------------------------------------
                    Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward
                    recurrence and K(ALPHA,X) from the wronskian
                    -----------------------------------------------------------*/
                    d2 = (estm[2] * ex + estm[3]).trunc();
                    m = d2 as usize;
                    c = nu.abs();
                    d3 = c + c;
                    d1 = d3 - 1.0;
                    f1 = 2.2250738585072014e-308;
                    f0 =
                        (2.0 * (c + d2) / ex + 0.5 * ex / (c + d2 + 1.0)) * 2.2250738585072014e-308;
                    i = 3;
                    while i <= m {
                        d2 -= 1.0;
                        f2 = (d3 + d2 + d2) * f0;
                        blpha = (1.0 + d1 / d2) * (f2 + blpha);
                        f2 = f2 / ex + f1;
                        f1 = f0;
                        f0 = f2;
                        i += 1
                    }
                    f1 = (d3 + 2.0) * f0 / ex + f1;
                    d1 = 0.0;
                    t1 = 1.0;
                    i = 1;
                    while i <= 7 {
                        d1 = c * d1 + p[i - 1];
                        t1 = c * t1 + q[i - 1];
                        i += 1
                    }
                    p0 = (c * (a + c * (p[7] - c * d1 / t1) - ex.ln())).exp() / ex;
                    f2 = (c + 0.5 - ratio) * f1 / ex;
                    bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
                    if *ize == 1 {
                        bk1 *= (-ex).exp()
                    }
                    wminf = estf[2] * ex + estf[3]
                } else {
                    /* ---------------------------------------------------------
                    Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
                    backward recurrence, for  X > 4.0
                    ----------------------------------------------------------*/
                    dm = (estm[4] / ex + estm[5]).trunc();
                    m = dm as usize;
                    d2 = dm - 0.5;
                    d2 *= d2;
                    d1 = dm + dm;
                    i = 2;
                    while i <= m {
                        dm -= 1.0;
                        d1 -= 2.0;
                        d2 -= d1;
                        ratio = (d3 + d2) / (twox + d1 - ratio);
                        blpha = (ratio + ratio * blpha) / dm;
                        i += 1
                    }
                    bk1 = 1.0
                        / ((strafe_consts::SQRT_2DPI + strafe_consts::SQRT_2DPI * blpha)
                            * ex.sqrt());
                    if *ize == 1 {
                        bk1 *= (-ex).exp()
                    }
                    wminf = estf[4] * (ex - (ex - estf[6]).abs()) + estf[5]
                }
                /* ---------------------------------------------------------
                Calculation of K(ALPHA+1,X)
                from K(ALPHA,X) and  K(ALPHA+1,X)/K(ALPHA,X)
                --------------------------------------------------------- */
                bk2 = bk1 + bk1 * (nu + 0.5 - ratio) / ex
            }
            current_block_262 = 16185292562584120790;
        }
        match current_block_262 {
            16185292562584120790 => {
                /*--------------------------------------------------------------------
                Calculation of 'NCALC', K(ALPHA+I,X), I  =  0, 1, ... , NCALC-1,
                &   K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1
                -------------------------------------------------------------------*/
                *ncalc = *nb as i32;
                bk[0] = bk1;
                if iend == 0 {
                    return;
                }
                j = 1 - k;
                bk[j] = bk2;
                if iend == 1 {
                    return;
                }
                m = if (wminf - nu) <= iend as f64 {
                    (wminf - nu) as usize
                } else {
                    iend as usize
                };
                i = 2;
                while i <= m {
                    t1 = bk1;
                    bk1 = bk2;
                    twonu += 2.0;
                    if ex < 1.0 {
                        if bk1 >= f64::MAX / twonu * ex {
                            break;
                        }
                    } else if bk1 / ex >= f64::MAX / twonu {
                        break;
                    }
                    bk2 = twonu / ex * bk1 + t1;
                    ii = i;
                    j += 1;
                    bk[j] = bk2;
                    i += 1
                }
                m = ii;
                if m == iend {
                    return;
                }
                ratio = bk2 / bk1;
                mplus1 = m + 1;
                *ncalc = -(1);
                i = mplus1;
                while i <= iend {
                    twonu += 2.0;
                    ratio = twonu / ex + 1.0 / ratio;
                    j += 1;
                    if j >= 1 {
                        bk[j] = ratio
                    } else {
                        if bk2 >= f64::MAX / ratio {
                            return;
                        }
                        bk2 *= ratio
                    }
                    i += 1
                }
                *ncalc = if 1 <= mplus1 - k {
                    ((mplus1) - k) as i32
                } else {
                    1
                };
                if *ncalc == 1 {
                    bk[0] = bk2
                }
                if *nb == 1 {
                    return;
                }
            }
            _ => {}
        }
        i = *ncalc as usize;
        while i < *nb {
            /* i == *ncalc */
            bk[i] *= bk[i - 1];
            *ncalc += 1;
            i += 1
        }
    };
}