r2rs-stats 0.1.1

use crate::sexprec::SEXP;
use ::libc;
extern "C" {
    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;
    #[no_mangle]
    fn pnorm(
        _: libc::c_double,
        _: libc::c_double,
        _: libc::c_double,
        _: libc::c_int,
        _: libc::c_int,
    ) -> libc::c_double;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn asReal(x: SEXP) -> libc::c_double;
    #[no_mangle]
    fn ScalarReal(_: libc::c_double) -> SEXP;
}
pub type C2RustUnnamed = libc::c_uint;
pub const TRUE: C2RustUnnamed = 1;
pub const FALSE: C2RustUnnamed = 0;
/* nil = NULL */
/* symbols */
/* lists of dotted pairs */
/* closures */
/* environments */
/* promises: [un]evaluated closure arguments */
/* language constructs (special lists) */
/* special forms */
/* builtin non-special forms */
/* "scalar" string type (internal only)*/
/* logical vectors */
/* 11 and 12 were factors and ordered factors in the 1990s */
/* integer vectors */
/* real variables */
/* complex variables */
/* string vectors */
/* dot-dot-dot object */
/* make "any" args work.
Used in specifying types for symbol
registration to mean anything is okay  */
/* generic vectors */
/* expressions vectors */
/* byte code */
/* external pointer */
/* weak reference */
/* raw bytes */
/* S4, non-vector */
/* used for detecting PROTECT issues in memory.c */
/* fresh node created in new page */
/* node released by GC */
/* Closure or Builtin or Special */
/* NOT YET */
/* These are also used with the write barrier on, in attrib.c and util.c */
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2000-2016     The R Core Team
 *  Copyright (C) 2003  The R Foundation
 *  based on AS 89 (C) 1975 Royal Statistical Society
 *
 *  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 2 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, a copy is available at
 *  https://www.R-project.org/Licenses/
 *
 */
/* Was
    double precision function prho(n, is, ifault)

 Changed to subroutine by KH to allow for .Fortran interfacing.
 Also, change `prho' argument in the code to `pv'.
 And, fix a bug.

 From R ver. 1.1.x [March, 2000] by MM:
 - Translate Fortran to C
 - use pnorm() instead of less precise alnorm().
 - new argument lower_tail --> potentially increased precision in extreme cases.
*/
#[no_mangle]
pub unsafe extern "C" fn prho(
    mut n: libc::c_int,
    mut is: libc::c_double,
    mut pv: *mut libc::c_double,
    mut ifault: *mut libc::c_int,
    mut lower_tail: libc::c_int,
) {
    /* Algorithm AS 89   Appl. Statist. (1975) Vol.24, No. 3, P377.

        To evaluate the probability  Pr[ S >= is ]
        {or Pr [ S < is]  if(lower_tail) }, where

        S   = (n^3 - n) * (1-R)/6,
        is  = (n^3 - n) * (1-r)/6,
        R,r = Spearman's rho (r.v. and observed),  and n >= 2
    */
    /* Edgeworth coefficients : */
    let c1: libc::c_double = 0.2274f64;
    let c2: libc::c_double = 0.2531f64;
    let c3: libc::c_double = 0.1745f64;
    let c4: libc::c_double = 0.0758f64;
    let c5: libc::c_double = 0.1033f64;
    let c6: libc::c_double = 0.3932f64;
    let c7: libc::c_double = 0.0879f64;
    let c8: libc::c_double = 0.0151f64;
    let c9: libc::c_double = 0.0072f64;
    let c10: libc::c_double = 0.0831f64;
    let c11: libc::c_double = 0.0131f64;
    let c12: libc::c_double = 4.6e-4f64;
    /* Local variables */
    let mut b: libc::c_double = 0.; /*, js */
    let mut u: libc::c_double = 0.;
    let mut x: libc::c_double = 0.;
    let mut y: libc::c_double = 0.;
    let mut n3: libc::c_double = 0.;
    /* originally: n_small = 6 (speed!);
     * even n_small = 10 (and n = 10) needs quite a bit longer than the approx!
     * larger than 12 ==> integer overflow in nfac and (probably) ifr
     */
    let mut l: [libc::c_int; 9] = [0; 9];
    let mut nfac: libc::c_int = 0;
    let mut i: libc::c_int = 0;
    let mut m: libc::c_int = 0;
    let mut mt: libc::c_int = 0;
    let mut ifr: libc::c_int = 0;
    let mut ise: libc::c_int = 0;
    let mut n1: libc::c_int = 0;
    /* Test admissibility of arguments and initialize */
    *pv = if lower_tail != 0 { 0.0f64 } else { 1.0f64 }; /* with p = 1 */
    if n <= 1 as libc::c_int {
        *ifault = 1 as libc::c_int; /* = (n^3 - n)/3 */
        return;
    }
    *ifault = 0 as libc::c_int;
    if is <= 0.0f64 {
        return;
    }
    n3 = n as libc::c_double;
    n3 *= (n3 * n3 - 1.0f64) / 3.0f64;
    if is > n3 {
        /* larger than maximal value */
        *pv = 1 as libc::c_int as libc::c_double - *pv;
        return;
    }
    /* NOT rounding to even anymore:  with ties, S, may even be non-integer!
     * js = is;
     * if(fmod(js, 2.) != 0.) ++js;
     */
    if n <= 9 as libc::c_int {
        nfac = 1.0f64 as libc::c_int; /* exact for n <= n_small */
        i = 1 as libc::c_int;
        while i <= n {
            nfac *= i;
            l[(i - 1 as libc::c_int) as usize] = i;
            i += 1
        }
        /* 2 <= n <= n_small :
         * Exact evaluation of probability */
        /* KH mod next line: was `!=' in the code but `.eq.' in the paper */
        if is == n3 {
            ifr = 1 as libc::c_int
        } else {
            ifr = 0 as libc::c_int;
            m = 0 as libc::c_int;
            while m < nfac {
                ise = 0 as libc::c_int;
                i = 0 as libc::c_int;
                while i < n {
                    n1 = i + 1 as libc::c_int - l[i as usize];
                    ise += n1 * n1;
                    i += 1
                }
                if is <= ise as libc::c_double {
                    ifr += 1
                }
                n1 = n;
                loop {
                    mt = l[0 as libc::c_int as usize];
                    i = 1 as libc::c_int;
                    while i < n1 {
                        l[(i - 1 as libc::c_int) as usize] = l[i as usize];
                        i += 1
                    }
                    n1 -= 1;
                    l[n1 as usize] = mt;
                    if !(mt == n1 + 1 as libc::c_int && n1 > 1 as libc::c_int) {
                        break;
                    }
                }
                m += 1
            }
        }
        *pv = (if lower_tail != 0 { (nfac) - ifr } else { ifr }) as libc::c_double
            / nfac as libc::c_double
    } else {
        /* n >= 7 : Evaluation by Edgeworth series expansion */
        y = n as libc::c_double;
        b = 1 as libc::c_int as libc::c_double / y;
        x = (6.0f64 * (is - 1 as libc::c_int as libc::c_double) * b
            / (y * y - 1 as libc::c_int as libc::c_double)
            - 1 as libc::c_int as libc::c_double)
            * sqrt(y - 1 as libc::c_int as libc::c_double);
        /* = rho * sqrt(n-1)  ==  rho / sqrt(var(rho))  ~  (0,1) */
        y = x * x;
        u = x
            * b
            * (c1
                + b * (c2 + c3 * b)
                + y * (-c4 + b * (c5 + c6 * b)
                    - y * b * (c7 + c8 * b - y * (c9 - c10 * b + y * b * (c11 - c12 * y)))));
        y = u / (y / 2.0f64).exp();
        *pv = (if lower_tail != 0 { -y } else { y })
            + pnorm(x, 0.0f64, 1.0f64, lower_tail, FALSE as libc::c_int);
        /* above was call to alnorm() [algorithm AS 66] */
        if *pv < 0 as libc::c_int as libc::c_double {
            *pv = 0.0f64
        }
        if *pv > 1 as libc::c_int as libc::c_double {
            *pv = 1.0f64
        }
    };
}
/* prho */
#[no_mangle]
pub unsafe extern "C" fn pRho(mut q: SEXP, mut sn: SEXP, mut lower: SEXP) -> SEXP {
    let mut s: libc::c_double = asReal(q);
    let mut p: libc::c_double = 0.;
    let mut n: libc::c_int = asInteger(sn);
    let mut ltail: libc::c_int = asInteger(lower);
    let mut ifault: libc::c_int = 0 as libc::c_int;
    prho(n, s, &mut p, &mut ifault, ltail);
    return ScalarReal(p);
}