r2rs-stats 0.1.1

use crate::{
    optim::getListElement,
    sexprec::{SEXP, SEXPREC, SEXPTYPE},
};
use ::libc;
extern "C" {
    #[no_mangle]
    fn __errno_location() -> *mut libc::c_int;
    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;

    /* NA_REAL: IEEE */
    #[no_mangle]
    static mut R_NaInt: libc::c_int;
    #[no_mangle]
    fn error(_: *const libc::c_char, _: ...) -> !;
    #[no_mangle]
    fn R_alloc(_: size_t, _: libc::c_int) -> *mut libc::c_char;
    #[no_mangle]
    fn TYPEOF(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn XLENGTH(x: SEXP) -> R_xlen_t;
    #[no_mangle]
    fn INTEGER(x: SEXP) -> *mut libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    #[no_mangle]
    fn SET_STRING_ELT(x: SEXP, i: R_xlen_t, v: SEXP);
    #[no_mangle]
    fn SET_VECTOR_ELT(x: SEXP, i: R_xlen_t, v: SEXP) -> SEXP;
    /* "name" */
    #[no_mangle]
    static mut R_NamesSymbol: SEXP;
    #[no_mangle]
    fn coerceVector(_: SEXP, _: SEXPTYPE) -> SEXP;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn mkChar(_: *const libc::c_char) -> SEXP;
    #[no_mangle]
    fn setAttrib(_: SEXP, _: SEXP, _: SEXP) -> SEXP;
    /* Defining NO_RINLINEDFUNS disables use to simulate platforms where
    this is not available */
    /* need remapped names here for use with R_NO_REMAP */
    /*
       These are the inlinable functions that are provided in Rinlinedfuns.h
       It is *essential* that these do not appear in any other header file,
       with or without the Rf_ prefix.
    */
    #[no_mangle]
    fn allocVector(_: SEXPTYPE, _: R_xlen_t) -> SEXP;
    #[no_mangle]
    fn isVectorAtomic(_: SEXP) -> Rboolean;
    #[no_mangle]
    fn ScalarInteger(_: libc::c_int) -> SEXP;
    #[no_mangle]
    fn ScalarReal(_: libc::c_double) -> SEXP;
    #[no_mangle]
    fn protect(_: SEXP) -> SEXP;
    #[no_mangle]
    fn unprotect(_: libc::c_int);
    /*
     *  R : A Computer Language for Statistical Data Analysis
     *  Copyright (C) 2012   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/
     */
    /* auxiliary */

    /* ../main/bind.c */
    /* ../main/errors.c : */
    /* ../main/devices.c, used in memory.c, gnuwin32/extra.c */
    /* ../../main/printutils.c : */
    /* main/sort.c */
    /* main/subset.c */
    /* main/subassign.c */
    /* main/util.c */
    #[no_mangle]
    fn UNIMPLEMENTED_TYPE(s: *const libc::c_char, x: SEXP) -> !;
}
pub type size_t = libc::c_ulong;
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2000, 2001 The R Core Team.
 *
 *  This header file is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation; either version 2.1 of the License, or
 *  (at your option) any later version.
 *
 *  This file is part of R. R is distributed under the terms of the
 *  GNU General Public License, either Version 2, June 1991 or Version 3,
 *  June 2007. See doc/COPYRIGHTS for details of the copyright status of R.
 *
 *  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 Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 */
/* Included by R.h: API */
pub type Rboolean = libc::c_uint;
pub const TRUE: Rboolean = 1;
pub const FALSE: Rboolean = 0;
pub type ptrdiff_t = libc::c_long;
/*, MAYBE */
/* both config.h and Rconfig.h set SIZEOF_SIZE_T, but Rconfig.h is
skipped if config.h has already been included. */
pub type R_xlen_t = ptrdiff_t;
/* 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"
 */
/* 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) 1998--2017  The R Core Team
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *
 *  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/
 */
/* Spline Interpolation
 * --------------------
 * C code to perform spline regression and interpolation.
 * There is code here for:
 *
 * 1. Natural splines.
 *
 * 2. Periodic splines
 *
 * 3. Splines with end-conditions determined by regression
 *    cubics in the start and end intervals (Forsythe et al).
 *
 *
 * Computational Techniques
 * ------------------------
 * A special LU decomposition for symmetric tridiagonal matrices
 * is used for all computations, except for periodic splines where
 * Choleski is more efficient.
 */
/*
 * Natural Splines
 * ---------------
 * Here the end-conditions are determined by setting the second
 * derivative of the spline at the end-points to equal to zero.
 *
 * There are n-2 unknowns (y[i]'' at x[2], ..., x[n-1]) and n-2
 * equations to determine them.  Either Choleski or Gaussian
 * elimination could be used.
 */
unsafe extern "C" fn natural_spline(
    mut n: R_xlen_t,
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut c: *mut libc::c_double,
    mut d: *mut libc::c_double,
) {
    if n < 2 as libc::c_int as libc::c_long {
        *__errno_location() = 33 as libc::c_int;
        return;
    }
    x = x.offset(-1);
    y = y.offset(-1);
    b = b.offset(-1);
    c = c.offset(-1);
    d = d.offset(-1);
    if n < 3 as libc::c_int as libc::c_long {
        let mut t: libc::c_double =
            *y.offset(2 as libc::c_int as isize) - *y.offset(1 as libc::c_int as isize);
        *b.offset(1 as libc::c_int as isize) =
            t / (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
        *b.offset(2 as libc::c_int as isize) = *b.offset(1 as libc::c_int as isize);
        let ref mut fresh0 = *d.offset(2 as libc::c_int as isize);
        *fresh0 = 0.0f64;
        let ref mut fresh1 = *d.offset(1 as libc::c_int as isize);
        *fresh1 = *fresh0;
        let ref mut fresh2 = *c.offset(2 as libc::c_int as isize);
        *fresh2 = *fresh1;
        *c.offset(1 as libc::c_int as isize) = *fresh2;
        return;
    }
    let nm1: R_xlen_t = n - 1 as libc::c_int as libc::c_long;
    let mut i: R_xlen_t = 0;
    /* Set up the tridiagonal system */
    /* b = diagonal, d = offdiagonal, c = right hand side */
    *d.offset(1 as libc::c_int as isize) =
        *x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize);
    *c.offset(2 as libc::c_int as isize) = (*y.offset(2 as libc::c_int as isize)
        - *y.offset(1 as libc::c_int as isize))
        / *d.offset(1 as libc::c_int as isize);
    i = 2 as libc::c_int as R_xlen_t;
    while i < n {
        *d.offset(i as isize) =
            *x.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *x.offset(i as isize);
        *b.offset(i as isize) = 2.0f64
            * (*d.offset((i - 1 as libc::c_int as libc::c_long) as isize) + *d.offset(i as isize));
        *c.offset((i + 1 as libc::c_int as libc::c_long) as isize) =
            (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *y.offset(i as isize))
                / *d.offset(i as isize);
        *c.offset(i as isize) =
            *c.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *c.offset(i as isize);
        i += 1
    }
    /* Gaussian elimination */
    i = 3 as libc::c_int as R_xlen_t;
    while i < n {
        let mut t_0: libc::c_double = *d.offset((i - 1 as libc::c_int as libc::c_long) as isize)
            / *b.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        *b.offset(i as isize) = *b.offset(i as isize)
            - t_0 * *d.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        *c.offset(i as isize) = *c.offset(i as isize)
            - t_0 * *c.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        i += 1
    }
    /* Backward substitution */
    *c.offset(nm1 as isize) = *c.offset(nm1 as isize) / *b.offset(nm1 as isize);
    i = n - 2 as libc::c_int as libc::c_long;
    while i > 1 as libc::c_int as libc::c_long {
        *c.offset(i as isize) = (*c.offset(i as isize)
            - *d.offset(i as isize) * *c.offset((i + 1 as libc::c_int as libc::c_long) as isize))
            / *b.offset(i as isize);
        i -= 1
    }
    /* End conditions */
    let ref mut fresh3 = *c.offset(n as isize);
    *fresh3 = 0.0f64;
    *c.offset(1 as libc::c_int as isize) = *fresh3;
    /* Get cubic coefficients */
    *b.offset(1 as libc::c_int as isize) = (*y.offset(2 as libc::c_int as isize)
        - *y.offset(1 as libc::c_int as isize))
        / *d.offset(1 as libc::c_int as isize)
        - *d.offset(i as isize) * *c.offset(2 as libc::c_int as isize);
    *c.offset(1 as libc::c_int as isize) = 0.0f64;
    *d.offset(1 as libc::c_int as isize) =
        *c.offset(2 as libc::c_int as isize) / *d.offset(1 as libc::c_int as isize);
    *b.offset(n as isize) = (*y.offset(n as isize) - *y.offset(nm1 as isize))
        / *d.offset(nm1 as isize)
        + *d.offset(nm1 as isize) * *c.offset(nm1 as isize);
    i = 2 as libc::c_int as R_xlen_t;
    while i < n {
        *b.offset(i as isize) = (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *y.offset(i as isize))
            / *d.offset(i as isize)
            - *d.offset(i as isize)
                * (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
                    + 2.0f64 * *c.offset(i as isize));
        *d.offset(i as isize) = (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *c.offset(i as isize))
            / *d.offset(i as isize);
        *c.offset(i as isize) = 3.0f64 * *c.offset(i as isize);
        i += 1
    }
    *c.offset(n as isize) = 0.0f64;
    *d.offset(n as isize) = 0.0f64;
}
/*
 * Splines a la Forsythe Malcolm and Moler
 * ---------------------------------------
 * In this case the end-conditions are determined by regression
 * cubic polynomials to the first and last 4 points and matching
 * the third derivitives of the spline at the end-points to the
 * third derivatives of these cubics at the end-points.
 */
unsafe extern "C" fn fmm_spline(
    mut n: R_xlen_t,
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut c: *mut libc::c_double,
    mut d: *mut libc::c_double,
) {
    /* Adjustment for 1-based arrays */
    x = x.offset(-1);
    y = y.offset(-1);
    b = b.offset(-1);
    c = c.offset(-1);
    d = d.offset(-1);
    if n < 2 as libc::c_int as libc::c_long {
        *__errno_location() = 33 as libc::c_int;
        return;
    }
    if n < 3 as libc::c_int as libc::c_long {
        let mut t: libc::c_double =
            *y.offset(2 as libc::c_int as isize) - *y.offset(1 as libc::c_int as isize);
        *b.offset(1 as libc::c_int as isize) =
            t / (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
        *b.offset(2 as libc::c_int as isize) = *b.offset(1 as libc::c_int as isize);
        let ref mut fresh4 = *d.offset(2 as libc::c_int as isize);
        *fresh4 = 0.0f64;
        let ref mut fresh5 = *d.offset(1 as libc::c_int as isize);
        *fresh5 = *fresh4;
        let ref mut fresh6 = *c.offset(2 as libc::c_int as isize);
        *fresh6 = *fresh5;
        *c.offset(1 as libc::c_int as isize) = *fresh6;
        return;
    }
    let nm1: R_xlen_t = n - 1 as libc::c_int as libc::c_long;
    let mut i: R_xlen_t = 0;
    /* Set up tridiagonal system */
    /* b = diagonal, d = offdiagonal, c = right hand side */
    *d.offset(1 as libc::c_int as isize) =
        *x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize); /* = +/- Inf for x[1]=x[2] -- problem? */
    *c.offset(2 as libc::c_int as isize) = (*y.offset(2 as libc::c_int as isize)
        - *y.offset(1 as libc::c_int as isize))
        / *d.offset(1 as libc::c_int as isize);
    i = 2 as libc::c_int as R_xlen_t;
    while i < n {
        *d.offset(i as isize) =
            *x.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *x.offset(i as isize);
        *b.offset(i as isize) = 2.0f64
            * (*d.offset((i - 1 as libc::c_int as libc::c_long) as isize) + *d.offset(i as isize));
        *c.offset((i + 1 as libc::c_int as libc::c_long) as isize) =
            (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *y.offset(i as isize))
                / *d.offset(i as isize);
        *c.offset(i as isize) =
            *c.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *c.offset(i as isize);
        i += 1
    }
    /* End conditions. */
    /* Third derivatives at x[0] and x[n-1] obtained */
    /* from divided differences */
    *b.offset(1 as libc::c_int as isize) = -*d.offset(1 as libc::c_int as isize);
    *b.offset(n as isize) = -*d.offset(nm1 as isize);
    let ref mut fresh7 = *c.offset(n as isize);
    *fresh7 = 0.0f64;
    *c.offset(1 as libc::c_int as isize) = *fresh7;
    if n > 3 as libc::c_int as libc::c_long {
        *c.offset(1 as libc::c_int as isize) = *c.offset(3 as libc::c_int as isize)
            / (*x.offset(4 as libc::c_int as isize) - *x.offset(2 as libc::c_int as isize))
            - *c.offset(2 as libc::c_int as isize)
                / (*x.offset(3 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
        *c.offset(n as isize) = *c.offset(nm1 as isize)
            / (*x.offset(n as isize) - *x.offset((n - 2 as libc::c_int as libc::c_long) as isize))
            - *c.offset((n - 2 as libc::c_int as libc::c_long) as isize)
                / (*x.offset(nm1 as isize)
                    - *x.offset((n - 3 as libc::c_int as libc::c_long) as isize));
        *c.offset(1 as libc::c_int as isize) = *c.offset(1 as libc::c_int as isize)
            * *d.offset(1 as libc::c_int as isize)
            * *d.offset(1 as libc::c_int as isize)
            / (*x.offset(4 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
        *c.offset(n as isize) = -*c.offset(n as isize)
            * *d.offset(nm1 as isize)
            * *d.offset(nm1 as isize)
            / (*x.offset(n as isize) - *x.offset((n - 3 as libc::c_int as libc::c_long) as isize))
    }
    /* Gaussian elimination */
    i = 2 as libc::c_int as R_xlen_t;
    while i <= n {
        let mut t_0: libc::c_double = *d.offset((i - 1 as libc::c_int as libc::c_long) as isize)
            / *b.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        *b.offset(i as isize) = *b.offset(i as isize)
            - t_0 * *d.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        *c.offset(i as isize) = *c.offset(i as isize)
            - t_0 * *c.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        i += 1
    }
    /* Backward substitution */
    *c.offset(n as isize) = *c.offset(n as isize) / *b.offset(n as isize);
    i = nm1;
    while i >= 1 as libc::c_int as libc::c_long {
        *c.offset(i as isize) = (*c.offset(i as isize)
            - *d.offset(i as isize) * *c.offset((i + 1 as libc::c_int as libc::c_long) as isize))
            / *b.offset(i as isize);
        i -= 1
    }
    /* c[i] is now the sigma[i-1] of the text */
    /* Compute polynomial coefficients */
    *b.offset(n as isize) = (*y.offset(n as isize)
        - *y.offset((n - 1 as libc::c_int as libc::c_long) as isize))
        / *d.offset((n - 1 as libc::c_int as libc::c_long) as isize)
        + *d.offset((n - 1 as libc::c_int as libc::c_long) as isize)
            * (*c.offset((n - 1 as libc::c_int as libc::c_long) as isize)
                + 2.0f64 * *c.offset(n as isize));
    i = 1 as libc::c_int as R_xlen_t;
    while i <= nm1 {
        *b.offset(i as isize) = (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *y.offset(i as isize))
            / *d.offset(i as isize)
            - *d.offset(i as isize)
                * (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
                    + 2.0f64 * *c.offset(i as isize));
        *d.offset(i as isize) = (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *c.offset(i as isize))
            / *d.offset(i as isize);
        *c.offset(i as isize) = 3.0f64 * *c.offset(i as isize);
        i += 1
    }
    *c.offset(n as isize) = 3.0f64 * *c.offset(n as isize);
    *d.offset(n as isize) = *d.offset(nm1 as isize);
}
/*
 * Periodic Spline
 * ---------------
 * The end conditions here match spline (and its derivatives)
 * at x[1] and x[n].
 *
 * Note: There is an explicit check that the user has supplied
 * data with y[1] equal to y[n].
 */
unsafe extern "C" fn periodic_spline(
    mut n: R_xlen_t,
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut c: *mut libc::c_double,
    mut d: *mut libc::c_double,
) {
    let mut e: *mut libc::c_double = R_alloc(
        n as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong as libc::c_int,
    ) as *mut libc::c_double;
    /* Adjustment for 1-based arrays */
    x = x.offset(-1);
    y = y.offset(-1);
    b = b.offset(-1);
    c = c.offset(-1);
    d = d.offset(-1);
    e = e.offset(-1);
    if n < 2 as libc::c_int as libc::c_long
        || *y.offset(1 as libc::c_int as isize) != *y.offset(n as isize)
    {
        *__errno_location() = 33 as libc::c_int;
        return;
    }
    if n == 2 as libc::c_int as libc::c_long {
        let ref mut fresh8 = *d.offset(2 as libc::c_int as isize);
        *fresh8 = 0.0f64;
        let ref mut fresh9 = *d.offset(1 as libc::c_int as isize);
        *fresh9 = *fresh8;
        let ref mut fresh10 = *c.offset(2 as libc::c_int as isize);
        *fresh10 = *fresh9;
        let ref mut fresh11 = *c.offset(1 as libc::c_int as isize);
        *fresh11 = *fresh10;
        let ref mut fresh12 = *b.offset(2 as libc::c_int as isize);
        *fresh12 = *fresh11;
        *b.offset(1 as libc::c_int as isize) = *fresh12;
        return;
    } else {
        if n == 3 as libc::c_int as libc::c_long {
            let ref mut fresh13 = *b.offset(3 as libc::c_int as isize);
            *fresh13 = -(*y.offset(1 as libc::c_int as isize)
                - *y.offset(2 as libc::c_int as isize))
                * (*x.offset(1 as libc::c_int as isize)
                    - 2 as libc::c_int as libc::c_double * *x.offset(2 as libc::c_int as isize)
                    + *x.offset(3 as libc::c_int as isize))
                / (*x.offset(3 as libc::c_int as isize) - *x.offset(2 as libc::c_int as isize))
                / (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
            let ref mut fresh14 = *b.offset(2 as libc::c_int as isize);
            *fresh14 = *fresh13;
            *b.offset(1 as libc::c_int as isize) = *fresh14;
            *c.offset(1 as libc::c_int as isize) = -(3 as libc::c_int) as libc::c_double
                * (*y.offset(1 as libc::c_int as isize) - *y.offset(2 as libc::c_int as isize))
                / (*x.offset(3 as libc::c_int as isize) - *x.offset(2 as libc::c_int as isize))
                / (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
            *c.offset(2 as libc::c_int as isize) = -*c.offset(1 as libc::c_int as isize);
            *c.offset(3 as libc::c_int as isize) = *c.offset(1 as libc::c_int as isize);
            *d.offset(1 as libc::c_int as isize) = -(2 as libc::c_int) as libc::c_double
                * *c.offset(1 as libc::c_int as isize)
                / 3 as libc::c_int as libc::c_double
                / (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize));
            *d.offset(2 as libc::c_int as isize) = -*d.offset(1 as libc::c_int as isize)
                * (*x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize))
                / (*x.offset(3 as libc::c_int as isize) - *x.offset(2 as libc::c_int as isize));
            *d.offset(3 as libc::c_int as isize) = *d.offset(1 as libc::c_int as isize);
            return;
        }
    }
    /* else --------- n >= 4 --------- */
    let mut s: libc::c_double = 0.;
    let nm1: R_xlen_t = n - 1 as libc::c_int as libc::c_long;
    let mut i: R_xlen_t = 0;
    /* Set up the matrix system */
    /* A = diagonal  B = off-diagonal  C = rhs */
    *d.offset(1 as libc::c_int as isize) =
        *x.offset(2 as libc::c_int as isize) - *x.offset(1 as libc::c_int as isize);
    *d.offset(nm1 as isize) = *x.offset(n as isize) - *x.offset(nm1 as isize);
    *b.offset(1 as libc::c_int as isize) =
        2.0f64 * (*d.offset(1 as libc::c_int as isize) + *d.offset(nm1 as isize));
    *c.offset(1 as libc::c_int as isize) = (*y.offset(2 as libc::c_int as isize)
        - *y.offset(1 as libc::c_int as isize))
        / *d.offset(1 as libc::c_int as isize)
        - (*y.offset(n as isize) - *y.offset(nm1 as isize)) / *d.offset(nm1 as isize);
    i = 2 as libc::c_int as R_xlen_t;
    while i < n {
        *d.offset(i as isize) =
            *x.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *x.offset(i as isize);
        *b.offset(i as isize) = 2.0f64
            * (*d.offset(i as isize) + *d.offset((i - 1 as libc::c_int as libc::c_long) as isize));
        *c.offset(i as isize) = (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *y.offset(i as isize))
            / *d.offset(i as isize)
            - (*y.offset(i as isize) - *y.offset((i - 1 as libc::c_int as libc::c_long) as isize))
                / *d.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        i += 1
    }
    /* Choleski decomposition */
    *b.offset(1 as libc::c_int as isize) = sqrt(*b.offset(1 as libc::c_int as isize));
    *e.offset(1 as libc::c_int as isize) =
        (*x.offset(n as isize) - *x.offset(nm1 as isize)) / *b.offset(1 as libc::c_int as isize);
    s = 0.0f64;
    i = 1 as libc::c_int as R_xlen_t;
    while i <= nm1 - 2 as libc::c_int as libc::c_long {
        *d.offset(i as isize) = *d.offset(i as isize) / *b.offset(i as isize);
        if i != 1 as libc::c_int as libc::c_long {
            *e.offset(i as isize) = -*e.offset((i - 1 as libc::c_int as libc::c_long) as isize)
                * *d.offset((i - 1 as libc::c_int as libc::c_long) as isize)
                / *b.offset(i as isize)
        }
        *b.offset((i + 1 as libc::c_int as libc::c_long) as isize) = sqrt(
            *b.offset((i + 1 as libc::c_int as libc::c_long) as isize)
                - *d.offset(i as isize) * *d.offset(i as isize),
        );
        s = s + *e.offset(i as isize) * *e.offset(i as isize);
        i += 1
    }
    *d.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize) = (*d
        .offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
        - *e.offset((nm1 - 2 as libc::c_int as libc::c_long) as isize)
            * *d.offset((nm1 - 2 as libc::c_int as libc::c_long) as isize))
        / *b.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize);
    *b.offset(nm1 as isize) = sqrt(
        *b.offset(nm1 as isize)
            - *d.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
                * *d.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
            - s,
    );
    /* Forward Elimination */
    *c.offset(1 as libc::c_int as isize) =
        *c.offset(1 as libc::c_int as isize) / *b.offset(1 as libc::c_int as isize);
    s = 0.0f64;
    i = 2 as libc::c_int as R_xlen_t;
    while i <= nm1 - 1 as libc::c_int as libc::c_long {
        *c.offset(i as isize) = (*c.offset(i as isize)
            - *d.offset((i - 1 as libc::c_int as libc::c_long) as isize)
                * *c.offset((i - 1 as libc::c_int as libc::c_long) as isize))
            / *b.offset(i as isize);
        s = s + *e.offset((i - 1 as libc::c_int as libc::c_long) as isize)
            * *c.offset((i - 1 as libc::c_int as libc::c_long) as isize);
        i += 1
    }
    *c.offset(nm1 as isize) = (*c.offset(nm1 as isize)
        - *d.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
            * *c.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
        - s)
        / *b.offset(nm1 as isize);
    *c.offset(nm1 as isize) = *c.offset(nm1 as isize) / *b.offset(nm1 as isize);
    *c.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize) = (*c
        .offset((nm1 - 1 as libc::c_int as libc::c_long) as isize)
        - *d.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize) * *c.offset(nm1 as isize))
        / *b.offset((nm1 - 1 as libc::c_int as libc::c_long) as isize);
    i = nm1 - 2 as libc::c_int as libc::c_long;
    while i >= 1 as libc::c_int as libc::c_long {
        *c.offset(i as isize) = (*c.offset(i as isize)
            - *d.offset(i as isize) * *c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *e.offset(i as isize) * *c.offset(nm1 as isize))
            / *b.offset(i as isize);
        i -= 1
    }
    /* Wrap around */
    *c.offset(n as isize) = *c.offset(1 as libc::c_int as isize);
    /* Compute polynomial coefficients */
    i = 1 as libc::c_int as R_xlen_t;
    while i <= nm1 {
        s = *x.offset((i + 1 as libc::c_int as libc::c_long) as isize) - *x.offset(i as isize);
        *b.offset(i as isize) = (*y.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *y.offset(i as isize))
            / s
            - s * (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
                + 2.0f64 * *c.offset(i as isize));
        *d.offset(i as isize) = (*c.offset((i + 1 as libc::c_int as libc::c_long) as isize)
            - *c.offset(i as isize))
            / s;
        *c.offset(i as isize) = 3.0f64 * *c.offset(i as isize);
        i += 1
    }
    *b.offset(n as isize) = *b.offset(1 as libc::c_int as isize);
    *c.offset(n as isize) = *c.offset(1 as libc::c_int as isize);
    *d.offset(n as isize) = *d.offset(1 as libc::c_int as isize);
}
/* These were/are the public interfaces */
unsafe extern "C" fn spline_coef(
    mut method: libc::c_int,
    mut n: R_xlen_t,
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut c: *mut libc::c_double,
    mut d: *mut libc::c_double,
) {
    match method {
        1 => {
            periodic_spline(n, x, y, b, c, d);
        }
        2 => {
            natural_spline(n, x, y, b, c, d);
        }
        3 => {
            fmm_spline(n, x, y, b, c, d);
        }
        _ => {}
    };
}
#[no_mangle]
pub unsafe extern "C" fn SplineCoef(mut method: SEXP, mut x: SEXP, mut y: SEXP) -> SEXP {
    x = protect(coerceVector(x, 14 as libc::c_int as SEXPTYPE));
    y = protect(coerceVector(y, 14 as libc::c_int as SEXPTYPE));
    let mut n: R_xlen_t = XLENGTH(x);
    let mut m: libc::c_int = asInteger(method);
    if XLENGTH(y) != n {
        error(b"inputs of different lengths\x00" as *const u8 as *const libc::c_char);
    }
    let mut b: SEXP = 0 as *mut SEXPREC;
    let mut c: SEXP = 0 as *mut SEXPREC;
    let mut d: SEXP = 0 as *mut SEXPREC;
    let mut ans: SEXP = 0 as *mut SEXPREC;
    let mut nm: SEXP = 0 as *mut SEXPREC;
    b = protect(allocVector(14 as libc::c_int as SEXPTYPE, n));
    c = protect(allocVector(14 as libc::c_int as SEXPTYPE, n));
    d = protect(allocVector(14 as libc::c_int as SEXPTYPE, n));
    let mut rb: *mut libc::c_double = REAL(b);
    let mut rc: *mut libc::c_double = REAL(c);
    let mut rd: *mut libc::c_double = REAL(d);
    let mut i: R_xlen_t = 0 as libc::c_int as R_xlen_t;
    while i < n {
        let ref mut fresh15 = *rd.offset(i as isize);
        *fresh15 = 0 as libc::c_int as libc::c_double;
        let ref mut fresh16 = *rc.offset(i as isize);
        *fresh16 = *fresh15;
        *rb.offset(i as isize) = *fresh16;
        i += 1
    }
    spline_coef(m, n, REAL(x), REAL(y), rb, rc, rd);
    ans = protect(allocVector(
        19 as libc::c_int as SEXPTYPE,
        7 as libc::c_int as R_xlen_t,
    ));
    SET_VECTOR_ELT(ans, 0 as libc::c_int as R_xlen_t, ScalarInteger(m));
    SET_VECTOR_ELT(
        ans,
        1 as libc::c_int as R_xlen_t,
        if n > 2147483647 as libc::c_int as libc::c_long {
            ScalarReal(n as libc::c_double)
        } else {
            ScalarInteger(n as libc::c_int)
        },
    );
    SET_VECTOR_ELT(ans, 2 as libc::c_int as R_xlen_t, x);
    SET_VECTOR_ELT(ans, 3 as libc::c_int as R_xlen_t, y);
    SET_VECTOR_ELT(ans, 4 as libc::c_int as R_xlen_t, b);
    SET_VECTOR_ELT(ans, 5 as libc::c_int as R_xlen_t, c);
    SET_VECTOR_ELT(ans, 6 as libc::c_int as R_xlen_t, d);
    nm = allocVector(16 as libc::c_int as SEXPTYPE, 7 as libc::c_int as R_xlen_t);
    setAttrib(ans, R_NamesSymbol, nm);
    SET_STRING_ELT(
        nm,
        0 as libc::c_int as R_xlen_t,
        mkChar(b"method\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        1 as libc::c_int as R_xlen_t,
        mkChar(b"n\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        2 as libc::c_int as R_xlen_t,
        mkChar(b"x\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        3 as libc::c_int as R_xlen_t,
        mkChar(b"y\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        4 as libc::c_int as R_xlen_t,
        mkChar(b"b\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        5 as libc::c_int as R_xlen_t,
        mkChar(b"c\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        6 as libc::c_int as R_xlen_t,
        mkChar(b"d\x00" as *const u8 as *const libc::c_char),
    );
    unprotect(6 as libc::c_int);
    return ans;
}
unsafe extern "C" fn spline_eval(
    mut method: libc::c_int,
    mut nu: R_xlen_t,
    mut u: *mut libc::c_double,
    mut v: *mut libc::c_double,
    mut n: R_xlen_t,
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut c: *mut libc::c_double,
    mut d: *mut libc::c_double,
) {
    /* Evaluate  v[l] := spline(u[l], ...),     l = 1,..,nu, i.e. 0:(nu-1)
     * Nodes x[i], coef (y[i]; b[i],c[i],d[i]); i = 1,..,n , i.e. 0:(*n-1)
     */
    let n_1: R_xlen_t = n - 1 as libc::c_int as libc::c_long;
    let mut i: R_xlen_t = 0;
    let mut l: R_xlen_t = 0;
    let mut dx: libc::c_double = 0.;
    if method == 1 as libc::c_int && n > 1 as libc::c_int as libc::c_long {
        /* periodic */
        dx = *x.offset(n_1 as isize) - *x.offset(0 as libc::c_int as isize);
        l = 0 as libc::c_int as R_xlen_t;
        while l < nu {
            *v.offset(l as isize) =
                *u.offset(l as isize) - *x.offset(0 as libc::c_int as isize) % dx;
            if *v.offset(l as isize) < 0.0f64 {
                *v.offset(l as isize) += dx
            }
            *v.offset(l as isize) += *x.offset(0 as libc::c_int as isize);
            l += 1
        }
    } else {
        l = 0 as libc::c_int as R_xlen_t;
        while l < nu {
            *v.offset(l as isize) = *u.offset(l as isize);
            l += 1
        }
    }
    l = 0 as libc::c_int as R_xlen_t;
    i = 0 as libc::c_int as R_xlen_t;
    while l < nu {
        let mut ul: libc::c_double = *v.offset(l as isize);
        if ul < *x.offset(i as isize)
            || i < n_1 && *x.offset((i + 1 as libc::c_int as libc::c_long) as isize) < ul
        {
            /* reset i  such that  x[i] <= ul <= x[i+1] : */
            i = 0 as libc::c_int as R_xlen_t;
            let mut j: R_xlen_t = n;
            loop {
                let mut k: R_xlen_t = (i + j) / 2 as libc::c_int as libc::c_long;
                if ul < *x.offset(k as isize) {
                    j = k
                } else {
                    i = k
                }
                if !(j > i + 1 as libc::c_int as libc::c_long) {
                    break;
                }
            }
        }
        dx = ul - *x.offset(i as isize);
        /* for natural splines extrapolate linearly left */
        let mut tmp: libc::c_double =
            if method == 2 as libc::c_int && ul < *x.offset(0 as libc::c_int as isize) {
                0.0f64
            } else {
                *d.offset(i as isize)
            };
        *v.offset(l as isize) = *y.offset(i as isize)
            + dx * (*b.offset(i as isize) + dx * (*c.offset(i as isize) + dx * tmp));
        l += 1
    }
}
// TODO: move to ../../../main/coerce.c
unsafe extern "C" fn asXlen(mut x: SEXP) -> R_xlen_t {
    if isVectorAtomic(x) as libc::c_uint != 0 && XLENGTH(x) >= 1 as libc::c_int as libc::c_long {
        match TYPEOF(x) {
            13 => return *INTEGER(x).offset(0 as libc::c_int as isize) as R_xlen_t,
            14 => return *REAL(x).offset(0 as libc::c_int as isize) as R_xlen_t,
            _ => {
                UNIMPLEMENTED_TYPE(b"asXlen\x00" as *const u8 as *const libc::c_char, x);
            }
        }
    }
    return R_NaInt as R_xlen_t;
}
#[no_mangle]
pub unsafe extern "C" fn SplineEval(mut xout: SEXP, mut z: SEXP) -> SEXP {
    xout = protect(coerceVector(xout, 14 as libc::c_int as SEXPTYPE));
    let mut nu: R_xlen_t = XLENGTH(xout);
    let mut nx: R_xlen_t = asXlen(getListElement(
        z,
        b"n\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    ));
    let mut yout: SEXP = protect(allocVector(14 as libc::c_int as SEXPTYPE, nu));
    let mut method: libc::c_int = asInteger(getListElement(
        z,
        b"method\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    ));
    let mut x: SEXP = getListElement(
        z,
        b"x\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    );
    let mut y: SEXP = getListElement(
        z,
        b"y\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    );
    let mut b: SEXP = getListElement(
        z,
        b"b\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    );
    let mut c: SEXP = getListElement(
        z,
        b"c\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    );
    let mut d: SEXP = getListElement(
        z,
        b"d\x00" as *const u8 as *const libc::c_char as *mut libc::c_char,
    );
    spline_eval(
        method,
        nu,
        REAL(xout),
        REAL(yout),
        nx,
        REAL(x),
        REAL(y),
        REAL(b),
        REAL(c),
        REAL(d),
    );
    unprotect(2 as libc::c_int);
    return yout;
}