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 R_alloc(_: size_t, _: libc::c_int) -> *mut libc::c_char;
    #[no_mangle]
    fn R_chk_free(_: *mut libc::c_void);
    #[no_mangle]
    fn R_chk_calloc(_: size_t, _: size_t) -> *mut libc::c_void;
    /* Vector Access Functions */
    #[no_mangle]
    fn LENGTH(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn asReal(x: SEXP) -> libc::c_double;
    #[no_mangle]
    fn duplicate(_: SEXP) -> SEXP;
    #[no_mangle]
    fn ScalarReal(_: libc::c_double) -> SEXP;
    #[no_mangle]
    fn R_pow_di(_: libc::c_double, _: libc::c_int) -> libc::c_double;
}
pub type size_t = libc::c_ulong;
/* 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 */
/* Two-sample two-sided asymptotic distribution */
unsafe extern "C" fn pkstwo(
    mut n: libc::c_int,
    mut x: *mut libc::c_double,
    mut tol: libc::c_double,
) {
    /* x[1:n] is input and output
     *
     * Compute
     *   \sum_{k=-\infty}^\infty (-1)^k e^{-2 k^2 x^2}
     *   = 1 + 2 \sum_{k=1}^\infty (-1)^k e^{-2 k^2 x^2}
     *   = \frac{\sqrt{2\pi}}{x} \sum_{k=1}^\infty \exp(-(2k-1)^2\pi^2/(8x^2))
     *
     * See e.g. J. Durbin (1973), Distribution Theory for Tests Based on the
     * Sample Distribution Function.  SIAM.
     *
     * The 'standard' series expansion obviously cannot be used close to 0;
     * we use the alternative series for x < 1, and a rather crude estimate
     * of the series remainder term in this case, in particular using that
     * ue^(-lu^2) \le e^(-lu^2 + u) \le e^(-(l-1)u^2 - u^2+u) \le e^(-(l-1))
     * provided that u and l are >= 1.
     *
     * (But note that for reasonable tolerances, one could simply take 0 as
     * the value for x < 0.2, and use the standard expansion otherwise.)
     *
     */
    let mut new: libc::c_double = 0.;
    let mut old: libc::c_double = 0.;
    let mut s: libc::c_double = 0.;
    let mut w: libc::c_double = 0.;
    let mut z: libc::c_double = 0.;
    let mut i: libc::c_int = 0;
    let mut k: libc::c_int = 0;
    let mut k_max: libc::c_int = 0;
    k_max = sqrt(2 as libc::c_int as libc::c_double - tol.ln()) as libc::c_int;
    i = 0 as libc::c_int;
    while i < n {
        if *x.offset(i as isize) < 1 as libc::c_int as libc::c_double {
            z = -(1.57079632679489661923f64 * 0.78539816339744830962f64)
                / (*x.offset(i as isize) * *x.offset(i as isize));
            w = (*x.offset(i as isize)).ln();
            s = 0 as libc::c_int as libc::c_double;
            k = 1 as libc::c_int;
            while k < k_max {
                s += ((k * k) as libc::c_double * z - w).exp();
                k += 2 as libc::c_int
            }
            *x.offset(i as isize) = s / 0.398942280401432677939946059934f64
        } else {
            z = -(2 as libc::c_int) as libc::c_double
                * *x.offset(i as isize)
                * *x.offset(i as isize);
            s = -(1 as libc::c_int) as libc::c_double;
            k = 1 as libc::c_int;
            old = 0 as libc::c_int as libc::c_double;
            new = 1 as libc::c_int as libc::c_double;
            while (old - new).abs() > tol {
                old = new;
                new += 2 as libc::c_int as libc::c_double
                    * s
                    * (z * k as libc::c_double * k as libc::c_double).exp();
                s *= -(1 as libc::c_int) as libc::c_double;
                k += 1
            }
            *x.offset(i as isize) = new
        }
        i += 1
    }
}
/* Two-sided two-sample */
unsafe extern "C" fn psmirnov2x(
    mut x: *mut libc::c_double,
    mut m: libc::c_int,
    mut n: libc::c_int,
) -> libc::c_double {
    let mut md: libc::c_double = 0.;
    let mut nd: libc::c_double = 0.;
    let mut q: libc::c_double = 0.;
    let mut u: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut w: libc::c_double = 0.;
    let mut i: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    if m > n {
        i = n;
        n = m;
        m = i
    }
    md = m as libc::c_double;
    nd = n as libc::c_double;
    /*
       q has 0.5/mn added to ensure that rounding error doesn't
       turn an equality into an inequality, eg abs(1/2-4/5)>3/10

    */
    q = (0.5f64 + (*x * md * nd - 1e-7f64).floor()) / (md * nd);
    u = R_alloc(
        (n + 1 as libc::c_int) as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong as libc::c_int,
    ) as *mut libc::c_double;
    j = 0 as libc::c_int;
    while j <= n {
        *u.offset(j as isize) = if j as libc::c_double / nd > q {
            0 as libc::c_int
        } else {
            1 as libc::c_int
        } as libc::c_double;
        j += 1
    }
    i = 1 as libc::c_int;
    while i <= m {
        w = i as libc::c_double / (i + n) as libc::c_double;
        if i as libc::c_double / md > q {
            *u.offset(0 as libc::c_int as isize) = 0 as libc::c_int as libc::c_double
        } else {
            *u.offset(0 as libc::c_int as isize) = w * *u.offset(0 as libc::c_int as isize)
        }
        j = 1 as libc::c_int;
        while j <= n {
            if (i as libc::c_double / md - j as libc::c_double / nd).abs() > q {
                *u.offset(j as isize) = 0 as libc::c_int as libc::c_double
            } else {
                *u.offset(j as isize) =
                    w * *u.offset(j as isize) + *u.offset((j - 1 as libc::c_int) as isize)
            }
            j += 1
        }
        i += 1
    }
    return *u.offset(n as isize);
}
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1999-2016   The R Core Team.
 *
 *  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/
 */
/* ks.c
   Compute the asymptotic distribution of the one- and two-sample
   two-sided Kolmogorov-Smirnov statistics, and the exact distributions
   in the two-sided one-sample and two-sample cases.
*/
unsafe extern "C" fn K(mut n: libc::c_int, mut d: libc::c_double) -> libc::c_double {
    /* Compute Kolmogorov's distribution.
       Code published in
     George Marsaglia and Wai Wan Tsang and Jingbo Wang (2003),
     "Evaluating Kolmogorov's distribution".
     Journal of Statistical Software, Volume 8, 2003, Issue 18.
     URL: http://www.jstatsoft.org/v08/i18/.
    */
    let mut k: libc::c_int = 0;
    let mut m: libc::c_int = 0;
    let mut i: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut g: libc::c_int = 0;
    let mut eH: libc::c_int = 0;
    let mut eQ: libc::c_int = 0;
    let mut h: libc::c_double = 0.;
    let mut s: libc::c_double = 0.;
    let mut H: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut Q: *mut libc::c_double = 0 as *mut libc::c_double;
    /*
       The faster right-tail approximation is omitted here.
       s = d*d*n;
       if(s > 7.24 || (s > 3.76 && n > 99))
           return 1-2*(-(2.000071+.331/sqrt(n)+1.409/n)*s).exp();
    */
    k = (n as libc::c_double * d) as libc::c_int + 1 as libc::c_int;
    m = 2 as libc::c_int * k - 1 as libc::c_int;
    h = k as libc::c_double - n as libc::c_double * d;
    H = R_chk_calloc(
        (m * m) as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong,
    ) as *mut libc::c_double;
    Q = R_chk_calloc(
        (m * m) as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong,
    ) as *mut libc::c_double;
    i = 0 as libc::c_int;
    while i < m {
        j = 0 as libc::c_int;
        while j < m {
            if (i - j + 1 as libc::c_int) < 0 as libc::c_int {
                *H.offset((i * m + j) as isize) = 0 as libc::c_int as libc::c_double
            } else {
                *H.offset((i * m + j) as isize) = 1 as libc::c_int as libc::c_double
            }
            j += 1
        }
        i += 1
    }
    i = 0 as libc::c_int;
    while i < m {
        *H.offset((i * m) as isize) -= R_pow_di(h, i + 1 as libc::c_int);
        *H.offset(((m - 1 as libc::c_int) * m + i) as isize) -= R_pow_di(h, m - i);
        i += 1
    }
    *H.offset(((m - 1 as libc::c_int) * m) as isize) += if 2 as libc::c_int as libc::c_double * h
        - 1 as libc::c_int as libc::c_double
        > 0 as libc::c_int as libc::c_double
    {
        R_pow_di(
            2 as libc::c_int as libc::c_double * h - 1 as libc::c_int as libc::c_double,
            m,
        )
    } else {
        0 as libc::c_int as libc::c_double
    };
    i = 0 as libc::c_int;
    while i < m {
        j = 0 as libc::c_int;
        while j < m {
            if i - j + 1 as libc::c_int > 0 as libc::c_int {
                g = 1 as libc::c_int;
                while g <= i - j + 1 as libc::c_int {
                    *H.offset((i * m + j) as isize) /= g as libc::c_double;
                    g += 1
                }
            }
            j += 1
        }
        i += 1
    }
    eH = 0 as libc::c_int;
    m_power(H, eH, Q, &mut eQ, m, n);
    s = *Q.offset(((k - 1 as libc::c_int) * m + k - 1 as libc::c_int) as isize);
    i = 1 as libc::c_int;
    while i <= n {
        s = s * i as libc::c_double / n as libc::c_double;
        if s < 1e-140f64 {
            s *= 1e140f64;
            eQ -= 140 as libc::c_int
        }
        i += 1
    }
    s *= R_pow_di(10.0f64, eQ);
    R_chk_free(H as *mut libc::c_void);
    H = 0 as *mut libc::c_double;
    R_chk_free(Q as *mut libc::c_void);
    Q = 0 as *mut libc::c_double;
    return s;
}
unsafe extern "C" fn m_multiply(
    mut A: *mut libc::c_double,
    mut B: *mut libc::c_double,
    mut C: *mut libc::c_double,
    mut m: libc::c_int,
) {
    /* Auxiliary routine used by K().
       Matrix multiplication.
    */
    let mut i: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut k: libc::c_int = 0;
    let mut s: libc::c_double = 0.;
    i = 0 as libc::c_int;
    while i < m {
        j = 0 as libc::c_int;
        while j < m {
            s = 0.0f64;
            k = 0 as libc::c_int;
            while k < m {
                s += *A.offset((i * m + k) as isize) * *B.offset((k * m + j) as isize);
                k += 1
            }
            *C.offset((i * m + j) as isize) = s;
            j += 1
        }
        i += 1
    }
}
unsafe extern "C" fn m_power(
    mut A: *mut libc::c_double,
    mut eA: libc::c_int,
    mut V: *mut libc::c_double,
    mut eV: *mut libc::c_int,
    mut m: libc::c_int,
    mut n: libc::c_int,
) {
    /* Auxiliary routine used by K().
       Matrix power.
    */
    let mut B: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut eB: libc::c_int = 0;
    let mut i: libc::c_int = 0;
    if n == 1 as libc::c_int {
        i = 0 as libc::c_int;
        while i < m * m {
            *V.offset(i as isize) = *A.offset(i as isize);
            i += 1
        }
        *eV = eA;
        return;
    }
    m_power(A, eA, V, eV, m, n / 2 as libc::c_int);
    B = R_chk_calloc(
        (m * m) as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong,
    ) as *mut libc::c_double;
    m_multiply(V, V, B, m);
    eB = 2 as libc::c_int * *eV;
    if n % 2 as libc::c_int == 0 as libc::c_int {
        i = 0 as libc::c_int;
        while i < m * m {
            *V.offset(i as isize) = *B.offset(i as isize);
            i += 1
        }
        *eV = eB
    } else {
        m_multiply(A, B, V, m);
        *eV = eA + eB
    }
    if *V.offset((m / 2 as libc::c_int * m + m / 2 as libc::c_int) as isize) > 1e140f64 {
        i = 0 as libc::c_int;
        while i < m * m {
            *V.offset(i as isize) = *V.offset(i as isize) * 1e-140f64;
            i += 1
        }
        *eV += 140 as libc::c_int
    }
    R_chk_free(B as *mut libc::c_void);
    B = 0 as *mut libc::c_double;
}
/* Two-sided two-sample */
#[no_mangle]
pub unsafe extern "C" fn pSmirnov2x(mut statistic: SEXP, mut snx: SEXP, mut sny: SEXP) -> SEXP {
    let mut nx: libc::c_int = asInteger(snx);
    let mut ny: libc::c_int = asInteger(sny);
    let mut st: libc::c_double = asReal(statistic);
    return ScalarReal(psmirnov2x(&mut st, nx, ny));
}
/* Two-sample two-sided asymptotic distribution */
#[no_mangle]
pub unsafe extern "C" fn pKS2(mut statistic: SEXP, mut stol: SEXP) -> SEXP {
    let mut n: libc::c_int = LENGTH(statistic);
    let mut tol: libc::c_double = asReal(stol);
    let mut ans: SEXP = duplicate(statistic);
    pkstwo(n, REAL(ans), tol);
    return ans;
}
/* The two-sided one-sample 'exact' distribution */
#[no_mangle]
pub unsafe extern "C" fn pKolmogorov2x(mut statistic: SEXP, mut sn: SEXP) -> SEXP {
    let mut n: libc::c_int = asInteger(sn);
    let mut st: libc::c_double = asReal(statistic);
    let mut p: libc::c_double = 0.;
    p = K(n, st);
    return ScalarReal(p);
}