r2rs-stats 0.1.1

use crate::sexprec::{SEXP, SEXPREC, SEXPTYPE};
use ::libc;
extern "C" {
    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;
    #[no_mangle]
    fn memcpy(_: *mut libc::c_void, _: *const libc::c_void, _: libc::c_ulong) -> *mut libc::c_void;
    #[no_mangle]
    fn memset(_: *mut libc::c_void, _: libc::c_int, _: libc::c_ulong) -> *mut libc::c_void;
    #[no_mangle]
    fn R_chk_free(_: *mut libc::c_void);
    #[no_mangle]
    fn R_chk_calloc(_: size_t, _: size_t) -> *mut libc::c_void;
    /* Random Number Generators */
    #[no_mangle]
    fn norm_rand() -> libc::c_double;
    #[no_mangle]
    fn PutRNGstate();
    #[no_mangle]
    fn GetRNGstate();
    #[no_mangle]
    fn error(_: *const libc::c_char, _: ...) -> !;
    #[no_mangle]
    fn isReal(s: SEXP) -> Rboolean;
    #[no_mangle]
    fn INTEGER(x: SEXP) -> *mut libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    /* "dimnames" */
    #[no_mangle]
    static mut R_DimSymbol: SEXP;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn asReal(x: SEXP) -> libc::c_double;
    #[no_mangle]
    fn alloc3DArray(_: SEXPTYPE, _: libc::c_int, _: libc::c_int, _: libc::c_int) -> SEXP;
    #[no_mangle]
    fn getAttrib(_: SEXP, _: SEXP) -> SEXP;
    #[no_mangle]
    fn isMatrix(_: SEXP) -> Rboolean;
    #[no_mangle]
    fn protect(_: SEXP) -> SEXP;
    #[no_mangle]
    fn unprotect(_: libc::c_int);
    #[no_mangle]
    fn rchisq(_: libc::c_double) -> libc::c_double;
    /* DPOTRF - compute the Cholesky factorization of a real */
    /* symmetric positive definite matrix A */
    #[no_mangle]
    fn dpotrf_(
        uplo: *const libc::c_char,
        n: *const libc::c_int,
        a: *mut libc::c_double,
        lda: *const libc::c_int,
        info: *mut libc::c_int,
    );
    #[no_mangle]
    fn dsyrk_(
        uplo: *const libc::c_char,
        trans: *const libc::c_char,
        n: *const libc::c_int,
        k: *const libc::c_int,
        alpha: *const libc::c_double,
        a: *const libc::c_double,
        lda: *const libc::c_int,
        beta: *const libc::c_double,
        c: *mut libc::c_double,
        ldc: *const libc::c_int,
    );
    #[no_mangle]
    fn dtrmm_(
        side: *const libc::c_char,
        uplo: *const libc::c_char,
        transa: *const libc::c_char,
        diag: *const libc::c_char,
        m: *const libc::c_int,
        n: *const libc::c_int,
        alpha: *const libc::c_double,
        a: *const libc::c_double,
        lda: *const libc::c_int,
        b: *mut libc::c_double,
        ldb: *const libc::c_int,
    );
    #[no_mangle]
    fn dcgettext(
        __domainname: *const libc::c_char,
        __msgid: *const libc::c_char,
        __category: libc::c_int,
    ) -> *mut libc::c_char;
}
pub type size_t = libc::c_ulong;
pub type Rboolean = libc::c_uint;
pub const TRUE: Rboolean = 1;
pub const FALSE: Rboolean = 0;
/* Fundamental Data Types:  These are largely Lisp
 * influenced structures, with the exception of LGLSXP,
 * INTSXP, REALSXP, CPLXSXP and STRSXP which are the
 * element types for S-like data objects.
 *
 *   --> TypeTable[] in ../main/util.c for  typeof()
 */
/* UUID identifying the internals version -- packages using compiled
code should be re-installed when this changes */
/*  These exact numeric values are seldom used, but they are, e.g., in
 *  ../main/subassign.c, and they are serialized.
*/
/* NOT YET using enum:
 *  1) The SEXPREC struct below has 'SEXPTYPE type : 5'
 * (making FUNSXP and CLOSXP equivalent in there),
 * giving (-Wall only ?) warnings all over the place
 * 2) Many switch(type) { case ... } statements need a final `default:'
 * added in order to avoid warnings like [e.g. l.170 of ../main/util.c]
 *   "enumeration value `FUNSXP' not handled in switch"
 */
/* 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) 2012-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/
 */
// memset, memcpy
/* *
* Simulate the Cholesky factor of a standardized Wishart variate with
* dimension p and nu degrees of freedom.
*
* @param nu degrees of freedom
* @param p dimension of the Wishart distribution
* @param upper if 0 the result is lower triangular, otherwise upper
               triangular
* @param ans array of size p * p to hold the result
*
* @return ans
*/
unsafe extern "C" fn std_rWishart_factor(
    mut nu: libc::c_double,
    mut p: libc::c_int,
    mut upper: libc::c_int,
    mut ans: *mut libc::c_double,
) -> *mut libc::c_double {
    let mut pp1: libc::c_int = p + 1 as libc::c_int;
    if nu < p as libc::c_double || p <= 0 as libc::c_int {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"inconsistent degrees of freedom and dimension\x00" as *const u8
                as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    memset(
        ans as *mut libc::c_void,
        0 as libc::c_int,
        ((p * p) as libc::c_ulong)
            .wrapping_mul(::std::mem::size_of::<libc::c_double>() as libc::c_ulong),
    );
    let mut j: libc::c_int = 0 as libc::c_int;
    while j < p {
        /* jth column */
        *ans.offset((j * pp1) as isize) = sqrt(rchisq(nu - j as libc::c_double)); /* lower triangle index */
        let mut i: libc::c_int = 0 as libc::c_int;
        while i < j {
            let mut uind: libc::c_int = i + j * p;
            let mut lind: libc::c_int = j + i * p;
            *ans.offset(if upper != 0 { uind } else { lind } as isize) = norm_rand();
            *ans.offset(if upper != 0 { lind } else { uind } as isize) =
                0 as libc::c_int as libc::c_double;
            i += 1
        }
        j += 1
    }
    return ans;
}
/* *
 * Simulate a sample of random matrices from a Wishart distribution
 *
 * @param ns Number of samples to generate
 * @param nuP Degrees of freedom
 * @param scal Positive-definite scale matrix
 *
 * @return
 */
#[no_mangle]
pub unsafe extern "C" fn rWishart(mut ns: SEXP, mut nuP: SEXP, mut scal: SEXP) -> SEXP {
    let mut ans: SEXP = 0 as *mut SEXPREC;
    let mut dims: *mut libc::c_int = INTEGER(getAttrib(scal, R_DimSymbol));
    let mut info: libc::c_int = 0;
    let mut n: libc::c_int = asInteger(ns);
    let mut psqr: libc::c_int = 0;
    let mut scCp: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut ansp: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut tmp: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut nu: libc::c_double = asReal(nuP);
    let mut one: libc::c_double = 1 as libc::c_int as libc::c_double;
    let mut zero: libc::c_double = 0 as libc::c_int as libc::c_double;
    if isMatrix(scal) as u64 == 0
        || isReal(scal) as u64 == 0
        || *dims.offset(0 as libc::c_int as isize) != *dims.offset(1 as libc::c_int as isize)
    {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"\'scal\' must be a square, real matrix\x00" as *const u8 as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    if n <= 0 as libc::c_int {
        n = 1 as libc::c_int
    }
    // allocate early to avoid memory leaks in Callocs below.
    ans = alloc3DArray(
        14 as libc::c_int as SEXPTYPE,
        *dims.offset(0 as libc::c_int as isize),
        *dims.offset(0 as libc::c_int as isize),
        n,
    );
    protect(ans);
    psqr = *dims.offset(0 as libc::c_int as isize) * *dims.offset(0 as libc::c_int as isize);
    tmp = R_chk_calloc(
        psqr as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong,
    ) as *mut libc::c_double;
    scCp = R_chk_calloc(
        psqr as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong,
    ) as *mut libc::c_double;
    memcpy(
        scCp as *mut libc::c_void,
        REAL(scal) as *const libc::c_void,
        (psqr as size_t).wrapping_mul(::std::mem::size_of::<libc::c_double>() as libc::c_ulong),
    );
    memset(
        tmp as *mut libc::c_void,
        0 as libc::c_int,
        (psqr as libc::c_ulong)
            .wrapping_mul(::std::mem::size_of::<libc::c_double>() as libc::c_ulong),
    );
    dpotrf_(
        b"U\x00" as *const u8 as *const libc::c_char,
        &mut *dims.offset(0 as libc::c_int as isize),
        scCp,
        &mut *dims.offset(0 as libc::c_int as isize),
        &mut info,
    );
    if info != 0 {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"\'scal\' matrix is not positive-definite\x00" as *const u8 as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    ansp = REAL(ans);
    GetRNGstate();
    let mut j: libc::c_int = 0 as libc::c_int;
    while j < n {
        let mut ansj: *mut libc::c_double = ansp.offset((j * psqr) as isize);
        std_rWishart_factor(
            nu,
            *dims.offset(0 as libc::c_int as isize),
            1 as libc::c_int,
            tmp,
        );
        dtrmm_(
            b"R\x00" as *const u8 as *const libc::c_char,
            b"U\x00" as *const u8 as *const libc::c_char,
            b"N\x00" as *const u8 as *const libc::c_char,
            b"N\x00" as *const u8 as *const libc::c_char,
            dims,
            dims,
            &mut one,
            scCp,
            dims,
            tmp,
            dims,
        );
        dsyrk_(
            b"U\x00" as *const u8 as *const libc::c_char,
            b"T\x00" as *const u8 as *const libc::c_char,
            &mut *dims.offset(1 as libc::c_int as isize),
            &mut *dims.offset(1 as libc::c_int as isize),
            &mut one,
            tmp,
            &mut *dims.offset(1 as libc::c_int as isize),
            &mut zero,
            ansj,
            &mut *dims.offset(1 as libc::c_int as isize),
        );
        let mut i: libc::c_int = 1 as libc::c_int;
        while i < *dims.offset(0 as libc::c_int as isize) {
            let mut k: libc::c_int = 0 as libc::c_int;
            while k < i {
                *ansj.offset((i + k * *dims.offset(0 as libc::c_int as isize)) as isize) =
                    *ansj.offset((k + i * *dims.offset(0 as libc::c_int as isize)) as isize);
                k += 1
            }
            i += 1
        }
        j += 1
    }
    PutRNGstate();
    R_chk_free(scCp as *mut libc::c_void);
    scCp = 0 as *mut libc::c_double;
    R_chk_free(tmp as *mut libc::c_void);
    tmp = 0 as *mut libc::c_double;
    unprotect(1 as libc::c_int);
    return ans;
}