r2rs-nmath 0.1.1

// Translation of nmath's qbeta
//
// 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, FloatConst};
use strafe_type::{FloatConstraint, LogProbability64, Positive64, Probability64, Real64};

use crate::{distribution::beta::pbeta_raw, func::lbeta, traits::DPQ};

const USE_LOG_X_CUTOFF: f64 = -5.0;
const N_NEWTON_FREE: i32 = 4;

/// The beta quantile function
pub fn qbeta<PR: Into<Probability64>, PO1: Into<Positive64>, PO2: Into<Positive64>>(
    alpha: PR,
    p: PO1,
    q: PO2,
    lower_tail: bool,
) -> Real64 {
    let alpha = alpha.into().unwrap();
    qbeta_inner(alpha, p, q, lower_tail, false)
}

/// The beta quantile function
pub fn log_qbeta<LP: Into<LogProbability64>, P1: Into<Positive64>, P2: Into<Positive64>>(
    alpha: LP,
    p: P1,
    q: P2,
    lower_tail: bool,
) -> Real64 {
    let alpha = alpha.into().unwrap();
    qbeta_inner(alpha, p, q, lower_tail, true)
}

pub fn qbeta_inner<P1: Into<Positive64>, P2: Into<Positive64>>(
    alpha: f64,
    p: P1,
    q: P2,
    lower_tail: bool,
    log: bool,
) -> Real64 {
    let p = p.into().unwrap();
    let q = q.into().unwrap();

    let mut qbet: [f64; 2] = [0.0; 2]; // = { qbeta(), 1 - qbeta() }
    qbeta_raw(
        alpha,
        p,
        q,
        lower_tail,
        log,
        None,
        USE_LOG_X_CUTOFF,
        N_NEWTON_FREE,
        &mut qbet,
    );
    qbet[0].into()
}

// CARE: assumes subnormal numbers, i.e., no underflow at DBL_MIN:
static DBL_very_MIN: f64 = 2.2250738585072014e-308 / 4.0;
static DBL_log_v_MIN: f64 = std::f64::consts::LN_2 * (f64::MIN_EXP as f64 - 2.0);
// Too extreme: inaccuracy in pbeta(); e.g for  qbeta(0.95, 1e-9, 20):
// -> in pbeta() --> bgrat(..... b*z == 0 underflow, hence inaccurate pbeta()
/* DBL_very_MIN  = 0x0.0000001p-1022, // = 2^-1050 = 2^(-1022 - 28) */
/* DBL_log_v_MIN = -1050. * M_LN2, // = log(DBL_very_MIN) ~= -727.8045 */
// the most extreme -- not ok, as pbeta() then behaves strangely,
// e.g., for  qbeta(0.95, 1e-8, 20):
/* DBL_very_MIN  = 0x0.0000000000001p-1022, // = 2^-1074 = 2^(-1022 -52) */
/* DBL_log_v_MIN = -1074. * M_LN2, // = log(DBL_very_MIN) */

static DBL_1__eps: f64 = 0.9999999999999999; // = 1 - 2^-53

/* set the exponent of acu to -2r-2 for r digits of accuracy */
/*---- NEW ---- -- still fails for p = 1e11, q=.5*/

const FPU: f64 = 3e-308;
/// acu_min:  Minimal value for accuracy 'acu' which will depend on (a,p);
/// acu_min >= fpu
const ACU_MIN: f64 = 1e-300;
const P_LO: f64 = FPU;
const P_HI: f64 = 1.0 - 2.22e-16;

const CONST1: f64 = 2.30753;
const CONST2: f64 = 0.27061;
const CONST3: f64 = 0.99229;
const CONST4: f64 = 0.04481;

#[derive(Copy, Clone, Default, Debug)]
pub struct QBetaState {
    swap_choose: bool,
    swap_tail: bool,
    log_: bool,
    give_log_q: bool,
    use_log_x: bool,
    warned: bool,
    add_N_step: bool,
    a: f64,
    la: f64,
    logbeta: f64,
    g: f64,
    h: f64,
    pp: f64,
    p_: f64,
    qq: f64,
    r: f64,
    s: f64,
    t: f64,
    w: f64,
    y: f64,
    u: f64,
    xinbta: f64,
    n_maybe_swaps: usize,
    u_n: f64,
    acu: f64,
    tx: f64,
    exit: bool,
}

/// Returns both qbeta() and its "mirror" 1-qbeta(). Useful notably when qbeta() ~= 1
///
/// * swap01 = {TRUE, NA, FALSE}: if NA, algorithm decides swap_tail
/// * log_q_cut = if == Inf: return log(qbeta(..)); otherwise,
/// if finite: the bound for switching to log(x)-scale; see use_log_x
/// * n_N = number of "unconstrained" Newton steps before switching to constrained
/// * qb\[0:1\] = { qbeta(), 1 - qbeta() }
fn qbeta_raw(
    alpha: f64,
    p: f64,
    q: f64,
    lower_tail: bool,
    log: bool,
    swap_01: Option<bool>,
    log_q_cut: f64,
    n_N: i32,
    qb: &mut [f64],
) {
    let mut state = QBetaState::default();
    state.swap_choose = swap_01.is_none();
    state.give_log_q = log_q_cut == f64::infinity();
    state.use_log_x = state.give_log_q; // or u < log_q_cut  below
    state.warned = false;
    state.add_N_step = true;
    state.y = -1.0;

    // Assuming p >= 0, q >= 0  here ...

    let return_q_0 = |qb: &mut [f64]| {
        if state.give_log_q {
            qb[0] = f64::neg_infinity();
            qb[1] = 0.0;
        } else {
            qb[0] = 0.0;
            qb[1] = 1.0;
        }
    };
    let return_q_1 = |qb: &mut [f64]| {
        if state.give_log_q {
            qb[0] = 0.0;
            qb[1] = f64::neg_infinity();
        } else {
            qb[0] = 1.0;
            qb[1] = 0.0;
        }
    };

    // Deal with boundary cases here:
    if alpha == f64::dt_0(lower_tail, log) {
        return_q_0(qb);
        return;
    }
    if alpha == f64::dt_1(lower_tail, log) {
        return_q_1(qb);
        return;
    }

    // check alpha {*before* transformation which may lose all accuracy}:
    if (log && alpha > 0.0) || (!log && (alpha < 0.0 || alpha > 1.0)) {
        // alpha is outside
        // R_ifDEBUG_printf("qbeta(alpha=%g, %g, %g, .., log=%d): %s%s\n",
        //                  alpha, p,q, log, "alpha not in ",
        //                  log ? "[-Inf, 0]" : "[0,1]");
        // ML_WARN_return_NAN :
        // ML_WARNING(ME_DOMAIN, "");
        qb[0] = f64::nan();
        qb[1] = f64::nan();
        return;
    }

    //  p==0, q==0, p = Inf, q = Inf  <==> treat as one- or two-point mass
    if p == 0.0 || q == 0.0 || !p.is_finite() || !q.is_finite() {
        // We know 0 < T(alpha) < 1 : pbeta() is constant and trivial in {0, 1/2, 1}
        // R_ifDEBUG_printf(
        //     "qbeta(%g, %g, %g, lower_t=%d, log=%d): (p,q)-boundary: trivial\n",
        //     alpha, p,q, lower_tail, log);
        let return_q_half = |qb: &mut [f64]| {
            if state.give_log_q {
                qb[0] = -f64::LN_2();
                qb[1] = -f64::LN_2();
            } else {
                qb[0] = 0.5;
                qb[1] = 0.5;
            }
        };
        if p == 0.0 && q == 0.0 {
            // point mass 1/2 at each of {0,1} :
            if alpha < f64::d_half(log) {
                return_q_0(qb);
            } else if alpha > f64::d_half(log) {
                return_q_1(qb);
            } else {
                // alpha == "1/2"
                return_q_half(qb);
            }
        } else if p == 0.0 || p / q == 0.0 {
            // point mass 1 at 0 - "flipped around"
            return_q_0(qb);
        } else if q == 0.0 || q / p == 0.0 {
            // point mass 1 at 0 - "flipped around"
            return_q_1(qb);
        } else {
            // p = q = Inf : point mass 1 at 1/2
            return_q_half(qb);
        }
        return;
    }

    /* initialize */
    state.p_ = alpha.dt_qiv(lower_tail, log); /* lower_tail prob (in any case) */
    // Conceptually,  0 < p_ < 1  (but can be 0 or 1 because of cancellation!)
    state.logbeta = lbeta(p, q).unwrap();

    state.swap_tail = if state.swap_choose {
        state.p_ > 0.5
    } else {
        swap_01.unwrap()
    };

    // R_ifDEBUG_printf(
    //     "qbeta(%g, %g, %g, lower_t=%d, log=%d, swap_01=%d, log_q_cut=%g, n_N=%d):%s\n"
    //     "  swap_tail=%s :",
    //     alpha, p,q, lower_tail, log, swap_01, log_q_cut, n_N,
    //     (log && (p_ == 0. || p_ == 1.)) ? (p_==0.?" p_=0":" p_=1") : "",
    //     (swap_tail ? "TRUE": "F"));

    maybe_swap(
        &mut state, alpha, p, q, lower_tail, log, swap_01, log_q_cut, n_N, qb,
    );
    if state.exit {
        return;
    }

    l_newton(
        &mut state, alpha, p, q, lower_tail, log, swap_01, log_q_cut, n_N, qb,
    );
    if state.exit {
        return;
    }

    l_converged(&mut state, log);
    if state.exit {
        return;
    }

    l_return(&mut state, log, qb);
}

fn maybe_swap(
    state: &mut QBetaState,
    alpha: f64,
    p: f64,
    q: f64,
    lower_tail: bool,
    log: bool,
    swap_01: Option<bool>,
    log_q_cut: f64,
    n_N: i32,
    qb: &mut [f64],
) {
    // change tail; default (swap_01 = NA): afterwards 0 < a <= 1/2
    if state.swap_tail {
        /* change tail, swap copies of {p,q}:  p <-> q :*/
        state.a = alpha.dt_civ(lower_tail, log); // = 1 - p_ , is < 1/2 if(swap_choose)
                                                 /* la := log(a), but without numerical cancellation: */
        state.la = alpha.dt_clog(lower_tail, log);
        state.pp = q;
        state.qq = p;
    } else {
        state.a = state.p_;
        state.la = alpha.dt_log(lower_tail, log);
        state.pp = p;
        state.qq = q;
    }
    state.n_maybe_swaps += 1;

    /* calculate the initial approximation */

    /* Desired accuracy for Newton iterations (below) should depend on  (a,p)
    * This is from Remark .. on AS 109, adapted.
    * However, it's not clear if this is "optimal" for IEEE double prec.

    * acu = fmax2(acu_min, pow(10., -25. - 5./(pp * pp) - 1./(a * a)));

    * NEW: 'acu' accuracy NOT for squared adjustment, but simple;
    * ---- i.e.,  "new acu" = sqrt(old acu)
    */
    state.acu =
        ACU_MIN.max(10.0.powf(-13.0 - 2.5 / (state.pp * state.pp) - 0.5 / (state.a * state.a)));
    // try to catch  "extreme left tail" early
    let u0 = (state.la + state.pp.ln() + state.logbeta) / state.pp; // = log(x_0)
    let mut rp = state.pp * (1.0 - state.qq) / (state.pp + 1.0);
    let log_eps_c = f64::LN_2() * (1.0 - f64::MANTISSA_DIGITS as f64); // = log(DBL_EPSILON) = -36.04..

    state.t = 0.2;
    // FIXME: Factor 0.2 is a bit arbitrary;  '1' is clearly much too much.
    let u0_maybe = (f64::LN_2() * f64::MIN_EXP as f64) < u0 && u0 < -0.01;
    /* 1. cannot allow exp(u0) = 0 ==> exp(u1) = exp(u0) = 0
     * 2. must: u0 < 0, but too close to 0 <==> x = exp(u0) = 0.99.. */
    // R_ifDEBUG_printf(
    //     "  n_maybe_swaps=%d, la=%#8g, u0=%#8g -> u0_maybe=%s (bnd: %g (%g)); ",
    //     n_maybe_swaps, la, u0, (u0_maybe ? "TRUE" : "F"),
    // u0_maybe ? (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2. : -0.,
    // u0_maybe ?  t*log_eps_c - log(fabs(rp)) : -0.
    // );

    state.u_n = 1.0; // to be  log(xinbta) <==> xinbta = exp(u_n).  1 is impossible
    if u0_maybe &&
        // qq <= 2 && // <--- "arbitrary"
        // u0 <  t*log_eps_c - log(fabs(rp)) &&
        u0 < (state.t*log_eps_c - ((state.pp*(1.0-state.qq)*(2.0-state.qq)/(2.*(state.pp+2.0))).abs()).ln())/2.0
    {
        // TODO: maybe jump here from below, when initial u "fails" ?
        // L_tail_u:
        // MM's one-step correction (cheaper than 1 Newton!)
        rp = rp * u0.exp(); // = rp*x0
        if rp > -1.0 {
            state.u = u0 - rp.ln_1p() / state.pp;
            // R_ifDEBUG_printf("u1-u0=%9.3g --> choosing u = u1 = %g\n", u-u0, u);
        } else {
            state.u = u0;
            // R_ifDEBUG_printf("cannot cheaply improve u0\n");
        }
        state.tx = state.u.exp();
        state.xinbta = state.u.exp();
        state.use_log_x = true; // or (u < log_q_cut)  ??
        l_newton(
            state, alpha, p, q, lower_tail, log, swap_01, log_q_cut, n_N, qb,
        );
        if state.exit {
            return;
        }
    }

    // y := y_\alpha in AS 64 := Hastings(1955) approximation of qnorm(1 - a) :
    state.r = (-2.0 * state.la).sqrt();
    state.y = state.r - (CONST1 + CONST2 * state.r) / (1. + (CONST3 + CONST4 * state.r) * state.r);
    // R_ifDEBUG_printf("y_a(AS64)=%g\n  ", y);
    if state.pp > 1.0 && state.qq > 1.0 {
        // use  Carter(1947), see AS 109, remark '5.'
        state.r = (state.y * state.y - 3.) / 6.;
        state.s = 1.0 / (state.pp + state.pp - 1.0);
        state.t = 1.0 / (state.qq + state.qq - 1.0);
        state.h = 2.0 / (state.s + state.t);
        state.w = state.y * (state.h + state.r).sqrt() / state.h
            - (state.t - state.s) * (state.r + 5.0 / 6.0 - 2.0 / (3.0 * state.h));
        // R_ifDEBUG_printf("p,q > 1 => w=%g <= 300 ? ", w);
        if state.w > 300.0 {
            // exp(w+w) is huge or overflows
            state.t = state.w + state.w + state.qq.ln() - state.pp.ln(); // = argument of log1pexp(.)
            state.u = // log(xinbta) = - log1p(qq/pp * exp(w+w)) = -log(1 + exp(t))
                if state.t <= 18.0 {-state.t.exp().ln_1p()} else {-state.t - (-state.t).exp()};
            state.xinbta = state.u.exp();
        } else {
            state.xinbta = state.pp / (state.pp + state.qq * (state.w + state.w).exp());
            state.u = // log(xinbta)
                - (state.qq/state.pp * (state.w+state.w).exp()).ln_1p();
        }
    } else {
        // use the original AS 64 proposal, Scheffé-Tukey (1944) and Wilson-Hilferty
        state.r = state.qq + state.qq;
        /* A slightly more stable version of  t := \chi^2_{alpha} of AS 64
         * t = 1. / (9. * qq); t = r * R_pow_di(1. - t + y * sqrt(t), 3);  */
        state.t = 1.0 / (3.0 * state.qq.sqrt()); // = sqrt(t) of formula above
        state.t = state.r * (1.0 + state.t * (-state.t + state.y)).pow_di(3); // = \chi^2_{alpha} of AS 64
        state.s = 4.0 * state.pp + state.r - 2.0; // 4p + 2q - 2 = numerator of new t' = s / t = s / chi^2
                                                  // R_ifDEBUG_printf("min(p,q) <= 1: t=%g, s=%g", t, s);
        if state.t == 0.0 || (state.t < 0.0 && state.s >= state.t) {
            // cannot use chisq approx
            // x0 = 1 - { (1-a)*q*B(p,q) } ^{1/q}    {AS 65}
            // xinbta = 1. - exp((log(1-a)+ log(qq) + logbeta) / qq);
            let l1ma = /* := log(1-a), directly from alpha (as 'la' above);
			   though only seen very small improvements */
                if state.swap_tail {alpha.dt_log(lower_tail, log) } else { alpha.dt_clog(lower_tail, log)};

            // R_ifDEBUG_printf(
            //     " t <= 0  => AS65: log1p(-a)=%7g, better l1ma=%.15g, relD=%g\n",
            //     log1p(-a), l1ma, log1p(-a)/l1ma - 1);

            let xx = (l1ma + state.qq.ln() + state.logbeta) / state.qq;
            // R_ifDEBUG_printf("  xx = (l1ma + log(qq) + logbeta) / qq = %.10g; ", xx);
            if xx <= 0.0 {
                state.xinbta = -xx.exp_m1();
                state.u = xx.log1_exp(); // =  log(xinbta) = log(1 - exp(...A...))
                                         // R_ifDEBUG_printf(" xx <= 0 useful ==> u, xinbta\n");
            } else {
                // xx > 0 ==> 1 - e^xx < 0 .. is nonsense
                // R_ifDEBUG_printf(" xx > 0: xinbta:= 1-e^xx < 0 'nonsense'; ");
                // Try MM's one-step correction (or else u0)
                let r_ = rp * u0.exp();
                if r_ > -1.0 {
                    state.u = u0 - r_.ln_1p() / state.pp; // R_ifDEBUG_printf("u:= u0 - log1p(r_)/pp = %g\n");
                } else {
                    state.u = u0; // R_ifDEBUG_printf("u:= u0 = %g\n");
                }
                state.xinbta = state.u.exp();
            }
        } else {
            state.t = state.s / state.t;
            // R_ifDEBUG_printf(" t > 0 or s < t < 0:  new t := s/t = %g ( > 1 ?)\n", t);
            if state.t <= 1.0 {
                // cannot use chisq, either
                state.u = u0;
                state.xinbta = state.u.exp();
            } else {
                // (1+x0)/(1-x0) = t,  solved for x0 :
                state.xinbta = 1.0 - 2.0 / (state.t + 1.0);
                state.u = (-2.0 / (state.t + 1.0)).ln_1p();
            }
        }
    }

    // Problem: If initial u is completely wrong, we make a wrong decision here
    if state.swap_choose &&  //   vvvv/ why -exp(*)? u  on log-x scale! Swapping can be very good, but needs smart if(..)
        (( state.swap_tail && state.u >= -log_q_cut.exp()) || // ==> "swap back"
            (!state.swap_tail && state.u >= -(4.0*log_q_cut).exp() && state.pp / state.qq < 1000.0) // ==> "swap now"
        )
    {
        // if(swap_tail)
        // R_ifDEBUG_printf("\"swap back\" as u = %g >= -exp(log_q_cut);", u);
        // else
        // R_ifDEBUG_printf("\"swap now\" as u = %g >= -exp(4*log_q_cut) && pp/qq = %g < 1000;",
        //                  u, pp/qq);
        // reverse swap (and typically use_log_x)
        state.swap_tail = !state.swap_tail;

        if state.swap_tail {
            // "swap now" (much less easily)
            state.a = alpha.dt_civ(lower_tail, log); // needed ?
            state.la = alpha.dt_clog(lower_tail, log);
            state.pp = q;
            state.qq = p;
        } else {
            // "swap back" :
            state.a = state.p_;
            state.la = alpha.dt_log(lower_tail, log);
            state.pp = p;
            state.qq = q;
        }
        // R_ifDEBUG_printf(" ==> la = %g\n", la);
        // we could redo computations above, but this should be stable
        state.u = state.u.log1_exp();
        state.xinbta = state.u.exp();

        /* Careful: "swap now"  should not fail if
           1) the above initial xinbta is "completely wrong"
           2) The correction step can go outside (u_n > 0 ==>  e^u > 1 is illegal)
           e.g., for  qbeta(0.2066, 0.143891, 0.05)
        */
    }

    if !state.use_log_x {
        state.use_log_x = state.u < log_q_cut; // <==> xinbta = e^u < exp(log_q_cut)
    }
    let bad_u = !state.u.is_finite();
    let bad_init = bad_u || state.xinbta > P_HI;

    // R_ifDEBUG_printf(" -> u = %g, e^u = xinbta = %.16g, (Newton acu=%g%s%s%s)\n",
    //                  u, xinbta, acu, (bad_u ? ", ** bad u **" : ""),
    // ((bad_init && !bad_u) ? ", ** bad_init **" : ""),
    // (use_log_x ? ", on u = LOG(x) SCALE" : ""));

    state.tx = state.xinbta; // keeping "original initial x" (for now)

    if bad_u || state.u < log_q_cut {
        /* e.g.
        qbeta(0.21, .001, 0.05)
        try "left border" quickly, i.e.,
        try at smallest positive number: */
        state.w = pbeta_raw(DBL_very_MIN, state.pp, state.qq, true, log);
        if state.w > if log { state.la } else { state.a } {
            // R_ifDEBUG_printf(
            // " quantile is left of %g: boundary \"convergence\"\n", DBL_very_MIN);
            if log || (state.w - state.a).abs() < (0.0 - state.a).abs() {
                // DBL_very_MIN is better than 0
                state.tx = DBL_very_MIN;
                state.u_n = DBL_log_v_MIN; // = log(DBL_very_MIN)
            } else {
                state.tx = 0.0;
                state.u_n = f64::neg_infinity();
            }
            state.use_log_x = log;
            state.add_N_step = false;
            l_return(state, log, qb);
            state.exit = true;
            return;
        } else {
            // R_ifDEBUG_printf(" pbeta(%g, %g, %g, T, log) = %g <= %g (= %s) --> continuing\n",
            // DBL_very_MIN, pp,qq, w,
            // (log ? la : a), (log ? "la" : "a"));
            if state.u < DBL_log_v_MIN {
                state.u = DBL_log_v_MIN; // = log(DBL_very_MIN)
                state.xinbta = DBL_very_MIN;
            }
        }
    }

    /* Sometimes the approximation is negative (and == 0 is also not "ok") */
    if bad_init && !(state.use_log_x && state.tx > 0.0) {
        if state.u == f64::neg_infinity() {
            // R_ifDEBUG_printf("  u = -Inf;");
            state.u = f64::LN_2() * f64::MIN_EXP as f64;
            state.xinbta = f64::min_value();
        } else {
            // R_ifDEBUG_printf(" bad_init: u=%g, xinbta=%g;", u,xinbta);
            state.xinbta = if state.xinbta > 1.1 {
                // i.e. "way off"
                0.5 // otherwise, keep the respective boundary:
            } else if state.xinbta < P_LO {
                state.u.exp()
            } else {
                P_HI
            };
            if bad_u {
                state.u = state.xinbta.ln();
            }
            // otherwise: not changing "potentially better" u than the above
        }
        // R_ifDEBUG_printf(" -> (partly)new u=%g, xinbta=%g\n", u,xinbta);
    }
}

/// Solve for x by a modified Newton-Raphson method, using pbeta_raw()
fn l_newton(
    state: &mut QBetaState,
    alpha: f64,
    p: f64,
    q: f64,
    lower_tail: bool,
    log: bool,
    swap_01: Option<bool>,
    log_q_cut: f64,
    n_N: i32,
    qb: &mut [f64],
) {
    state.r = 1.0 - state.pp;
    state.t = 1.0 - state.qq;
    let mut wprev = 0.0;
    let mut prev = 1.0;
    let mut adj = 1.0; // -Wall

    if state.use_log_x {
        // find  log(xinbta) -- work in  u := log(x) scale
        // if(bad_init && tx > 0) xinbta = tx;// may have been better

        let mut i_pb = 0;
        while i_pb < 1000 {
            // using log == TRUE  unconditionally here
            /* FIXME: if exp(u) = xinbta underflows to 0,
             *  want different formula pbeta_log(u, ..) */
            state.y = pbeta_raw(
                state.xinbta,
                state.pp,
                state.qq,
                /*lower_tail = */ true,
                true,
            );

            /* w := Newton step size for   L(u) = log F(e^u)  =!= 0;   u := log(x)
             *   =  (L(.) - la) / L'(.);  L'(u)= (F'(e^u) * e^u ) / F(e^u)
             *   =  (L(.) - la)*F(.) / {F'(e^u) * e^u } =
             *   =  (L(.) - la) * e^L(.) * e^{-log F'(e^u) - u}
             *   =  ( y   - la) * e^{ y - u -log F'(e^u)}
             and  -log F'(x)= -log f(x) = - -logbeta + (1-p) log(x) + (1-q) log(1-x)
                                        = logbeta + (1-p) u + (1-q) log(1-e^u)
            */
            state.w = if state.y == f64::neg_infinity() {
                // y = -Inf  well possible: we are on log scale!
                0.0
            } else {
                (state.y - state.la)
                    * (state.y - state.u
                        + state.logbeta
                        + state.r * state.u
                        + state.t * state.u.log1_exp())
                    .exp()
            };
            // R_ifDEBUG_printf("N(i=%2d): u=%#20.16g, lnpb(e^u)=%#15.11g, w=%#12.7g",
            //                  i_pb, u, y, w);
            if !state.w.is_finite() {
                // what we should do is ==> go back, get better starting value
                // R_ifDEBUG_printf(" -- not finite --> %s\n",
                //                  (n_maybe_swaps <= 1) ? "goto maybe_swap" : "give up");
                if state.n_maybe_swaps <= 1 {
                    maybe_swap(
                        state, alpha, p, q, lower_tail, log, swap_01, log_q_cut, n_N, qb,
                    );
                    if state.exit {
                        return;
                    }
                }
                /* else  was 'break;' ...
                but rather give up returning NaN directly as in "normal scale" Newton */
                // ML_WARNING(ME_DOMAIN, "");
                qb[0] = f64::nan();
                qb[1] = f64::nan();
                state.exit = true;
                return;
            }
            if i_pb >= n_N && state.w * wprev <= 0.0 {
                prev = adj.abs().max(FPU);
            }
            // R_ifDEBUG_printf(", %s prev=%g,",
            //                  (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev);
            state.g = 1.0;
            let mut i_inn = 0;
            while i_inn < 1000 {
                adj = state.g * state.w;
                // safe guard (here, from the very beginning)
                if adj.abs() < prev {
                    state.u_n = state.u - adj; // u_{n+1} = u_n - g*w
                    if state.u_n <= 0.0 {
                        // <==> 0 <  xinbta := e^u  <= 1
                        if prev <= state.acu || state.w.abs() <= state.acu {
                            // R_ifDEBUG_printf(
                            //     " it{in}=%d, -adj=%g, %s <= acu  ==> convergence\n",
                            //     i_inn, -adj, (prev <= acu) ? "prev" : "|w|");
                            l_converged(state, log);
                            if state.exit {
                                return;
                            }
                        }
                        // if (u_n != ML_NEGINF && u_n != 1)
                        break;
                    }
                }
                state.g /= 3.0;
                i_inn += 1;
            }
            // (cancellation in (u_n -u) => may differ from adj:
            let D = adj.abs().min((state.u_n - state.u).abs());
            /* R_ifDEBUG_printf(" delta(u)=%g\n", u_n - u); */
            // R_ifDEBUG_printf(" it{in}=%d, d.(u)=%6.3g, D/|.|=%.3g\n",
            //                  i_inn, u_n - u, D/fabs(u_n + u));
            if D <= 4e-16 * (state.u_n + state.u).abs() {
                l_converged(state, log);
                if state.exit {
                    return;
                }
            }
            state.u = state.u_n;
            state.xinbta = state.u.exp();
            wprev = state.w;
            i_pb += 1;
        } // for(i )
    } else {
        // "normal scale" Newton

        let mut i_pb = 0;
        while i_pb < 1000 {
            state.y = pbeta_raw(
                state.xinbta,
                state.pp,
                state.qq,
                /*lower_tail = */ true,
                log,
            );
            // delta{y} :   d_y = y - (log ? la : a);

            /* w := Newton step size  (F(.) - a) / F'(.)  or,
             * --   log: (lF - la) / (F' / F) = exp(lF) * (lF - la) / F'
             */
            state.w = if log {
                (state.y - state.la)
                    * (state.y
                        + state.logbeta
                        + state.r * state.xinbta.ln()
                        + state.t * (-state.xinbta).ln_1p())
                    .exp()
            } else {
                (state.y - state.a)
                    * (state.logbeta
                        + state.r * state.xinbta.ln()
                        + state.t * (-state.xinbta).ln_1p())
                    .exp()
            };
            if !state.w.is_finite() {
                // what we should do is ==> go back, get better starting value
                // R_ifDEBUG_printf(
                //     "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g -- is not finite --> %s\n",
                //     i_pb, xinbta, y, w, (n_maybe_swaps <= 2) ? "goto maybe_swap" : "give up");
                if state.n_maybe_swaps <= 2 {
                    if !log && state.n_maybe_swaps == 2 {
                        state.use_log_x = true;
                    } // try now
                    if !log || state.n_maybe_swaps <= 1 {
                        maybe_swap(
                            state, alpha, p, q, lower_tail, log, swap_01, log_q_cut, n_N, qb,
                        );
                        if state.exit {
                            return;
                        }
                    }
                }
                /* else  was 'break;' ...
                but rather give up returning NaN directly as in "normal scale" Newton */
                // ML_WARNING(ME_DOMAIN, "");
                qb[0] = f64::nan();
                qb[1] = f64::nan();
                state.exit = true;
                return;
            }
            if i_pb >= n_N && state.w * wprev <= 0.0 {
                prev = adj.abs().max(FPU);
            }
            // R_ifDEBUG_printf(
            //     "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g, %s prev=%g,",
            //     i_pb, xinbta, y, w,
            //     (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev);
            state.g = 1.0;
            let mut i_inn = 0;
            while i_inn < 1000 {
                adj = state.g * state.w;
                // take full Newton steps at the beginning; only then safe guard:
                if i_pb < n_N || adj.abs() < prev {
                    state.tx = state.xinbta - adj; // x_{n+1} = x_n - g*w
                    if 0.0 <= state.tx && state.tx <= 1.0 {
                        if prev <= state.acu || state.w.abs() <= state.acu {
                            // R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g, %s <= acu  ==> convergence\n",
                            //                  i_inn, -adj, (prev <= acu) ? "prev" : "|w|");
                            l_converged(state, log);
                            if state.exit {
                                return;
                            }
                        }
                        if state.tx != 0.0 && state.tx != 1.0 {
                            break;
                        }
                    }
                }
                state.g /= 3.0;
                i_inn += 1;
            }
            // R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g\n", i_inn, tx - xinbta);
            if (state.tx - state.xinbta).abs() <= 4e-16 * (state.tx + state.xinbta) {
                // "<=" : (.) == 0
                l_converged(state, log);
                if state.exit {
                    return;
                }
            }
            state.xinbta = state.tx;
            if state.tx == 0.0 {
                // "we have lost"
                break;
            }
            wprev = state.w;
            i_pb += 1;
        } // for( i_pb ..)
    } // end{else : normal scale Newton}

    /*-- NOT converged: Iteration count --*/
    state.warned = true;
    // ML_WARNING(ME_PRECISION, "qbeta");
}

fn l_converged(state: &mut QBetaState, log: bool) {
    state.log_ = log || state.use_log_x;
    // R_ifDEBUG_printf(" %s: Final delta(y) = %g%s\n",
    //                  warned ? "_NO_ convergence" : "converged",
    //                  y - (log_ ? la : a), (log_ ? " (log_)" : ""));
    if (state.log_ && state.y == f64::neg_infinity()) || (!state.log_ && state.y == 0.0) {
        // stuck at left, try if smallest positive number is "better"
        state.w = pbeta_raw(DBL_very_MIN, state.pp, state.qq, true, state.log_);
        if state.log_ || (state.w - state.a).abs() <= (state.y - state.a).abs() {
            state.tx = DBL_very_MIN;
            state.u_n = DBL_log_v_MIN; // = log(DBL_very_MIN)
        }
        state.add_N_step = false; // not trying to do better anymore
    } else if !state.warned
        && if state.log_ {
            (state.y - state.la).abs() > 3.0
        } else {
            (state.y - state.a).abs() > 1e-4
        }
        && !(state.log_ && state.y == f64::neg_infinity() &&
    // e.g. qbeta(-1e-10, .2, .03, log=TRUE) cannot get accurate ==> do NOT warn
    pbeta_raw(DBL_1__eps, // = 1 - eps
              state.pp, state.qq, true, true) > state.la + 2.0)
    {
        // MATHLIB_WARNING2( // low accuracy for more platform independent output:
        // "qbeta(a, *) =: x0 with |pbeta(x0,*%s) - alpha| = %.5g is not accurate",
        // (log_ ? ", log_" : ""), fabs(y - (log_ ? la : a)));
    }
}

fn l_return(state: &mut QBetaState, log: bool, qb: &mut [f64]) {
    // use  if (use_log_x) u_n else tx   {and u_n is on log scale}
    if state.give_log_q {
        // {currently not used from R's qbeta()} ==> use_log_x , too
        if !state.use_log_x { // (see if claim above is true)
             // MATHLIB_WARNING(
             //     "qbeta() L_return, u_n=%g;  give_log_q=TRUE but use_log_x=FALSE -- please report!",
             //     u_n);
        }
        let r = state.u_n.log1_exp();
        if state.swap_tail {
            qb[0] = r;
            qb[1] = state.u_n;
        } else {
            qb[0] = state.u_n;
            qb[1] = r;
        }
    } else {
        if state.use_log_x {
            if state.add_N_step {
                /* add one last Newton step on original x scale, e.g., for
                qbeta(2^-98, 0.125, 2^-96) */
                if state.u_n != 1.0 {
                    // u_n has been computed above
                    state.xinbta = state.u_n.exp();
                }
                state.y = pbeta_raw(
                    state.xinbta,
                    state.pp,
                    state.qq,
                    /*lower_tail = */ true,
                    log,
                );
                state.w = if log {
                    (state.y - state.la)
                        * (state.y
                            + state.logbeta
                            + state.r * state.xinbta.ln()
                            + state.t * (-state.xinbta).ln_1p())
                        .exp()
                } else {
                    (state.y - state.a)
                        * (state.logbeta
                            + state.r * state.xinbta.ln()
                            + state.t * (-state.xinbta).ln_1p())
                        .exp()
                };
                if state.w.is_finite() {
                    state.tx = state.xinbta - state.w;
                    // R_ifDEBUG_printf(" Final Newton correction(non-log scale):\n"
                    //                  "  xinbta=%.16g, y=%g, w=-Delta(x)=%g. \n=> new x=%.16g\n",
                    //                  xinbta, y, w, tx);
                } else {
                    // Newton step w  cannot be used
                    // R_ifDEBUG_printf(" Final Newton correction(non-log scale), canNOT use step 'w':\n"
                    //                  "  xinbta=%.16g, y=%g, w=-Delta(x)=%g\n=> keep x=xinbta\n", xinbta, y, w);
                    state.tx = state.xinbta;
                }
            } else {
                if state.swap_tail {
                    qb[0] = -state.u_n.exp_m1();
                    qb[1] = state.u_n.exp();
                } else {
                    qb[0] = state.u_n.exp();
                    qb[1] = -state.u_n.exp_m1();
                }
                return;
            }
        }
        if state.swap_tail {
            qb[0] = 1.0 - state.tx;
            qb[1] = state.tx;
        } else {
            qb[0] = state.tx;
            qb[1] = 1.0 - state.tx;
        }
    }
    return;
}