r2rs-nmath 0.1.1

// Translation of nmath's toms708
//
// 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_consts::{LN_SQRT_2TPI, SQRT_PI};

use crate::{distribution::func::logspace_add, traits::DPQ};

#[derive(Copy, Clone, Debug, Default)]
struct BRatioState {
    do_swap: bool,
    did_bup: bool,
    n: isize,
    ierr1: usize,
    z: f64,
    a0: f64,
    b0: f64,
    x0: f64,
    y0: f64,
    lambda: f64,
    eps: f64,
    exit: bool,
}

fn SET_0_noswap(state: &mut BRatioState, a: f64, x: f64, b: f64, y: f64) {
    state.a0 = a;
    state.x0 = x;
    state.b0 = b;
    state.y0 = y;
}

fn SET_0_swap(state: &mut BRatioState, a: f64, x: f64, b: f64, y: f64) {
    state.a0 = b;
    state.x0 = y;
    state.b0 = a;
    state.y0 = x;
}

/// Evaluation of the Incomplete Beta function $I_x(a,b)$
///
/// It is assumed that a and b are nonnegative, and that $x <= 1$
/// and $y = 1 - x$.  Bratio assigns $w$ and $w1$ the values
///
/// $w  = I_x(a,b)$
///
/// $w1 = 1 - I_x(a,b)$
///
/// ierr is a variable that reports the status of the results.
/// If no input errors are detected then ierr is set to 0 and
/// w and w1 are computed. otherwise, if an error is detected,
/// then w and w1 are assigned the value 0 and ierr is set to
/// one of the following values ...
///
/// ierr = 1  if a or b is negative
/// ierr = 2  if a = b = 0
/// ierr = 3  if x < 0 or x > 1
/// ierr = 4  if y < 0 or y > 1
/// ierr = 5  if x + y != 1
/// ierr = 6  if x = a = 0
/// ierr = 7  if y = b = 0
/// ierr = 8	(not used currently)
/// ierr = 9  NaN in a, b, x, or y
/// ierr = 10     (not used currently)
/// ierr = 11  bgrat() error code 1 \[+ warning in bgrat()\]
/// ierr = 12  bgrat() error code 2   (no warning here)
/// ierr = 13  bgrat() error code 3   (no warning here)
/// ierr = 14  bgrat() error code 4 \[+ WARNING in bgrat()\]
///
/// Based on C translation of ACM TOMS 708
/// Please do not change this, e.g. to use R's versions of the
/// ancillary routines, without investigating the error analysis as we
/// do need very high relative accuracy.  This version has about
/// 14 digits accuracy.
///
/// More specifically,  Brown & Levy (1994) "Certification of Algorithm 708" write
/// "
/// The number of significant digits of accuracy \[..\] was calculated \[..\] as
///
/// -log10 (2 RelativeError),
///
/// \[....\]
/// Accuracy ranged from 9.64 significant digits to 15.65 with a
/// median of 14.65 and a lower quartile of 13.81.
/// \[...\]
/// ... overall accuracy increases slightly as a/b moves away from 1.
/// Linear regression indicates that
/// (1) an average of 13.71 significant digits are obtained in cases in which a = b and
/// (2) the number increases 0.14 significant digits for each unit change in log10(a/b).
/// "
///
/// Written by Alfred H. Morris, Jr.
/// Naval Surface Warfare Center
/// Dahlgren, Virginia
/// Revised ... Nov 1991
pub fn bratio(
    a: f64,
    b: f64,
    x: f64,
    y: f64,
    w: &mut f64,
    w1: &mut f64,
    ierr: &mut usize,
    log: bool,
) {
    let mut state = BRatioState::default();
    state.ierr1 = 0;

    /*  eps is a machine dependent constant: the smallest
     *      floating point number for which   1. + eps > 1.
     * NOTE: for almost all purposes it is replaced by 1e-15 (~= 4.5 times larger) below */
    state.eps = 2.0 * (0.5 * f64::EPSILON); /* == DBL_EPSILON (in R, Rmath) */

    /* ----------------------------------------------------------------------- */
    *w = f64::d_0(log);
    *w1 = f64::d_0(log);

    // safeguard, preventing infinite loops further down
    if x.is_nan() || y.is_nan() || a.is_nan() || b.is_nan() {
        *ierr = 9;
        return;
    }

    if a < 0.0 || b < 0.0 {
        *ierr = 1;
        return;
    }
    if a == 0.0 && b == 0.0 {
        *ierr = 2;
        return;
    }
    if x < 0.0 || x > 1.0 {
        *ierr = 3;
        return;
    }
    if y < 0.0 || y > 1.0 {
        *ierr = 4;
        return;
    }

    /* check that  'y == 1 - x' : */
    state.z = x + y - 0.5 - 0.5;

    if state.z.abs() > state.eps * 3.0 {
        *ierr = 5;
        return;
    }

    // R_ifDEBUG_printf("bratio(a=%g, b=%g, x=%9g, y=%9g, .., log_p=%d): ",
    // a,b,x,y, log_p);
    *ierr = 0;
    if x == 0.0 {
        l200(&mut state, a, ierr);
        if state.exit {
            return;
        }
    }
    if y == 0.0 {
        l210(&mut state, b, ierr);
        if state.exit {
            return;
        }
    }

    if a == 0.0 {
        l211(&mut state, w, w1, log);
        if state.exit {
            return;
        }
    }
    if b == 0.0 {
        l201(&mut state, w, w1, log);
        if state.exit {
            return;
        }
    }

    state.eps = state.eps.max(1e-15);
    let a_lt_b = a < b;
    if
    /* max(a,b) */
    if a_lt_b { b } else { a } < state.eps * 0.001 {
        /* procedure for a and b < 0.001 * eps */
        // L230:  -- result *independent* of x (!)
        // *w  = a/(a+b)  and  w1 = b/(a+b) :
        if log {
            if a_lt_b {
                *w = (-a / (a + b)).ln_1p(); // notably if a << b
                *w1 = (a / (a + b)).ln();
            } else {
                // b <= a
                *w = (b / (a + b)).ln();
                *w1 = (-b / (a + b)).ln_1p();
            }
        } else {
            *w = b / (a + b);
            *w1 = a / (a + b);
        }

        // R_ifDEBUG_printf("a & b very small -> simple ratios (%g,%g)\n", *w,*w1);
        return;
    }

    if a.min(b) <= 1.0 {
        /*------------------------ a <= 1  or  b <= 1 ---- */

        state.do_swap = x > 0.5;
        if state.do_swap {
            SET_0_swap(&mut state, a, x, b, y);
        } else {
            SET_0_noswap(&mut state, a, x, b, y);
        }
        /* now have  x0 <= 1/2 <= y0  (still  x0+y0 == 1) */

        // R_ifDEBUG_printf(" min(a,b) <= 1, do_swap=%d;", do_swap);

        if state.b0 < state.eps.min(state.eps * state.a0) {
            /* L80: */
            *w = fpser(state.a0, state.b0, state.x0, state.eps, log);
            *w1 = if log { w.log1_exp() } else { 0.5 - *w + 0.5 };
            // R_ifDEBUG_printf("  b0 small -> w := fpser(*) = %.15g\n", *w);
            l_end(&mut state, w, w1);
            if state.exit {
                return;
            }
        }

        if state.a0 < state.eps.min(state.eps * state.b0) && state.b0 * state.x0 <= 1.0 {
            /* L90: */
            *w1 = apser(state.a0, state.b0, state.x0, state.eps);
            // R_ifDEBUG_printf("  a0 small -> w1 := apser(*) = %.15g\n", *w1);
            l_end_from_w1(&mut state, w, w1, log);
            if state.exit {
                return;
            }
        }

        state.did_bup = false;
        if state.a0.max(state.b0) > 1.0 {
            /* L20:  min(a,b) <= 1 < max(a,b)  */
            // R_ifDEBUG_printf("\n L20:  min(a,b) <= 1 < max(a,b); ");
            if state.b0 <= 1.0 {
                l_w_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }

            if state.x0 >= 0.29
            /* was 0.3, PR#13786 */
            {
                l_w1_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }

            if state.x0 < 0.1 && state.x0 * state.b0.powf(state.a0) <= 0.7 {
                l_w_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }

            if state.b0 > 15.0 {
                *w1 = 0.0;
                l131(&mut state, w, w1, ierr, log);
            }
        } else {
            /*  a, b <= 1 */
            // R_ifDEBUG_printf("\n      both a,b <= 1; ");
            if state.a0 >= 0.2.min(state.b0) {
                l_w_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }

            if state.x0.powf(state.a0) <= 0.9 {
                l_w_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }

            if state.x0 >= 0.3 {
                l_w1_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }
        }
        state.n = 20; /* goto L130; */
        *w1 = bup(
            state.b0, state.a0, state.y0, state.x0, state.n, state.eps, false,
        );
        state.did_bup = true;
        // R_ifDEBUG_printf("  ... n=20 and *w1 := bup(*) = %.15g; ", *w1);
        state.b0 += state.n as f64;
        l131(&mut state, w, w1, ierr, log);
        if state.exit {
            return;
        }
    } else {
        /* L30: -------------------- both  a, b > 1  {a0 > 1  &  b0 > 1} ---*/

        /* lambda := a y - b x  =  (a + b)y  =  a - (a+b)x    {using x + y == 1},
         * ------ using the numerically best version : */
        state.lambda = if (a + b).is_finite() {
            if a > b {
                (a + b) * y - b
            } else {
                a - (a + b) * x
            }
        } else {
            a * y - b * x
        };
        state.do_swap = state.lambda < 0.0;
        if state.do_swap {
            state.lambda = -state.lambda;
            SET_0_swap(&mut state, a, x, b, y);
        } else {
            SET_0_noswap(&mut state, a, x, b, y);
        }

        // R_ifDEBUG_printf("  L30:  both  a, b > 1; |lambda| = %#g, do_swap = %d\n",
        //                  lambda, do_swap);

        if state.b0 < 40.0 {
            // R_ifDEBUG_printf("  b0 < 40;");
            if state.b0 * state.x0 <= 0.7 || (log && state.lambda > 650.0) {
                // << added 2010-03; svn r51327
                l_w_bpser(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            } else {
                l140(&mut state, w, w1, ierr, log);
                if state.exit {
                    return;
                }
            }
        } else if state.a0 > state.b0 {
            /* ----  a0 > b0 >= 40  ---- */
            // R_ifDEBUG_printf("  a0 > b0 >= 40;");
            if state.b0 <= 100. || state.lambda > state.b0 * 0.03 {
                l_bfrac(&mut state, w, w1, log);
                if state.exit {
                    return;
                }
            }
        } else if state.a0 <= 100.0 {
            // R_ifDEBUG_printf("  a0 <= 100; a0 <= b0 >= 40;");
            l_bfrac(&mut state, w, w1, log);
            if state.exit {
                return;
            }
        } else if state.lambda > state.a0 * 0.03 {
            // R_ifDEBUG_printf("  b0 >= a0 > 100; lambda > a0 * 0.03 ");
            l_bfrac(&mut state, w, w1, log);
            if state.exit {
                return;
            }
        }

        /* else if none of the above    L180: */
        *w = basym(state.a0, state.b0, state.lambda, state.eps * 100.0, log);
        *w1 = if log { w.log1_exp() } else { 0.5 - *w + 0.5 };
        // R_ifDEBUG_printf("  b0 >= a0 > 100; lambda <= a0 * 0.03: *w:= basym(*) =%.15g\n",
        //                  *w);
        l_end(&mut state, w, w1);
        if state.exit {
            return;
        }
    } /* else: a, b > 1 */
    /*            EVALUATION OF THE APPROPRIATE ALGORITHM */

    l_w_bpser(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_w1_bpser(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_bfrac(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l140(&mut state, w, w1, ierr, log);
    if state.exit {
        return;
    }

    /* TERMINATION OF THE PROCEDURE */

    l200(&mut state, a, ierr);
    if state.exit {
        return;
    }
    // else:
    l201(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l210(&mut state, b, ierr);
    if state.exit {
        return;
    }
    // else:
    l211(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_end_from_w(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_end_from_w1(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_end_from_w1_log(&mut state, w, w1, log);
    if state.exit {
        return;
    }

    l_end(&mut state, w, w1);
    if state.exit {
        return;
    }
} /* bratio */

#[inline]
fn l131(state: &mut BRatioState, w: &mut f64, w1: &mut f64, ierr: &mut usize, log: bool) {
    // R_ifDEBUG_printf(" L131: bgrat(*, w1=%.15g) ", *w1);
    bgrat(
        state.b0,
        state.a0,
        state.y0,
        state.x0,
        w1,
        15.0 * state.eps,
        &mut state.ierr1,
        false,
    );
    // REprintf(" ==> new w1=%.15g", *w1);
    // if(ierr1) REprintf(" ERROR(code=%d)\n", ierr1) ; else REprintf("\n");
    if *w1 == 0.0 || (0.0 < *w1 && *w1 < f64::MIN) {
        // w1=0 or very close:
        // "almost surely" from underflow, try more: [2013-03-04]
        // FIXME: it is even better to do this in bgrat *directly* at least for the case
        //  !did_bup, i.e., where *w1 = (0 or -Inf) on entry
        //         R_ifDEBUG_printf(" denormalized or underflow (?) -> retrying: ");
        if state.did_bup {
            // re-do that part on log scale:
            *w1 = bup(
                state.b0 - state.n as f64,
                state.a0,
                state.y0,
                state.x0,
                state.n,
                state.eps,
                true,
            );
        } else {
            *w1 = f64::neg_infinity(); // = 0 on log-scale}
            bgrat(
                state.b0,
                state.a0,
                state.y0,
                state.x0,
                w1,
                15.0 * state.eps,
                &mut state.ierr1,
                true,
            );
            if state.ierr1 != 0 {
                *ierr = 10 + state.ierr1;
            }
            // REprintf(" ==> new log(w1)=%.15g", *w1);
            // if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n");
            l_end_from_w1_log(state, w, w1, log);
            if state.exit {
                return;
            }
        }
        // else
        if state.ierr1 != 0 {
            *ierr = 10 + state.ierr1;
        }
        // if *w1 < 0.0 {
        // MATHLIB_WARNING4("bratio(a=%g, b=%g, x=%g): bgrat() -> w1 = %g",
        //                  a,b,x, *w1);
    }
    l_end_from_w1(state, w, w1, log);
    if state.exit {
        return;
    }
}

#[inline]
fn l_w_bpser(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    *w = bpser(state.a0, state.b0, state.x0, state.eps, log);
    *w1 = if log { w.log1_exp() } else { 0.5 - *w + 0.5 };
    // R_ifDEBUG_printf(" L_w_bpser: *w := bpser(*) = %.15g\n", *w);
    l_end(state, w, w1);
}

#[inline]
fn l_w1_bpser(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    // was L110
    *w1 = bpser(state.b0, state.a0, state.y0, state.eps, log);
    *w = if log { w1.log1_exp() } else { 0.5 - *w1 + 0.5 };
    // R_ifDEBUG_printf(" L_w1_bpser: *w1 := bpser(*) = %.15g\n", *w1);
    l_end(state, w, w1);
}

#[inline]
fn l_bfrac(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    *w = bfrac(
        state.a0,
        state.b0,
        state.x0,
        state.y0,
        state.lambda,
        state.eps * 15.,
        log,
    );
    *w1 = if log { w.log1_exp() } else { 0.5 - *w + 0.5 };
    // R_ifDEBUG_printf(" L_bfrac: *w := bfrac(*) = %g\n", *w);
    l_end(state, w, w1);
}

#[inline]
fn l140(state: &mut BRatioState, w: &mut f64, w1: &mut f64, ierr: &mut usize, log: bool) {
    /* b0 := fractional_part( b0 )  in (0, 1]  */
    state.n = state.b0 as isize;
    state.b0 -= state.n as f64;
    if state.b0 == 0.0 {
        state.n -= 1;
        state.b0 = 1.0;
    }

    *w = bup(
        state.b0, state.a0, state.y0, state.x0, state.n, state.eps, false,
    );

    if *w < f64::min_value() && log {
        /* do not believe it; try bpser() : */
        // R_ifDEBUG_printf(" L140: bup(b0=%g,..)=%.15g < DBL_MIN - not used; ", b0, *w);
        /*revert: */
        state.b0 += state.n as f64;
        /* which is only valid if b0 <= 1 || b0*x0 <= 0.7 */
        l_w_bpser(state, w, w1, log);
        if state.exit {
            return;
        }
    }
    // R_ifDEBUG_printf(" L140: *w := bup(b0=%g,..) = %.15g; ", b0, *w);
    if state.x0 <= 0.7 {
        /* log_p :  TODO:  w = bup(.) + bpser(.)  -- not so easy to use log-scale */
        *w += bpser(
            state.a0, state.b0, state.x0, state.eps, /* log_p = */ false,
        );
        // R_ifDEBUG_printf(" x0 <= 0.7: *w := *w + bpser(*) = %.15g\n", *w);
        l_end_from_w(state, w, w1, log);
        if state.exit {
            return;
        }
    }
    /* L150: */
    if state.a0 <= 15.0 {
        state.n = 20;
        *w += bup(
            state.a0, state.b0, state.x0, state.y0, state.n, state.eps, false,
        );
        // R_ifDEBUG_printf("\n a0 <= 15: *w := *w + bup(*) = %.15g;", *w);
        state.a0 += state.n as f64;
    }
    // R_ifDEBUG_printf(" bgrat(*, w=%.15g) ", *w);
    bgrat(
        state.a0,
        state.b0,
        state.x0,
        state.y0,
        w,
        15.0 * state.eps,
        &mut state.ierr1,
        false,
    );
    if state.ierr1 != 0 {
        *ierr = 10 + state.ierr1;
    }

    // REprintf("==> new w=%.15g", *w);
    // if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n");

    l_end_from_w(state, w, w1, log);
}

#[inline]
fn l200(state: &mut BRatioState, a: f64, ierr: &mut usize) {
    if a == 0.0 {
        *ierr = 6;
        state.exit = true;
        return;
    }
}

#[inline]
fn l201(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    *w = f64::d_0(log);
    *w1 = f64::d_1(log);
    state.exit = true;
    return;
}

#[inline]
fn l210(state: &mut BRatioState, b: f64, ierr: &mut usize) {
    if b == 0.0 {
        *ierr = 7;
        state.exit = true;
        return;
    }
}

#[inline]
fn l211(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    *w = f64::d_1(log);
    *w1 = f64::d_0(log);
    state.exit = true;
    return;
}

#[inline]
fn l_end_from_w(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    if log {
        *w1 = (-*w).ln_1p();
        *w = w.ln();
    } else {
        *w1 = 0.5 - *w + 0.5;
    }
    l_end(state, w, w1)
}

#[inline]
fn l_end_from_w1(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    if log {
        *w = (-*w1).ln_1p();
        *w1 = (*w1).ln();
    } else {
        *w = 0.5 - *w1 + 0.5;
    }
    l_end(state, w, w1);
}

#[inline]
fn l_end_from_w1_log(state: &mut BRatioState, w: &mut f64, w1: &mut f64, log: bool) {
    // *w1 = log(w1) already; w = 1 - w1  ==> log(w) = log(1 - w1) = log(1 - exp(*w1))
    if log {
        *w = (*w1).log1_exp();
    } else {
        *w  = /* 1 - exp(*w1) */ -w1.exp_m1();
        *w1 = w1.exp();
    }
    l_end(state, w, w1);
}

#[inline]
fn l_end(state: &mut BRatioState, w: &mut f64, w1: &mut f64) {
    if state.do_swap {
        /* swap */
        let t = *w;
        *w = *w1;
        *w1 = t;
    }
    state.exit = true;
    return;
}

/// Evaluation of $I_{X}(A,B)$
///
/// For $B < \text{min}(eps, eps*A)$ and $X <= 0.5$
fn fpser(a: f64, b: f64, x: f64, eps: f64, log: bool) -> f64 {
    let mut ans;
    let mut c;
    let mut s;
    let mut t;
    let mut an;
    let tol;

    /* SET  ans := x^a : */
    if log {
        ans = a * x.ln();
    } else if a > eps * 0.001 {
        t = a * x.ln();
        if t < exparg(1) {
            /* exp(t) would underflow */
            return 0.0;
        }
        ans = t.exp();
    } else {
        ans = 1.0;
    }

    /*                NOTE THAT 1/B(A,B) = B */

    if log {
        ans += b.ln() - a.ln();
    } else {
        ans *= b / a;
    }

    tol = eps / a;
    an = a + 1.;
    t = x;
    s = t / an;
    loop {
        an += 1.;
        t = x * t;
        c = t / an;
        s += c;

        if !(c.abs() > tol) {
            break;
        }
    }

    if log {
        ans += (a * s).ln_1p();
    } else {
        ans *= a * s + 1.;
    }
    return ans;
} /* fpser */

///     apser() yields the incomplete beta ratio  I_{1-x}(b,a)  for
///     a <= min(eps,eps*b), b*x <= 1, and x <= 0.5,  i.e., a is very small.
///     Use only if above inequalities are satisfied.
fn apser(a: f64, b: f64, x: f64, eps: f64) -> f64 {
    let g = 0.577215664901533;

    let tol;
    let c;
    let mut j;
    let mut s;
    let mut t;
    let mut aj;
    let bx = b * x;

    t = x - bx;
    if b * eps <= 0.02 {
        c = x.ln() + psi(b) + g + t;
    } else {
        // b > 2e13 : psi(b) ~= log(b)
        c = bx.ln() + g + t;
    }

    tol = eps * 5. * c.abs();
    j = 1.;
    s = 0.;
    loop {
        j += 1.;
        t *= x - bx / j;
        aj = t / j;
        s += aj;

        if !(aj.abs() > tol) {
            break;
        }
    }

    return -a * (c + s);
} /* apser */

/// Power SERies expansion for evaluating I_x(a,b) when
///	       b <= 1 or b*x <= 0.7.   eps is the tolerance used.
/// NB: if log_p is true, also use it if   (b < 40  & lambda > 650)
fn bpser(a: f64, b: f64, x: f64, eps: f64, log: bool) -> f64 {
    let mut i;
    let m;
    let mut ans;
    let mut c;
    let t;
    let mut u;
    let z;
    let a0;
    let mut b0;
    let apb;

    if x == 0.0 {
        return f64::d_0(log);
    }
    /* ----------------------------------------------------------------------- */
    /*	      compute the factor  x^a/(a*Beta(a,b)) */
    /* ----------------------------------------------------------------------- */
    a0 = a.min(b);
    if a0 >= 1.0 {
        /*		 ------	 1 <= a0 <= b0  ------ */
        z = a * x.ln() - betaln(a, b);
        ans = if log { z - a.ln() } else { z.exp() / a };
    } else {
        b0 = a.max(b);

        if b0 < 8.0 {
            if b0 <= 1.0 {
                /*	 ------	 a0 < 1	 and  b0 <= 1  ------ */

                if log {
                    ans = a * x.ln();
                } else {
                    ans = x.powf(a);
                    if ans == 0.0 {
                        /* once underflow, always underflow .. */
                        return ans;
                    }
                }
                apb = a + b;
                if apb > 1.0 {
                    u = a + b - 1.0;
                    z = (gam1(u) + 1.0) / apb;
                } else {
                    z = gam1(apb) + 1.0;
                }
                c = (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;

                if log {
                    /* FIXME ? -- improve quite a bit for c ~= 1 */
                    ans += (c * (b / apb)).ln();
                } else {
                    ans *= c * (b / apb);
                }
            } else {
                /* 	------	a0 < 1 < b0 < 8	 ------ */

                u = gamln1(a0);
                m = (b0 - 1.0) as isize;
                if m >= 1 {
                    c = 1.0;

                    i = 1;
                    while i <= m as i32 {
                        b0 += -1.0;
                        c *= b0 / (a0 + b0);

                        i += 1;
                    }
                    u += c.ln();
                }

                z = a * x.ln() - u;
                b0 += -1.; // => b0 in (0, 7)
                apb = a0 + b0;
                if apb > 1.0 {
                    u = a0 + b0 - 1.;
                    t = (gam1(u) + 1.) / apb;
                } else {
                    t = gam1(apb) + 1.;
                }

                if log {
                    /* FIXME? potential for improving log(t) */
                    ans = z + (a0 / a).ln() + gam1(b0).ln_1p() - t.ln();
                } else {
                    ans = z.exp() * (a0 / a) * (gam1(b0) + 1.) / t;
                }
            }
        } else {
            /* 		------  a0 < 1 < 8 <= b0  ------ */

            u = gamln1(a0) + algdiv(a0, b0);
            z = a * x.ln() - u;

            if log {
                ans = z + (a0 / a).ln();
            } else {
                ans = a0 / a * z.exp();
            }
        }
    }
    // R_ifDEBUG_printf(" bpser(a=%g, b=%g, x=%g, log=%d): prelim.ans = %.14g;\n",
    // a,b,x, log_p, ans);
    if ans == f64::d_0(log) || (!log && a <= eps * 0.1) {
        return ans;
    }

    /* ----------------------------------------------------------------------- */
    /*		       COMPUTE THE SERIES */
    /* ----------------------------------------------------------------------- */
    let tol = eps / a;
    let mut n = 0.0;
    let mut sum = 0.0;
    let mut w;
    c = 1.;
    loop {
        // sum is alternating as long as n < b (<==> 1 - b/n < 0)
        n += 1.0;
        c *= (0.5 - b / n + 0.5) * x;
        w = c / (a + n);
        sum += w;

        if !(n < 1e7 && w.abs() > tol) {
            break;
        }
    }

    // if fabs(w) > tol) { // the series did not converge (in time)
    // // warn only when the result seems to matter:
    // if(( log_p && !(a*sum > -1. && fabs(log1p(a * sum)) < eps*fabs(ans))) ||
    // (!log_p && fabs(a*sum + 1.) != 1.))
    // MATHLIB_WARNING5(
    // " bpser(a=%g, b=%g, x=%g,...) did not converge (n=1e7, |w|/tol=%g > 1; A=%g)",
    // a,b,x, fabs(w)/tol, ans);
    // }

    // R_ifDEBUG_printf("  -> n=%.0f iterations, |w|=%g %s %g=tol:=eps/a ==> a*sum=%g\n",
    // n, fabs(w), (fabs(w) > tol) ? ">!!>" : "<=",
    // tol, a*sum);

    if log {
        if a * sum > -1.0 {
            ans += (a * sum).ln_1p();
        } else {
            if ans > f64::neg_infinity() {
                // MATHLIB_WARNING3(
                // "pbeta(*, log.p=true) -> bpser(a=%g, b=%g, x=%g,...) underflow to -Inf",
                // a,b,x);
                ans = f64::neg_infinity();
            }
        }
    } else if a * sum > -1.0 {
        ans *= a * sum + 1.0;
    } else {
        // underflow to
        ans = 0.;
    }
    return ans;
} /* bpser */

///    EVALUATION OF I_x(A,B) - I_x(A+N,B) WHERE N IS A POSITIVE INT.
///     EPS IS THE TOLERANCE USED.
fn bup(a: f64, b: f64, x: f64, y: f64, n: isize, eps: f64, give_log: bool) -> f64 {
    let mut ret_val;
    let mut i;
    let mut k;
    let mut mu;
    let mut d;
    let mut l;

    // Obtain the scaling factor exp(-mu) and exp(mu)*(x^a * y^b / beta(a,b))/a

    let apb = a + b;
    let ap1 = a + 1.;
    if n > 1 && a >= 1. && apb >= ap1 * 1.1 {
        mu = (exparg(1)).abs() as isize;
        k = exparg(0) as isize;
        if mu > k {
            mu = k;
        }
        d = (-(mu as f64)).exp();
    } else {
        mu = 0;
        d = 1.;
    }

    /* L10: */
    ret_val = if give_log {
        brcmp1(mu, a, b, x, y, true) - a.ln()
    } else {
        brcmp1(mu, a, b, x, y, false) / a
    };
    if n == 1 || (give_log && ret_val == f64::neg_infinity()) || (!give_log && ret_val == 0.) {
        return ret_val;
    }

    let nm1 = n - 1;
    let mut w = d;

    /*          LET K BE THE INDEX OF THE MAXIMUM TERM */

    k = 0;
    if b > 1.0 {
        if y > 1e-4 {
            let r = (b - 1.0) * x / y - a;
            if r >= 1.0 {
                k = if r < nm1 as f64 { r as isize } else { nm1 };
            }
        } else {
            k = nm1;
        }

        //          ADD THE INCREASING TERMS OF THE SERIES - if k > 0
        /* L30: */
        i = 0.0;
        while i < k as f64 {
            l = i;
            d *= (apb + l) / (ap1 + l) * x;
            w += d;

            i += 1.0
        }
    }

    // L40:     ADD THE REMAINING TERMS OF THE SERIES
    i = k as f64;
    while i < nm1 as f64 {
        l = i as f64;
        d *= (apb + l) / (ap1 + l) * x;
        w += d;
        if d <= eps * w {
            /* relativ convergence (eps) */
            break;
        }

        i += 1.0
    }

    // L50: TERMINATE THE PROCEDURE
    if give_log {
        ret_val += w.ln();
    } else {
        ret_val *= w;
    }

    return ret_val;
} /* bup */

/// Continued fraction expansion for I_x(a,b) when a, b > 1.
/// It is assumed that  lambda = (a + b)*y - b.
fn bfrac(a: f64, b: f64, x: f64, y: f64, lambda: f64, eps: f64, log: bool) -> f64 {
    let c;
    let mut e;
    let mut n;
    let mut p;
    let mut r;
    let mut s;
    let mut t;
    let mut w;
    let c0;
    let c1;
    let mut r0;
    let mut an;
    let mut bn;
    let yp1;
    let mut anp1;
    let mut bnp1;
    let mut beta;
    let mut alpha;
    let brc;

    if !lambda.is_finite() {
        return f64::nan();
    } // TODO: can return 0 or 1 (?)
      // R_ifDEBUG_printf(" bfrac(a=%g, b=%g, x=%g, y=%g, lambda=%g, eps=%g, log_p=%d):",
      // a,b,x,y, lambda, eps, log_p);
    brc = brcomp(a, b, x, y, log);
    if brc.is_nan() {
        // e.g. from   L <- 1e308; pnbinom(L, L, mu = 5)
        // R_ifDEBUG_printf(" --> brcomp(a,b,x,y) = NaN\n");
        // ML_WARN_return_NAN; // TODO: could we know better?
        return f64::nan();
    }
    if !log && brc == 0.0 {
        // R_ifDEBUG_printf(" --> brcomp(a,b,x,y) underflowed to 0.\n");
        return 0.0;
    }

    // else
    // REprintf("\n");

    c = lambda + 1.0;
    c0 = b / a;
    c1 = 1.0 / a + 1.0;
    yp1 = y + 1.0;

    n = 0.0;
    p = 1.0;
    s = a + 1.0;
    an = 0.0;
    bn = 1.0;
    anp1 = 1.0;
    bnp1 = c / c1;
    r = c1 / c;

    /*        CONTINUED FRACTION CALCULATION */

    loop {
        n += 1.0;
        t = n / a;
        w = n * (b - n) * x;
        e = a / s;
        alpha = p * (p + c0) * e * e * (w * x);
        e = (t + 1.0) / (c1 + t + t);
        beta = n + w / s + e * (c + n * yp1);
        p = t + 1.0;
        s += 2.0;

        /* update an, bn, anp1, and bnp1 */

        t = alpha * an + beta * anp1;
        an = anp1;
        anp1 = t;
        t = alpha * bn + beta * bnp1;
        bn = bnp1;
        bnp1 = t;

        r0 = r;
        r = anp1 / bnp1;

        // R_ifDEBUG_printf(" n=%5.0f, a_{n,n+1}= (%12g,%12g),  b_{n,n+1} = (%12g,%12g) => r0,r = (%14g,%14g)\n",
        // n, an,anp1, bn,bnp1, r0, r);

        if (r - r0).abs() <= eps * r {
            break;
        }

        /* rescale an, bn, anp1, and bnp1 */

        an /= bnp1;
        bn /= bnp1;
        anp1 = r;
        bnp1 = 1.0;

        if !(n < 10000.0) {
            break;
        }
    } // arbitrary; had '1' --> infinite loop for  lambda = Inf
      // R_ifDEBUG_printf("  in bfrac(): n=%.0f terms cont.frac.; brc=%g, r=%g\n",
      // n, brc, r);
      // if(n >= 10000 && fabs(r - r0) > eps * r)
      // MATHLIB_WARNING5(
      // " bfrac(a=%g, b=%g, x=%g, y=%g, lambda=%g) did *not* converge (in 10000 steps)\n",
      // a,b,x,y, lambda);
    return if log { brc + r.ln() } else { brc * r };
} /* bfrac */

fn brcomp(a: f64, b: f64, x: f64, y: f64, log: bool) -> f64 {
    /* -----------------------------------------------------------------------
     *		 Evaluation of x^a * y^b / Beta(a,b)
     * ----------------------------------------------------------------------- */

    let const__ = 0.398942280401433; /* == 1/sqrt(2*pi); */
    /* R has  M_1_SQRT_2PI , and M_LN_SQRT_2PI = ln(sqrt(2*pi)) = 0.918938.. */
    let mut i;
    let n;
    let mut c;
    let mut e;
    let mut u;
    let v;
    let mut z;
    let a0;
    let mut b0;
    let apb;

    if x == 0.0 || y == 0.0 {
        return f64::d_0(log);
    }
    a0 = a.min(b);
    if a0 < 8.0 {
        let lnx;
        let lny;
        if x <= 0.375 {
            lnx = x.ln();
            lny = alnrel(-x);
        } else {
            if y > 0.375 {
                lnx = x.ln();
                lny = y.ln();
            } else {
                lnx = alnrel(-y);
                lny = y.ln();
            }
        }

        z = a * lnx + b * lny;
        if a0 >= 1.0 {
            z -= betaln(a, b);
            return z.d_exp(log);
        }

        /* ----------------------------------------------------------------------- */
        /*		PROCEDURE FOR a < 1 OR b < 1 */
        /* ----------------------------------------------------------------------- */

        b0 = a.max(b);
        if b0 >= 8.0 {
            /* L80: */
            u = gamln1(a0) + algdiv(a0, b0);

            return if log {
                a0.ln() + (z - u)
            } else {
                a0 * (z - u).exp()
            };
        }
        /* else : */

        if b0 <= 1.0 {
            /*		algorithm for max(a,b) = b0 <= 1 */

            let e_z = z.d_exp(log);

            if !log && e_z == 0.0 {
                /* exp() underflow */
                return 0.0;
            }

            apb = a + b;
            if apb > 1.0 {
                u = a + b - 1.0;
                z = (gam1(u) + 1.0) / apb;
            } else {
                z = gam1(apb) + 1.0;
            }

            c = (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;
            /* FIXME? log(a0*c)= log(a0)+ log(c) and that is improvable */
            return if log {
                e_z + (a0 * c).ln() - (a0 / b0).ln_1p()
            } else {
                e_z * (a0 * c) / (a0 / b0 + 1.0)
            };
        }

        /* else : 		  ALGORITHM FOR 1 < b0 < 8 */

        u = gamln1(a0);
        n = (b0 - 1.0) as usize;
        if n >= 1 {
            c = 1.0;

            i = 1;
            while i <= n {
                b0 += -1.0;
                c *= b0 / (a0 + b0);

                i += 1;
            }
            u = c.ln() + u;
        }
        z -= u;
        b0 += -1.0;
        apb = a0 + b0;
        let t;
        if apb > 1.0 {
            u = a0 + b0 - 1.0;
            t = (gam1(u) + 1.0) / apb;
        } else {
            t = gam1(apb) + 1.0;
        }

        return if log {
            a0.ln() + z + (gam1(b0)).ln_1p() - t.ln()
        } else {
            a0 * z.exp() * (gam1(b0) + 1.) / t
        };
    } else {
        /* ----------------------------------------------------------------------- */
        /*		PROCEDURE FOR A >= 8 AND B >= 8 */
        /* ----------------------------------------------------------------------- */
        let h;
        let x0;
        let y0;
        let lambda;
        if a <= b {
            h = a / b;
            x0 = h / (h + 1.0);
            y0 = 1.0 / (h + 1.0);
            lambda = a - (a + b) * x;
        } else {
            h = b / a;
            x0 = 1.0 / (h + 1.0);
            y0 = h / (h + 1.0);
            lambda = (a + b) * y - b;
        }

        e = -lambda / a;
        if e.abs() > 0.6 {
            u = e - (x / x0).ln();
        } else {
            u = rlog1(e);
        }

        e = lambda / b;
        if e.abs() <= 0.6 {
            v = rlog1(e);
        } else {
            v = e - (y / y0).ln();
        }

        z = if log {
            -(a * u + b * v)
        } else {
            -(a * u + b * v).exp()
        };

        return if log {
            -LN_SQRT_2TPI + 0.5 * (b * x0).ln() + z - bcorr(a, b)
        } else {
            const__ * (b * x0).sqrt() * z * (-bcorr(a, b)).exp()
        };
    }
} /* brcomp */

/// Evaluation of    exp(mu) * x^a * y^b / beta(a,b)
/// called only once from  bup(),  as   r = brcmp1(mu, a, b, x, y, false) / a;
fn brcmp1(mu: isize, a: f64, b: f64, x: f64, y: f64, give_log: bool) -> f64 {
    let const__ = 0.398942280401433; /* == 1/sqrt(2*pi); */
    /* R has  M_1_SQRT_2PI */

    /* Local variables */
    let mut c;
    let t;
    let mut u;
    let v;
    let mut z;
    let a0;
    let mut b0;
    let apb;

    a0 = a.min(b);
    if a0 < 8.0 {
        let lnx;
        let lny;
        if x <= 0.375 {
            lnx = x.ln();
            lny = alnrel(-x);
        } else if y > 0.375 {
            // L11:
            lnx = x.ln();
            lny = y.ln();
        } else {
            lnx = alnrel(-y);
            lny = y.ln();
        }

        // L20:
        z = a * lnx + b * lny;
        if a0 >= 1.0 {
            z -= betaln(a, b);
            return esum(mu, z, give_log);
        }
        // else :
        /* ----------------------------------------------------------------------- */
        /*              PROCEDURE FOR A < 1 OR B < 1 */
        /* ----------------------------------------------------------------------- */
        // L30:
        b0 = a.max(b);
        if b0 >= 8.0 {
            /* L80:                  ALGORITHM FOR b0 >= 8 */
            u = gamln1(a0) + algdiv(a0, b0);
            // R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1, b0 >= 8;  z=%.15g\n", z);
            return if give_log {
                a0.ln() + esum(mu, z - u, true)
            } else {
                a0 * esum(mu, z - u, false)
            };
        } else if b0 <= 1.0 {
            //                   a0 < 1, b0 <= 1
            let ans = esum(mu, z, give_log);
            if ans == if give_log { f64::neg_infinity() } else { 0.0 } {
                return ans;
            }

            apb = a + b;
            if apb > 1.0 {
                // L40:
                u = a + b - 1.0;
                z = (gam1(u) + 1.0) / apb;
            } else {
                z = gam1(apb) + 1.0;
            }
            // L50:
            c = if give_log {
                gam1(a).ln_1p() + gam1(b).ln_1p() - z.ln()
            } else {
                (gam1(a) + 1.0) * (gam1(b) + 1.0) / z
            };
            // R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1, b0 <= 1;  c=%.15g\n", c);
            return if give_log {
                ans + a0.ln() + c - (a0 / b0).ln_1p()
            } else {
                ans * (a0 * c) / (a0 / b0 + 1.)
            };
        }
        // else:               algorithm for	a0 < 1 < b0 < 8
        // L60:
        u = gamln1(a0);
        let n = (b0 - 1.0) as isize;
        if n >= 1 {
            c = 1.0;

            let mut i = 1;
            while i <= n {
                b0 += -1.0;
                c *= b0 / (a0 + b0);
                /* L61: */

                i += 1;
            }
            u += c.ln(); // TODO?: log(c) = log( prod(...) ) =  sum( log(...) )
        }
        // L70:
        z -= u;
        b0 += -1.0;
        apb = a0 + b0;
        if apb > 1.0 {
            // L71:
            t = (gam1(apb - 1.0) + 1.0) / apb;
        } else {
            t = gam1(apb) + 1.0;
        }
        // R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1 < b0 < 8;  t=%.15g\n", t);
        // L72:
        return if give_log {
            a0.ln() + esum(mu, z, true) + gam1(b0).ln_1p() - t.ln() // TODO? log(t) = log1p(..)
        } else {
            a0 * esum(mu, z, false) * (gam1(b0) + 1.0) / t
        };
    } else {
        /* ----------------------------------------------------------------------- */
        /*              PROCEDURE FOR A >= 8 AND B >= 8 */
        /* ----------------------------------------------------------------------- */
        // L100:
        let h;
        let x0;
        let y0;
        let lambda;
        if a > b {
            // L101:
            h = b / a;
            x0 = 1.0 / (h + 1.0); // => lx0 := log(x0) = 0 - log1p(h)
            y0 = h / (h + 1.0);
            lambda = (a + b) * y - b;
        } else {
            h = a / b;
            x0 = h / (h + 1.0); // => lx0 := log(x0) = - log1p(1/h)
            y0 = 1.0 / (h + 1.0);
            lambda = a - (a + b) * x;
        }
        let lx0 = -(b / a).ln_1p(); // in both cases

        // R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a,b >= 8;	x0=%.15g, lx0=log(x0)=%.15g\n",
        // x0, lx0);
        // L110:
        let mut e = -lambda / a;
        if e.abs() > 0.6 {
            // L111:
            u = e - (x / x0).ln();
        } else {
            u = rlog1(e);
        }

        // L120:
        e = lambda / b;
        if e.abs() > 0.6 {
            // L121:
            v = e - (y / y0).ln();
        } else {
            v = rlog1(e);
        }

        // L130:
        z = esum(mu, -(a * u + b * v), give_log);
        return if give_log {
            const__.ln() + (b.ln() + lx0) / 2.0 + z - bcorr(a, b)
        } else {
            const__ * (b * x0).sqrt() * z * (-bcorr(a, b)).exp()
        };
    }
} /* brcmp1 */

/// Asymptotic Expansion for I_x(a,b)  when a is larger than b.
///     Compute   w := w + I_x(a,b)
///     It is assumed a >= 15 and b <= 1.
///     eps is the tolerance used.
///     ierr is a variable that reports the status of the results.
///
/// if(log_w),  *w  itself must be in log-space;
///     compute   w := w + I_x(a,b)  but return *w = log(w):
///          *w := log(exp(*w) + I_x(a,b)) = logspace_add(*w, log( I_x(a,b) ))
fn bgrat(a: f64, b: f64, x: f64, y: f64, w: &mut f64, eps: f64, ierr: &mut usize, log: bool) {
    let n_terms_bgrat = 30;
    let mut c = vec![0.0; n_terms_bgrat];
    let mut d = vec![0.0; n_terms_bgrat];
    let bm1 = b - 0.5 - 0.5;
    let nu = a + bm1 * 0.5; /* nu = a + (b-1)/2 =: T, in (9.1) of
                             * Didonato & Morris(1992), p.362 */
    let lnx = if y > 0.375 { x.ln() } else { alnrel(-y) };
    let z = -nu * lnx; // z =: u in (9.1) of D.&M.(1992)

    if b * z == 0.0 {
        // should not happen, but does, e.g.,
        // for  pbeta(1e-320, 1e-5, 0.5)  i.e., _subnormal_ x,
        // Warning ... bgrat(a=20.5, b=1e-05, x=1, y=9.99989e-321): ..
        // MATHLIB_WARNING5(
        // "bgrat(a=%g, b=%g, x=%g, y=%g): z=%g, b*z == 0 underflow, hence inaccurate pbeta()",
        // a,b,x,y, z);
        /* L_Error:    THE EXPANSION CANNOT BE COMPUTED */
        *ierr = 1;
        return;
    }

    /*                 COMPUTATION OF THE EXPANSION */

    /* r1 = b * (gam1(b) + 1.) * exp(b * log(z)),// = b/gamma(b+1) z^b = z^b / gamma(b)
     * set r := exp(-z) * z^b / gamma(b) ;
     *          gam1(b) = 1/gamma(b+1) - 1 , b in [-1/2, 3/2] */
    // exp(a*lnx) underflows for large (a * lnx); e.g. large a ==> using log_r := log(r):
    // r = r1 * exp(a * lnx) * exp(bm1 * 0.5 * lnx);
    // log(r)=log(b) + log1p(gam1(b)) + b * log(z) + (a * lnx) + (bm1 * 0.5 * lnx),
    let log_r = b.ln() + gam1(b).ln_1p() + b * z.ln() + nu * lnx;
    // FIXME work with  log_u = log(u)  also when log_p=false  (??)
    // u is 'factored out' from the expansion {and multiplied back, at the end}:
    let log_u = log_r - (algdiv(b, a) + b * nu.ln()); // algdiv(b,a) = log(gamma(a)/gamma(a+b))
                                                      /* u = (log_p) ? log_r - u : exp(log_r-u); // =: M  in (9.2) of {reference above} */
    /* u = algdiv(b, a) + b * log(nu);// algdiv(b,a) = log(gamma(a)/gamma(a+b)) */
    // u = (log_p) ? log_u : exp(log_u); // =: M  in (9.2) of {reference above}
    let u = log_u.exp();

    if log_u == f64::neg_infinity() {
        // R_ifDEBUG_printf(" bgrat(*): underflow log_u = -Inf  = log_r -u', log_r = %g ",
        // log_r);
        /* L_Error:    THE EXPANSION CANNOT BE COMPUTED */
        *ierr = 2;
        return;
    }

    let u_0 = u == 0.0; // underflow --> do work with log(u) == log_u !
    let l = // := *w/u .. but with care: such that it also works when u underflows to 0:
if log
{
    if *w == f64::neg_infinity() {
    0.0} else {
        (*w - log_u).exp()
    }
}else{ if *w == 0.0 { 0.0 }else{ (w.ln() - log_u).exp()}};

    // R_ifDEBUG_printf(" bgrat(a=%g, b=%g, x=%g, *)\n -> u=%g, l='w/u'=%g, ",
    // a,b,x, u, l);

    let q_r = grat_r(b, z, log_r, eps); // = q/r of former grat1(b,z, r, &p, &q)
    let v = 0.25 / (nu * nu);
    let t2 = lnx * 0.25 * lnx;
    let mut j = q_r;
    let mut sum = j;
    let mut t = 1.0;
    let mut cn = 1.0;
    let mut n2 = 0.0;

    let mut n = 1;
    while n <= n_terms_bgrat {
        let bp2n = b + n2;
        j = (bp2n * (bp2n + 1.0) * j + (z + bp2n + 1.0) * t) * v;
        n2 += 2.0;
        t *= t2;
        cn /= n2 * (n2 + 1.0);
        let nm1 = n - 1;
        c[nm1] = cn;
        let mut s = 0.0;
        if n > 1 {
            let mut coef = b - n as f64;

            let mut i = 1;
            while i <= nm1 {
                s += coef * c[i - 1] * d[nm1 - i];
                coef += b;

                i += 1;
            }
        }
        d[nm1] = bm1 * cn + s / n as f64;
        let dj = d[nm1] * j;
        sum += dj;
        if sum <= 0.0 {
            // R_ifDEBUG_printf(" bgrat(*): sum_n(..) <= 0; should not happen (n=%d)\n", n);
            /* L_Error:    THE EXPANSION CANNOT BE COMPUTED */
            *ierr = 3;
            return;
        }
        if dj.abs() <= eps * (sum + l) {
            *ierr = 0;
            break;
        } else if n == n_terms_bgrat {
            // never? ; please notify R-core if seen:
            *ierr = 4;
            // MATHLIB_WARNING5(
            // "bgrat(a=%g, b=%g, x=%g) *no* convergence: NOTIFY R-core!\n dj=%g, rel.err=%g\n",
            // a,b,x, dj, fabs(dj) /(sum + l));
        }

        n += 1;
    } // for(n .. n_terms..)

    /*                    ADD THE RESULTS TO W */

    if log {
        // *w is in log space already:
        *w = logspace_add(*w, log_u + sum.ln());
    } else {
        *w += if u_0 {
            (log_u + sum.ln()).exp()
        } else {
            u * sum
        };
    }
    return;
} /* bgrat */

/// called only from bgrat() , as   q_r = grat_r(b, z, log_r, eps)  :
///        Scaled complement of incomplete gamma ratio function
///                   grat_r(a,x,r) :=  Q(a,x) / r
/// where
///               Q(a,x) = pgamma(x,a, lower.tail=false)
///     and            r = e^(-x)* x^a / Gamma(a) ==  exp(log_r)
///
/// It is assumed that a <= 1.  eps is the tolerance to be used.
fn grat_r(a: f64, x: f64, log_r: f64, eps: f64) -> f64 {
    if a * x == 0.0 {
        /* L130: */
        if x <= a {
            /* L100: */
            return (-log_r).exp();
        } else {
            /* L110:*/
            return 0.0;
        }
    } else if a == 0.5 {
        // e.g. when called from pt()
        /* L120: */
        if x < 0.25 {
            let p = erf__(x.sqrt());
            // R_ifDEBUG_printf(" grat_r(a=%g, x=%g ..)): a=1/2 --> p=erf__(.)= %g\n",
            // a, x, p);
            return (0.5 - p + 0.5) * (-log_r).exp();
        } else {
            // 2013-02-27: improvement for "large" x: direct computation of q/r:
            let sx = x.sqrt();
            let q_r = erfc1(1, sx) / sx * SQRT_PI;
            // R_ifDEBUG_printf(" grat_r(a=%g, x=%g ..)): a=1/2 --> q_r=erfc1(..)/r= %g\n",
            // a,x, q_r);
            return q_r;
        }
    } else if x < 1.1 {
        /* L10:  Taylor series for  P(a,x)/x^a */

        let mut an = 3.0;
        let mut c = x;
        let mut sum = x / (a + 3.0);
        let tol = eps * 0.1 / (a + 1.0);
        let mut t;
        loop {
            an += 1.0;
            c *= -(x / an);
            t = c / (a + an);
            sum += t;

            if !(t.abs() > tol) {
                break;
            }
        }

        // R_ifDEBUG_printf(" grat_r(a=%g, x=%g, log_r=%g): sum=%g; Taylor w/ %.0f terms",
        // a,x,log_r, sum, an-3.);
        let j = a * x * ((sum / 6.0 - 0.5 / (a + 2.0)) * x + 1.0 / (a + 1.0));
        let z = a * x.ln();
        let h = gam1(a);
        let g = h + 1.0;

        if (x >= 0.25 && (a < x / 2.59)) || (z > -0.13394) {
            // L40:
            let l = rexpm1(z);
            let q = ((l + 0.5 + 0.5) * j - l) * g - h;
            if q <= 0.0 {
                // R_ifDEBUG_printf(" => q_r= 0.\n");
                /* L110:*/
                return 0.;
            } else {
                // R_ifDEBUG_printf(" => q_r=%.15g\n", q * exp(-log_r));
                return q * (-log_r).exp();
            }
        } else {
            let p = z.exp() * g * (0.5 - j + 0.5);
            // R_ifDEBUG_printf(" => q_r=%.15g\n", (0.5 - p + 0.5) * exp(-log_r));
            return /* q/r = */ (0.5 - p + 0.5) * (-log_r).exp();
        }
    } else {
        /* L50: ----  (x >= 1.1)  ---- Continued Fraction Expansion */

        let mut a2n_1 = 1.0;
        let mut a2n = 1.0;
        let mut b2n_1 = x;
        let mut b2n = x + (1.0 - a);
        let mut c = 1.0;
        let mut am0;
        let mut an0;

        loop {
            a2n_1 = x * a2n + c * a2n_1;
            b2n_1 = x * b2n + c * b2n_1;
            am0 = a2n_1 / b2n_1;
            c += 1.0;
            let c_a = c - a;
            a2n = a2n_1 + c_a * a2n;
            b2n = b2n_1 + c_a * b2n;
            an0 = a2n / b2n;

            if !((an0 - am0).abs() >= eps * an0) {
                break;
            }
        }

        // R_ifDEBUG_printf(" grat_r(a=%g, x=%g, log_r=%g): Cont.frac. %.0f terms => q_r=%.15g\n",
        // a,x, log_r, c-1., an0);
        return /* q/r = (r * an0)/r = */ an0;
    }
} /* grat_r */

///     ASYMPTOTIC EXPANSION FOR I_x(A,B) FOR LARGE A AND B. */
///     LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED. */
///     IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT */
///     A AND B ARE GREATER THAN OR EQUAL TO 15. */
fn basym(a: f64, b: f64, lambda: f64, eps: f64, log: bool) -> f64 {
    /* ------------------------ */
    /*     ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP */
    /*            ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. */
    let num_IT = 20;
    /*            THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1. */

    let e0 = 1.12837916709551; /* e0 == 2/sqrt(pi) */
    let e1 = 0.353553390593274; /* e1 == 2^(-3/2)   */
    let ln_e0 = 0.120782237635245; /* == ln(e0) */

    let mut a0 = vec![0.0; num_IT + 1];
    let mut b0 = vec![0.0; num_IT + 1];
    let mut c = vec![0.0; num_IT + 1];
    let mut d = vec![0.0; num_IT + 1];

    let f = a * rlog1(-lambda / a) + b * rlog1(lambda / b);
    let t;
    if log {
        t = -f;
    } else {
        t = (-f).exp();
        if t == 0.0 {
            return 0.0; /* once underflow, always underflow .. */
        }
    }
    let z0 = f.sqrt();
    let z = z0 / e1 * 0.5;
    let z2 = f + f;
    let h;
    let r0;
    let r1;
    let w0;

    if a < b {
        h = a / b;
        r0 = 1.0 / (h + 1.0);
        r1 = (b - a) / b;
        w0 = 1.0 / (a * (h + 1.0)).sqrt();
    } else {
        h = b / a;
        r0 = 1.0 / (h + 1.0);
        r1 = (b - a) / a;
        w0 = 1.0 / (b * (h + 1.0)).sqrt();
    }

    a0[0] = r1 * 0.66666666666666663;
    c[0] = a0[0] * -0.5;
    d[0] = -c[0];
    let mut j0 = 0.5 / e0 * erfc1(1, z0);
    let mut j1 = e1;
    let mut sum = j0 + d[0] * w0 * j1;

    let mut s = 1.0;
    let h2 = h * h;
    let mut hn = 1.0;
    let mut w = w0;
    let mut znm1 = z;
    let mut zn = z2;

    let mut n = 2;
    while n <= num_IT {
        hn *= h2;
        a0[n - 1] = r0 * 2.0 * (h * hn + 1.0) / (n as f64 + 2.0);
        let np1 = n + 1;
        s += hn;
        a0[np1 - 1] = r1 * 2.0 * s / (n as f64 + 3.0);

        let mut i = n;
        while i <= np1 {
            let r = (i as f64 + 1.0) * -0.5;
            b0[0] = r * a0[0];

            let mut m = 2;
            while m <= i {
                let mut bsum = 0.0;

                let mut j = 1;
                while j <= m - 1 {
                    let mmj = m - j;
                    bsum += (j as f64 * r - mmj as f64) * a0[j - 1] * b0[mmj - 1];

                    j += 1;
                }
                b0[m - 1] = r * a0[m - 1] + bsum / m as f64;

                m += 1;
            }
            c[i - 1] = b0[i - 1] / (i as f64 + 1.0);

            let mut dsum = 0.0;
            let mut j = 1;
            while j <= i - 1 {
                dsum += d[i - j - 1] * c[j - 1];

                j += 1;
            }
            d[i - 1] = -(dsum + c[i - 1]);

            i += 1;
        }

        j0 = e1 * znm1 + (n as f64 - 1.0) * j0;
        j1 = e1 * zn + n as f64 * j1;
        znm1 = z2 * znm1;
        zn = z2 * zn;
        w *= w0;
        let t0 = d[n - 1] * w * j0;
        w *= w0;
        let t1 = d[np1 - 1] * w * j1;
        sum += t0 + t1;
        if t0.abs() + t1.abs() <= eps * sum {
            break;
        }
        n += 2;
    }

    if log {
        return ln_e0 + t - bcorr(a, b) + sum.ln();
    } else {
        let u = (-bcorr(a, b)).exp();
        return e0 * t * u * sum;
    }
} /* basym_ */

///     If l = 0 then  exparg(l) = The largest positive W for which
///     exp(W) can be computed. With 0.99999 fuzz  ==> exparg(0) =   709.7756  nowadays
///
///     if l = 1 (nonzero) then  exparg(l) = the largest negative W for
///     which the computed value of exp(W) is nonzero.
///     With 0.99999 fuzz			  ==> exparg(1) =  -709.0825  nowadays
///
///     Note... only an approximate value for exparg(L) is needed.
fn exparg(l: usize) -> f64 {
    let lnb = 0.69314718055995;
    let m = if l == 0 {
        f64::MAX_EXP
    } else {
        f64::MIN_EXP - 1
    };

    return m as f64 * lnb * 0.99999;
} /* exparg */

/// EVALUATION OF EXP(MU + X) */
fn esum(mu: isize, x: f64, give_log: bool) -> f64 {
    if give_log {
        return x + mu as f64;
    }

    // else :
    let w;
    if x > 0.0 {
        /* L10: */
        if mu > 0 {
            return (mu as f64).exp() * x.exp();
        }
        w = mu as f64 + x;
        if w < 0.0 {
            return (mu as f64).exp() * x.exp();
        }
    } else {
        /* x <= 0 */
        if mu < 0 {
            return (mu as f64).exp() * x.exp();
        }
        w = mu as f64 + x;
        if w > 0.0 {
            return (mu as f64).exp() * x.exp();
        }
    }
    return w.exp();
} /* esum */

/// EVALUATION OF THE FUNCTION EXP(X) - 1
fn rexpm1(x: f64) -> f64 {
    let p1 = 9.14041914819518e-10;
    let p2 = 0.0238082361044469;
    let q1 = -0.499999999085958;
    let q2 = 0.107141568980644;
    let q3 = -0.0119041179760821;
    let q4 = 5.95130811860248e-4;

    if x.abs() <= 0.15 {
        return x * (((p2 * x + p1) * x + 1.0) / ((((q4 * x + q3) * x + q2) * x + q1) * x + 1.0));
    } else {
        /* |x| > 0.15 : */
        let w = x.exp();
        return if x > 0.0 {
            w * (0.5 - 1.0 / w + 0.5)
        } else {
            w - 0.5 - 0.5
        };
    }
} /* rexpm1 */

/// Evaluation of the function ln(1 + a)
fn alnrel(a: f64) -> f64 {
    if a.abs() > 0.375 {
        return (1.0 + a).ln();
    }
    // else : |a| <= 0.375
    let p1 = -1.29418923021993;
    let p2 = 0.405303492862024;
    let p3 = -0.0178874546012214;
    let q1 = -1.62752256355323;
    let q2 = 0.747811014037616;
    let q3 = -0.0845104217945565;

    let t = a / (a + 2.0);
    let t2 = t * t;
    let w = (((p3 * t2 + p2) * t2 + p1) * t2 + 1.0) / (((q3 * t2 + q2) * t2 + q1) * t2 + 1.0);
    return t * 2.0 * w;
} /* alnrel */

/// Evaluation of the function  x - ln(1 + x)
fn rlog1(x: f64) -> f64 {
    let a = 0.0566749439387324;
    let b = 0.0456512608815524;
    let p0 = 0.333333333333333;
    let p1 = -0.224696413112536;
    let p2 = 0.00620886815375787;
    let q1 = -1.27408923933623;
    let q2 = 0.354508718369557;

    let mut h;
    let r;
    let t;
    let w;
    let w1;
    if x < -0.39 || x > 0.57 {
        /* direct evaluation */
        w = x + 0.5 + 0.5;
        return x - w.ln();
    }
    /* else */
    if x < -0.18 {
        /* L10: */
        h = x + 0.3;
        h /= 0.7;
        w1 = a - h * 0.3;
    } else if x > 0.18 {
        /* L20: */
        h = x * 0.75 - 0.25;
        w1 = b + h / 3.0;
    } else {
        /*		Argument Reduction */
        h = x;
        w1 = 0.0;
    }

    /* L30:              	Series Expansion */

    r = h / (h + 2.0);
    t = r * r;
    w = ((p2 * t + p1) * t + p0) / ((q2 * t + q1) * t + 1.0);
    return t * 2.0 * (1.0 / (1.0 - r) - r * w) + w1;
} /* rlog1 */

/// EVALUATION OF THE REAL ERROR FUNCTION
fn erf__(x: f64) -> f64 {
    /* Initialized data */

    let c = 0.564189583547756;
    let a = [
        7.7105849500132e-5,
        -0.00133733772997339,
        0.0323076579225834,
        0.0479137145607681,
        0.128379167095513,
    ];
    let b = [0.00301048631703895, 0.0538971687740286, 0.375795757275549];
    let p = [
        -1.36864857382717e-7,
        0.564195517478974,
        7.21175825088309,
        43.1622272220567,
        152.98928504694,
        339.320816734344,
        451.918953711873,
        300.459261020162,
    ];
    let q = [
        1.,
        12.7827273196294,
        77.0001529352295,
        277.585444743988,
        638.980264465631,
        931.35409485061,
        790.950925327898,
        300.459260956983,
    ];
    let r = [
        2.10144126479064,
        26.2370141675169,
        21.3688200555087,
        4.6580782871847,
        0.282094791773523,
    ];
    let s = [
        94.153775055546,
        187.11481179959,
        99.0191814623914,
        18.0124575948747,
    ];

    /* Local variables */
    let mut t;
    let x2;
    let ax;
    let bot;
    let top;

    ax = x.abs();
    if ax <= 0.5 {
        t = x * x;
        top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1.0;
        bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.0;

        return x * (top / bot);
    }

    // else:  |x| > 0.5

    if ax <= 4.0 {
        //  |x| in (0.5, 4]
        top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) * ax + p[5]) * ax
            + p[6])
            * ax
            + p[7];
        bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) * ax + q[5]) * ax
            + q[6])
            * ax
            + q[7];
        let R = 0.5 - (-x * x).exp() * top / bot + 0.5;
        return if x < 0.0 { -R } else { R };
    }

    // else:  |x| > 4

    if ax >= 5.8 {
        return if x > 0.0 { 1.0 } else { -1.0 };
    }

    // else:  4 < |x| < 5.8
    x2 = x * x;
    t = 1.0 / x2;
    top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4];
    bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.0;
    t = (c - top / (x2 * bot)) / ax;
    let R = 0.5 - (-x2).exp() * t + 0.5;
    return if x < 0.0 { -R } else { R };
} /* erf */

/// EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
///
/// ERFC1(IND,X) = ERFC(X)            IF IND = 0
/// ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
fn erfc1(ind: usize, x: f64) -> f64 {
    /* Initialized data */

    let c = 0.564189583547756;
    let a = [
        7.7105849500132e-5,
        -0.00133733772997339,
        0.0323076579225834,
        0.0479137145607681,
        0.128379167095513,
    ];
    let b = [0.00301048631703895, 0.0538971687740286, 0.375795757275549];
    let p = [
        -1.36864857382717e-7,
        0.564195517478974,
        7.21175825088309,
        43.1622272220567,
        152.98928504694,
        339.320816734344,
        451.918953711873,
        300.459261020162,
    ];
    let q = [
        1.,
        12.7827273196294,
        77.0001529352295,
        277.585444743988,
        638.980264465631,
        931.35409485061,
        790.950925327898,
        300.459260956983,
    ];
    let r = [
        2.10144126479064,
        26.2370141675169,
        21.3688200555087,
        4.6580782871847,
        0.282094791773523,
    ];
    let s = [
        94.153775055546,
        187.11481179959,
        99.0191814623914,
        18.0124575948747,
    ];

    let mut ret_val;
    let e;
    let mut t;
    let w;
    let bot;
    let top;

    let ax = x.abs();
    //				|X| <= 0.5 */
    if ax <= 0.5 {
        let t = x * x;
        let top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1.;
        let bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.;
        ret_val = 0.5 - x * (top / bot) + 0.5;
        if ind != 0 {
            ret_val = t.exp() * ret_val;
        }
        return ret_val;
    }
    // else (L10:):		0.5 < |X| <= 4
    if ax <= 4. {
        top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) * ax + p[5]) * ax
            + p[6])
            * ax
            + p[7];
        bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) * ax + q[5]) * ax
            + q[6])
            * ax
            + q[7];
        ret_val = top / bot;
    } else {
        //			|X| > 4
        // L20:
        if x <= -5.6 {
            // L50:            	LIMIT VALUE FOR "LARGE" NEGATIVE X
            ret_val = 2.0;
            if ind != 0 {
                ret_val = (x * x).exp() * 2.0;
            }
            return ret_val;
        }
        if ind == 0 && (x > 100.0 || x * x > -exparg(1)) {
            // LIMIT VALUE FOR LARGE POSITIVE X   WHEN IND = 0
            // L60:
            return 0.0;
        }

        // L30:
        t = 1.0 / (x * x);
        top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4];
        bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.0;
        ret_val = (c - t * top / bot) / ax;
    }

    // L40:                 FINAL ASSEMBLY
    if ind != 0 {
        if x < 0.0 {
            ret_val = (x * x).exp() * 2.0 - ret_val;
        }
    } else {
        // L41:  ind == 0 :
        w = x * x;
        t = w;
        e = w - t;
        ret_val = (0.5 - e + 0.5) * (-t).exp() * ret_val;
        if x < 0.0 {
            ret_val = 2.0 - ret_val;
        }
    }
    return ret_val;
} /* erfc1 */

/// COMPUTATION OF 1/GAMMA(A+1) - 1  FOR -0.5 <= A <= 1.5
fn gam1(a: f64) -> f64 {
    let d;
    let mut t;
    let w;
    let bot;
    let top;

    t = a;
    d = a - 0.5;
    // t := if(a > 1/2)  a-1  else  a
    if d > 0.0 {
        t = d - 0.5;
    }
    if t < 0.0 {
        /* L30: */
        let r = [
            -0.422784335098468,
            -0.771330383816272,
            -0.244757765222226,
            0.118378989872749,
            9.30357293360349e-4,
            -0.0118290993445146,
            0.00223047661158249,
            2.66505979058923e-4,
            -1.32674909766242e-4,
        ];
        let s1 = 0.273076135303957;
        let s2 = 0.0559398236957378;

        top = (((((((r[8] * t + r[7]) * t + r[6]) * t + r[5]) * t + r[4]) * t + r[3]) * t + r[2])
            * t
            + r[1])
            * t
            + r[0];
        bot = (s2 * t + s1) * t + 1.;
        w = top / bot;
        // R_ifDEBUG_printf("  gam1(a = %.15g): t < 0: w=%.15g\n", a, w);
        if d > 0.0 {
            return t * w / a;
        } else {
            return a * (w + 0.5 + 0.5);
        }
    } else if t == 0.0 {
        // L10: a in {0, 1}
        return 0.0;
    } else {
        /* t > 0;  L20: */
        let p = [
            0.577215664901533,
            -0.409078193005776,
            -0.230975380857675,
            0.0597275330452234,
            0.0076696818164949,
            -0.00514889771323592,
            5.89597428611429e-4,
        ];
        let q = [
            1.,
            0.427569613095214,
            0.158451672430138,
            0.0261132021441447,
            0.00423244297896961,
        ];

        top = (((((p[6] * t + p[5]) * t + p[4]) * t + p[3]) * t + p[2]) * t + p[1]) * t + p[0];
        bot = (((q[4] * t + q[3]) * t + q[2]) * t + q[1]) * t + 1.;
        w = top / bot;
        // R_ifDEBUG_printf("  gam1(a = %.15g): t > 0: (is a < 1.5 ?)  w=%.15g\n",
        // a, w);
        if d > 0.0 {
            /* L21: */
            return t / a * (w - 0.5 - 0.5);
        } else {
            return a * w;
        }
    }
} /* gam1 */

/// EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 <= A <= 1.25
fn gamln1(a: f64) -> f64 {
    let w;
    if a < 0.6 {
        let p0 = 0.577215664901533;
        let p1 = 0.844203922187225;
        let p2 = -0.168860593646662;
        let p3 = -0.780427615533591;
        let p4 = -0.402055799310489;
        let p5 = -0.0673562214325671;
        let p6 = -0.00271935708322958;
        let q1 = 2.88743195473681;
        let q2 = 3.12755088914843;
        let q3 = 1.56875193295039;
        let q4 = 0.361951990101499;
        let q5 = 0.0325038868253937;
        let q6 = 6.67465618796164e-4;
        w = ((((((p6 * a + p5) * a + p4) * a + p3) * a + p2) * a + p1) * a + p0)
            / ((((((q6 * a + q5) * a + q4) * a + q3) * a + q2) * a + q1) * a + 1.);
        return -(a) * w;
    } else {
        /* 0.6 <= a <= 1.25 */
        let r0 = 0.422784335098467;
        let r1 = 0.848044614534529;
        let r2 = 0.565221050691933;
        let r3 = 0.156513060486551;
        let r4 = 0.017050248402265;
        let r5 = 4.97958207639485e-4;
        let s1 = 1.24313399877507;
        let s2 = 0.548042109832463;
        let s3 = 0.10155218743983;
        let s4 = 0.00713309612391;
        let s5 = 1.16165475989616e-4;
        let x = a - 0.5 - 0.5;
        w = (((((r5 * x + r4) * x + r3) * x + r2) * x + r1) * x + r0)
            / (((((s5 * x + s4) * x + s3) * x + s2) * x + s1) * x + 1.);
        return x * w;
    }
} /* gamln1 */

/// Evaluation of the Digamma function psi(x)
///
/// -----------
///
/// Psi(xx) is assigned the value 0 when the digamma function cannot
/// be computed.
///
/// The main computation involves evaluation of rational Chebyshev
/// approximations published in Math. Comp. 27, 123-127(1973) by
/// Cody, Strecok and Thacher.
fn psi(mut x: f64) -> f64 {
    /* --------------------------------------------------------------------- */
    /*     Psi was written at Argonne National Laboratory for the FUNPACK */
    /*     package of special function subroutines. Psi was modified by */
    /*     A.H. Morris (NSWC). */
    /* --------------------------------------------------------------------- */

    let piov4 = 0.785398163397448; /* == pi / 4 */
    /*     dx0 = zero of psi() to extended precision : */
    let dx0 = 1.461632144968362341262659542325721325;

    /* --------------------------------------------------------------------- */
    /*     COEFFICIENTS FOR RATIONAL APPROXIMATION OF */
    /*     PSI(X) / (X - X0),  0.5 <= X <= 3. */
    let p1 = [
        0.0089538502298197,
        4.77762828042627,
        142.441585084029,
        1186.45200713425,
        3633.51846806499,
        4138.10161269013,
        1305.60269827897,
    ];
    let q1 = [
        44.8452573429826,
        520.752771467162,
        2210.0079924783,
        3641.27349079381,
        1908.310765963,
        6.91091682714533e-6,
    ];
    /* --------------------------------------------------------------------- */

    /* --------------------------------------------------------------------- */
    /*     COEFFICIENTS FOR RATIONAL APPROXIMATION OF */
    /*     PSI(X) - LN(X) + 1 / (2*X),  X > 3. */

    let p2 = [
        -2.12940445131011,
        -7.01677227766759,
        -4.48616543918019,
        -0.648157123766197,
    ];
    let q2 = [
        32.2703493791143,
        89.2920700481861,
        54.6117738103215,
        7.77788548522962,
    ];
    /* --------------------------------------------------------------------- */

    let mut i;
    let mut m;
    let mut n;
    let mut nq;
    let d2;
    let mut w;
    let z;
    let mut den;
    let mut aug;
    let mut sgn;
    let xmx0;
    let mut xmax1;
    let mut upper;
    let xsmall;

    /* --------------------------------------------------------------------- */

    /*     MACHINE DEPENDENT CONSTANTS ... */

    /* --------------------------------------------------------------------- */
    /*	  XMAX1	 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
              WITH ENTIRELY INT REPRESENTATION.  ALSO USED
              AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
              ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
              PSI MAY BE REPRESENTED AS LOG(X).
    * Originally:  xmax1 = amin1(ipmpar(3), 1./spmpar(1))  */
    xmax1 = i32::MAX;
    d2 = 0.5 / (0.5 * f64::EPSILON); /*= 0.5 / (0.5 * DBL_EPS) = 1/DBL_EPSILON = 2^52 */
    if xmax1 > d2 as i32 {
        xmax1 = d2 as i32
    };

    /* --------------------------------------------------------------------- */
    /*        XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) */
    /*                 MAY BE REPRESENTED BY 1/X. */
    xsmall = 1e-9;
    /* --------------------------------------------------------------------- */
    aug = 0.;
    if x < 0.5 {
        /* --------------------------------------------------------------------- */
        /*     X < 0.5,  USE REFLECTION FORMULA */
        /*     PSI(1-X) = PSI(X) + PI * COTAN(PI*X) */
        /* --------------------------------------------------------------------- */
        if x.abs() <= xsmall {
            if x == 0.0 {
                return 0.0;
            }
            /* --------------------------------------------------------------------- */
            /*     0 < |X| <= XSMALL.  USE 1/X AS A SUBSTITUTE */
            /*     FOR  PI*COTAN(PI*X) */
            /* --------------------------------------------------------------------- */
            aug = -1.0 / x;
        } else {
            /* |x| > xsmall */
            /* --------------------------------------------------------------------- */
            /*     REDUCTION OF ARGUMENT FOR COTAN */
            /* --------------------------------------------------------------------- */
            /* L100: */
            w = -x;
            sgn = piov4;
            if w <= 0.0 {
                w = -w;
                sgn = -sgn;
            }
            /* --------------------------------------------------------------------- */
            /*     MAKE AN ERROR EXIT IF |X| >= XMAX1 */
            /* --------------------------------------------------------------------- */
            if w >= xmax1 as f64 {
                return 0.0;
            }
            nq = w as isize;
            w -= nq as f64;
            nq = (w * 4.0) as isize;
            w = (w - nq as f64 * 0.25) * 4.0;
            /* --------------------------------------------------------------------- */
            /*     W IS NOW RELATED TO THE FRACTIONAL PART OF  4. * X. */
            /*     ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST */
            /*     QUADRANT AND DETERMINE SIGN */
            /* --------------------------------------------------------------------- */
            n = nq / 2;
            if n + n != nq {
                w = 1. - w;
            }
            z = piov4 * w;
            m = n / 2;
            if m + m != n {
                sgn = -sgn;
            }
            /* --------------------------------------------------------------------- */
            /*     DETERMINE FINAL VALUE FOR  -PI*COTAN(PI*X) */
            /* --------------------------------------------------------------------- */
            n = (nq + 1) / 2;
            m = n / 2;
            m += m;
            if m == n {
                /* --------------------------------------------------------------------- */
                /*     CHECK FOR SINGULARITY */
                /* --------------------------------------------------------------------- */
                if z == 0.0 {
                    return 0.0;
                }
                /* --------------------------------------------------------------------- */
                /*     USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND */
                /*     SIN/COS AS A SUBSTITUTE FOR TAN */
                /* --------------------------------------------------------------------- */
                aug = sgn * (z.cos() / z.sin() * 4.0);
            } else {
                /* L140: */
                aug = sgn * (z.sin() / z.cos() * 4.0);
            }
        }

        x = 1.0 - x;
    }
    /* L200: */
    if x <= 3.0 {
        /* --------------------------------------------------------------------- */
        /*     0.5 <= X <= 3. */
        /* --------------------------------------------------------------------- */
        den = x;
        upper = p1[0] * x;

        i = 1;
        while i <= 5 {
            den = (den + q1[i - 1]) * x;
            upper = (upper + p1[i]) * x;

            i += 1;
        }

        den = (upper + p1[6]) / (den + q1[5]);
        xmx0 = x - dx0;
        return den * xmx0 + aug;
    }

    /* --------------------------------------------------------------------- */
    /*     IF X >= XMAX1, PSI = LN(X) */
    /* --------------------------------------------------------------------- */
    if x < xmax1 as f64 {
        /* --------------------------------------------------------------------- */
        /*     3. < X < XMAX1 */
        /* --------------------------------------------------------------------- */
        w = 1.0 / (x * x);
        den = w;
        upper = p2[0] * w;

        i = 1;
        while i <= 3 {
            den = (den + q2[i - 1]) * w;
            upper = (upper + p2[i]) * w;

            i += 1;
        }

        aug = upper / (den + q2[3]) - 0.5 / x + aug;
    }
    return aug + x.ln();
} /* psi */

#[derive(Copy, Clone, Debug, Default)]
struct BetaLnState {
    a: f64,
    b: f64,
    w: f64,
    ret: f64,
    exit: bool,
}

/// Evaluation of the logarithm of the beta function  ln(beta(a0,b0))
fn betaln(a0: f64, b0: f64) -> f64 {
    let mut state = BetaLnState::default();
    let e = 0.918938533204673; /* e == 0.5*LN(2*PI) */

    state.a = a0.min(b0);
    state.b = a0.max(b0);

    if state.a < 8.0 {
        if state.a < 1.0 {
            /* ----------------------------------------------------------------------- */
            //                    		A < 1
            /* ----------------------------------------------------------------------- */
            if state.b < 8.0 {
                return gamln(state.a) + (gamln(state.b) - gamln(state.a + state.b));
            } else {
                return gamln(state.a) + algdiv(state.a, state.b);
            }
        }
        /* else */
        /* ----------------------------------------------------------------------- */
        //				1 <= A < 8
        /* ----------------------------------------------------------------------- */
        if state.a < 2.0 {
            if state.b <= 2.0 {
                return gamln(state.a) + gamln(state.b) - gsumln(state.a, state.b);
            }
            /* else */

            if state.b < 8.0 {
                state.w = 0.0;
                l40(&mut state);
                if state.exit {
                    return state.ret;
                } else {
                    unreachable!()
                }
            }
            return gamln(state.a) + algdiv(state.a, state.b);
        }
        // else L30:    REDUCTION OF A WHEN B <= 1000

        if state.b <= 1e3 {
            let n = state.a - 1.0;
            state.w = 1.0;
            let mut i = 1;
            while i as f64 <= n {
                state.a += -1.0;
                let h = state.a / state.b;
                state.w *= h / (h + 1.0);

                i += 1;
            }
            state.w = state.w.ln();

            if state.b >= 8.0 {
                return state.w + gamln(state.a) + algdiv(state.a, state.b);
            }

            // else
            l40(&mut state);
            if state.exit {
                return state.ret;
            } else {
                unreachable!()
            }
        } else {
            // L50:	reduction of A when  B > 1000
            let n = (state.a - 1.0) as usize;
            state.w = 1.0;
            let mut i = 1;
            while i <= n {
                state.a += -1.0;
                state.w *= state.a / (state.a / state.b + 1.0);

                i += 1;
            }
            return state.w.ln() - n as f64 * state.b.ln()
                + (gamln(state.a) + algdiv(state.a, state.b));
        }
    } else {
        /* ----------------------------------------------------------------------- */
        // L60:			A >= 8
        /* ----------------------------------------------------------------------- */

        let w = bcorr(state.a, state.b);
        let h = state.a / state.b;
        let u = -(state.a - 0.5) * (h / (h + 1.0)).ln();
        let v = state.b * alnrel(h);
        if u > v {
            return state.b.ln() * -0.5 + e + w - v - u;
        } else {
            return state.b.ln() * -0.5 + e + w - u - v;
        }
    }
} /* betaln */

#[inline]
fn l40(state: &mut BetaLnState) {
    // 	1 < A <= B < 8 :  reduction of B
    let n = (state.b - 1.0) as usize;
    let mut z = 1.0;
    let mut i = 1;
    while i <= n {
        state.b += -1.0;
        z *= state.b / (state.a + state.b);

        i += 1;
    }
    state.ret = state.w + z.ln() + (gamln(state.a) + (gamln(state.b) - gsumln(state.a, state.b)));
    state.exit = true;
}

/// EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
/// FOR 1 <= A <= 2  AND  1 <= B <= 2
fn gsumln(a: f64, b: f64) -> f64 {
    let x = a + b - 2.0; /* in [0, 2] */

    if x <= 0.25 {
        return gamln1(x + 1.0);
    }

    /* else */
    if x <= 1.25 {
        return gamln1(x) + alnrel(x);
    }
    /* else x > 1.25 : */
    return gamln1(x - 1.0) + (x * (x + 1.0)).ln();
} /* gsumln */

/// EVALUATION OF  DEL(A0) + DEL(B0) - DEL(A0 + B0)  WHERE
/// LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
/// IT IS ASSUMED THAT A0 >= 8 AND B0 >= 8.
fn bcorr(a0: f64, b0: f64) -> f64 {
    /* Initialized data */

    let c0 = 0.0833333333333333;
    let c1 = -0.00277777777760991;
    let c2 = 7.9365066682539e-4;
    let c3 = -5.9520293135187e-4;
    let c4 = 8.37308034031215e-4;
    let c5 = -0.00165322962780713;

    /* System generated locals */
    let ret_val;
    let mut r1;

    /* Local variables */
    let a;
    let b;
    let c;
    let h;
    let mut t;
    let mut w;
    let x;
    let s3;
    let s5;
    let x2;
    let s7;
    let s9;
    let s11;
    /* ------------------------ */
    a = a0.min(b0);
    b = a0.max(b0);

    h = a / b;
    c = h / (h + 1.0);
    x = 1.0 / (h + 1.0);
    x2 = x * x;

    /*                SET SN = (1 - X^N)/(1 - X) */

    s3 = x + x2 + 1.0;
    s5 = x + x2 * s3 + 1.0;
    s7 = x + x2 * s5 + 1.0;
    s9 = x + x2 * s7 + 1.0;
    s11 = x + x2 * s9 + 1.0;

    /*                SET W = DEL(B) - DEL(A + B) */

    /* Computing 2nd power */
    r1 = 1.0 / b;
    t = r1 * r1;
    w = ((((c5 * s11 * t + c4 * s9) * t + c3 * s7) * t + c2 * s5) * t + c1 * s3) * t + c0;
    w *= c / b;

    /*                   COMPUTE  DEL(A) + W */

    /* Computing 2nd power */
    r1 = 1.0 / a;
    t = r1 * r1;
    ret_val = (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a + w;
    return ret_val;
} /* bcorr */

/// COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B >= 8
///
/// --------
///
/// IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
/// LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
fn algdiv(a: f64, b: f64) -> f64 {
    /* Initialized data */

    let c0 = 0.0833333333333333;
    let c1 = -0.00277777777760991;
    let c2 = 7.9365066682539e-4;
    let c3 = -5.9520293135187e-4;
    let c4 = 8.37308034031215e-4;
    let c5 = -0.00165322962780713;

    let c;
    let d;
    let h;
    let t;
    let u;
    let v;
    let mut w;
    let x;
    let s3;
    let s5;
    let x2;
    let s7;
    let s9;
    let s11;

    /* ------------------------ */
    if a > b {
        h = b / a;
        c = 1.0 / (h + 1.0);
        x = h / (h + 1.0);
        d = a + (b - 0.5);
    } else {
        h = a / b;
        c = h / (h + 1.0);
        x = 1.0 / (h + 1.0);
        d = b + (a - 0.5);
    }

    /* Set s<n> = (1 - x^n)/(1 - x) : */

    x2 = x * x;
    s3 = x + x2 + 1.0;
    s5 = x + x2 * s3 + 1.0;
    s7 = x + x2 * s5 + 1.0;
    s9 = x + x2 * s7 + 1.0;
    s11 = x + x2 * s9 + 1.0;

    /* w := Del(b) - Del(a + b) */

    t = 1.0 / (b * b);
    w = ((((c5 * s11 * t + c4 * s9) * t + c3 * s7) * t + c2 * s5) * t + c1 * s3) * t + c0;
    w *= c / b;

    /*                    COMBINE THE RESULTS */

    u = d * alnrel(a / b);
    v = a * (b.ln() - 1.0);
    if u > v {
        return w - v - u;
    } else {
        return w - u - v;
    }
} /* algdiv */

/// Evaluation of  ln(gamma(a))  for positive a
/// -----------------------------------------------------------------------
/// Written by Alfred H. Morris
/// Naval Surface Warfare Center
/// Dahlgren, Virginia
fn gamln(a: f64) -> f64 {
    let d = 0.418938533204673; /* d == 0.5*(LN(2*PI) - 1) */

    let c0 = 0.0833333333333333;
    let c1 = -0.00277777777760991;
    let c2 = 7.9365066682539e-4;
    let c3 = -5.9520293135187e-4;
    let c4 = 8.37308034031215e-4;
    let c5 = -0.00165322962780713;

    if a <= 0.8 {
        return gamln1(a) - a.ln(); /* ln(G(a+1)) - ln(a) == ln(G(a+1)/a) = ln(G(a)) */
    } else if a <= 2.25 {
        return gamln1(a - 0.5 - 0.5);
    } else if a < 10.0 {
        let mut i;
        let n = (a - 1.25) as usize;
        let mut t = a;
        let mut w = 1.0;

        i = 1;
        while i <= n {
            t += -1.0;
            w *= t;

            i += 1;
        }
        return gamln1(t - 1.0) + w.ln();
    } else {
        /* a >= 10 */
        let t = 1.0 / (a * a);
        let w = (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a;
        return d + w + (a - 0.5) * (a.ln() - 1.0);
    }
} /* gamln */