r2rs-stats 0.1.1

use ::libc;
extern "C" {

    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;
    /* General Support Functions */
    #[no_mangle]
    fn imax2(_: libc::c_int, _: libc::c_int) -> libc::c_int;
    #[no_mangle]
    fn imin2(_: libc::c_int, _: libc::c_int) -> libc::c_int;
}
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;
/*, MAYBE */
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995, 1996, 1997  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1998--2000, 2013  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/
 */
/* for INT_MAX */
/* for size_t */
/* for abs */
/*  Fast Fourier Transform
 *
 *  These routines are based on code by Richard Singleton in the
 *  book "Programs for Digital Signal Processing" put out by IEEE.
 *
 *  I have translated them to C and moved the memory allocation
 *  so that it takes place under the control of the algorithm
 *  which calls these; for R, see ../main/fourier.c
 *
 *  void fft_factor(int n, int *maxf, int *maxp)
 *
 * This factorizes the series length and computes the values of
 * maxf and maxp which determine the amount of scratch storage
 * required by the algorithm.
 *
 * If maxf is zero on return, an error occured during factorization.
 * The nature of the error can be determined from the value of maxp.
 * If maxp is zero, an invalid (zero) parameter was passed and
 * if maxp is one,  the internal nfac array was too small.  This can only
 * happen for series lengths which exceed 12,754,584.
 *
 * PROBLEM (see fftmx  below): nfac[] is overwritten by fftmx() in fft_work()
 * -------  Consequence:  fft_factor() must be called way too often,
 * at least from  do_mvfft() [ ../main/fourier.c ]
 *
 * The following arrays need to be allocated following the call to
 * fft_factor and preceding the call to fft_work.
 *
 *  work double[4*maxf]
 *  iwork int[maxp]
 *
 *  int fft_work(double *a, double *b, int nseg, int n, int nspn,
 *   int isn, double *work, int *iwork)
 *
 * The routine returns 1 if the transform was completed successfully and
 *       0 if invalid values of the parameters were supplied.
 *
 *  Ross Ihaka
 *  University of Auckland
 *  February 1997
 *  ==========================================================================
 *
 *  Header from the original Singleton algorithm:
 *
 * --------------------------------------------------------------
 * subroutine:  fft
 * multivariate complex fourier transform, computed in place
 * using mixed-radix fast fourier transform algorithm.
 * --------------------------------------------------------------
 *
 * arrays a and b originally hold the real and imaginary
 * components of the data, and return the real and
 * imaginary components of the resulting fourier coefficients.
 * multivariate data is indexed according to the fortran
 * array element successor function, without limit
 * on the number of implied multiple subscripts.
 * the subroutine is called once for each variate.
 * the calls for a multivariate transform may be in any order.
 *
 * n is the dimension of the current variable.
 * nspn is the spacing of consecutive data values
 * while indexing the current variable.
 * nseg nseg*n*nspn is the total number of complex data values.
 * isn the sign of isn determines the sign of the complex
 * exponential, and the magnitude of isn is normally one.
 * the magnitude of isn determines the indexing increment for a&b.
 *
 * if fft is called twice, with opposite signs on isn, an
 * identity transformation is done...calls can be in either order.
 * the results are scaled by 1/n when the sign of isn is positive.
 *
 * a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
 * is computed by
 * call fft(a,b,n2*n3,n1, 1,   -1)
 * call fft(a,b,n3   ,n2,n1,   -1)
 * call fft(a,b,1,    n3,n1*n2,-1)
 *
 * a single-variate transform of n complex data values is computed by
 * call fft(a,b,1,n,1,-1)
 *
 * the data may alternatively be stored in a single complex
 * array a, then the magnitude of isn changed to two to
 * give the correct indexing increment and a(2) used to
 * pass the initial address for the sequence of imaginary
 * values, e.g.
 *    call fft(a,a(2),nseg,n,nspn,-2)
 *
 * nfac[15] (array) is working storage for factoring n.  the smallest
 * number exceeding the 15 locations provided is 12,754,584.
 *
 * Update in R 3.1.0: nfac[20], increased array size. It is now possible to
 * factor any positive int n, up to 2^31 - 1.
 */
unsafe extern "C" fn fftmx(
    mut a: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut ntot: libc::c_int,
    mut n: libc::c_int,
    mut nspan: libc::c_int,
    mut isn: libc::c_int,
    mut m: libc::c_int,
    mut kt_0: libc::c_int,
    mut at: *mut libc::c_double,
    mut ck: *mut libc::c_double,
    mut bt: *mut libc::c_double,
    mut sk: *mut libc::c_double,
    mut np: *mut libc::c_int,
    mut nfac_0: *mut libc::c_int,
) {
    let mut current_block: u64;
    /* called from fft_work() */
    /* Design BUG: One purpose of fft_factor() would be to compute
     * ---------- nfac[] once and for all; and fft_work() [i.e. fftmx ]
     *  could reuse the factorization.
     * However: nfac[] is `destroyed' currently in the code below
     */
    let mut aa: libc::c_double = 0.; /*the global one!*/
    let mut aj: libc::c_double = 0.; /* = pi/4 =^= 45 degrees */
    let mut ajm: libc::c_double = 0.; /* 72 = 45 / .625  degrees */
    let mut ajp: libc::c_double = 0.; /* 120.sin() = sqrt(3)/2 */
    let mut ak: libc::c_double = 0.;
    let mut akm: libc::c_double = 0.;
    let mut akp: libc::c_double = 0.;
    let mut bb: libc::c_double = 0.;
    let mut bj: libc::c_double = 0.;
    let mut bjm: libc::c_double = 0.;
    let mut bjp: libc::c_double = 0.;
    let mut bk: libc::c_double = 0.;
    let mut bkm: libc::c_double = 0.;
    let mut bkp: libc::c_double = 0.;
    let mut c1: libc::c_double = 0.;
    let mut c2: libc::c_double = 0 as libc::c_int as libc::c_double;
    let mut c3: libc::c_double = 0 as libc::c_int as libc::c_double;
    let mut c72: libc::c_double = 0.;
    let mut cd: libc::c_double = 0.;
    let mut dr: libc::c_double = 0.;
    let mut rad: libc::c_double = 0.;
    let mut s1: libc::c_double = 0.;
    let mut s120: libc::c_double = 0.;
    let mut s2: libc::c_double = 0 as libc::c_int as libc::c_double;
    let mut s3: libc::c_double = 0 as libc::c_int as libc::c_double;
    let mut s72: libc::c_double = 0.;
    let mut sd: libc::c_double = 0.;
    let mut i: libc::c_int = 0;
    let mut inc: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut jc: libc::c_int = 0;
    let mut jf: libc::c_int = 0;
    let mut jj: libc::c_int = 0;
    let mut k: libc::c_int = 0;
    let mut k1: libc::c_int = 0;
    let mut k2: libc::c_int = 0;
    let mut k3: libc::c_int = 0 as libc::c_int;
    let mut k4: libc::c_int = 0;
    let mut kk: libc::c_int = 0;
    let mut klim: libc::c_int = 0;
    let mut ks: libc::c_int = 0;
    let mut kspan: libc::c_int = 0;
    let mut kspnn: libc::c_int = 0;
    let mut lim: libc::c_int = 0;
    let mut maxf_0: libc::c_int = 0;
    let mut mm: libc::c_int = 0;
    let mut nn: libc::c_int = 0;
    let mut nt: libc::c_int = 0;
    a = a.offset(-1);
    b = b.offset(-1);
    at = at.offset(-1);
    ck = ck.offset(-1);
    bt = bt.offset(-1);
    sk = sk.offset(-1);
    np = np.offset(-1);
    nfac_0 = nfac_0.offset(-1);
    inc = isn.abs();
    nt = inc * ntot;
    ks = inc * nspan;
    rad = 0.78539816339744830962f64;
    s72 = rad / 0.625f64;
    c72 = s72.cos();
    s72 = s72.sin();
    s120 = 0.5f64 * 1.732050807568877293527446341506f64;
    if isn <= 0 as libc::c_int {
        s72 = -s72;
        s120 = -s120;
        rad = -rad
    }
    kspan = ks;
    nn = nt - inc;
    jc = ks / n;
    /* sin, cos values are re-initialized each lim steps */
    lim = 32 as libc::c_int;
    klim = lim * jc;
    i = 0 as libc::c_int;
    jf = 0 as libc::c_int;
    maxf_0 = *nfac_0.offset((m - kt_0) as isize);
    if kt_0 > 0 as libc::c_int {
        maxf_0 = imax2(*nfac_0.offset(kt_0 as isize), maxf_0)
    }
    loop
    /* compute fourier transform */
    {
        dr = 8.0f64 * jc as libc::c_double / kspan as libc::c_double;
        cd = (0.5f64 * dr * rad).sin();
        cd = 2.0f64 * cd * cd;
        sd = (dr * rad).sin();
        kk = 1 as libc::c_int;
        i += 1;
        if *nfac_0.offset(i as isize) != 2 as libc::c_int {
            /* transform for factor of 4 */
            if *nfac_0.offset(i as isize) != 4 as libc::c_int {
                /* transform for odd factors */
                k = *nfac_0.offset(i as isize);
                kspnn = kspan;
                kspan /= k;
                if k == 3 as libc::c_int {
                    loop
                    /* transform for factor of 3 (optional code) */
                    {
                        k1 = kk + kspan;
                        k2 = k1 + kspan;
                        ak = *a.offset(kk as isize);
                        bk = *b.offset(kk as isize);
                        aj = *a.offset(k1 as isize) + *a.offset(k2 as isize);
                        bj = *b.offset(k1 as isize) + *b.offset(k2 as isize);
                        *a.offset(kk as isize) = ak + aj;
                        *b.offset(kk as isize) = bk + bj;
                        ak = -0.5f64 * aj + ak;
                        bk = -0.5f64 * bj + bk;
                        aj = (*a.offset(k1 as isize) - *a.offset(k2 as isize)) * s120;
                        bj = (*b.offset(k1 as isize) - *b.offset(k2 as isize)) * s120;
                        *a.offset(k1 as isize) = ak - bj;
                        *b.offset(k1 as isize) = bk + aj;
                        *a.offset(k2 as isize) = ak + bj;
                        *b.offset(k2 as isize) = bk - aj;
                        kk = k2 + kspan;
                        if kk < nn {
                            continue;
                        }
                        kk = kk - nn;
                        if !(kk <= kspan) {
                            break;
                        }
                    }
                } else if k == 5 as libc::c_int {
                    /* transform for factor of 5 (optional code) */
                    c2 = c72 * c72 - s72 * s72;
                    s2 = 2.0f64 * c72 * s72;
                    loop {
                        k1 = kk + kspan;
                        k2 = k1 + kspan;
                        k3 = k2 + kspan;
                        k4 = k3 + kspan;
                        akp = *a.offset(k1 as isize) + *a.offset(k4 as isize);
                        akm = *a.offset(k1 as isize) - *a.offset(k4 as isize);
                        bkp = *b.offset(k1 as isize) + *b.offset(k4 as isize);
                        bkm = *b.offset(k1 as isize) - *b.offset(k4 as isize);
                        ajp = *a.offset(k2 as isize) + *a.offset(k3 as isize);
                        ajm = *a.offset(k2 as isize) - *a.offset(k3 as isize);
                        bjp = *b.offset(k2 as isize) + *b.offset(k3 as isize);
                        bjm = *b.offset(k2 as isize) - *b.offset(k3 as isize);
                        aa = *a.offset(kk as isize);
                        bb = *b.offset(kk as isize);
                        *a.offset(kk as isize) = aa + akp + ajp;
                        *b.offset(kk as isize) = bb + bkp + bjp;
                        ak = akp * c72 + ajp * c2 + aa;
                        bk = bkp * c72 + bjp * c2 + bb;
                        aj = akm * s72 + ajm * s2;
                        bj = bkm * s72 + bjm * s2;
                        *a.offset(k1 as isize) = ak - bj;
                        *a.offset(k4 as isize) = ak + bj;
                        *b.offset(k1 as isize) = bk + aj;
                        *b.offset(k4 as isize) = bk - aj;
                        ak = akp * c2 + ajp * c72 + aa;
                        bk = bkp * c2 + bjp * c72 + bb;
                        aj = akm * s2 - ajm * s72;
                        bj = bkm * s2 - bjm * s72;
                        *a.offset(k2 as isize) = ak - bj;
                        *a.offset(k3 as isize) = ak + bj;
                        *b.offset(k2 as isize) = bk + aj;
                        *b.offset(k3 as isize) = bk - aj;
                        kk = k4 + kspan;
                        if kk < nn {
                            continue;
                        }
                        kk = kk - nn;
                        if !(kk <= kspan) {
                            break;
                        }
                    }
                } else {
                    if !(k == jf) {
                        jf = k;
                        s1 = rad / (k as libc::c_double / 8.0f64);
                        c1 = s1.cos();
                        s1 = s1.sin();
                        *ck.offset(jf as isize) = 1.0f64;
                        *sk.offset(jf as isize) = 0.0f64;
                        j = 1 as libc::c_int;
                        while j < k {
                            /* k is changing as well */
                            *ck.offset(j as isize) =
                                *ck.offset(k as isize) * c1 + *sk.offset(k as isize) * s1;
                            *sk.offset(j as isize) =
                                *ck.offset(k as isize) * s1 - *sk.offset(k as isize) * c1;
                            k -= 1;
                            *ck.offset(k as isize) = *ck.offset(j as isize);
                            *sk.offset(k as isize) = -*sk.offset(j as isize);
                            j += 1
                        }
                    }
                    loop {
                        k1 = kk;
                        k2 = kk + kspnn;
                        aa = *a.offset(kk as isize);
                        bb = *b.offset(kk as isize);
                        ak = aa;
                        bk = bb;
                        j = 1 as libc::c_int;
                        k1 = k1 + kspan;
                        loop {
                            k2 = k2 - kspan;
                            j += 1;
                            *at.offset(j as isize) =
                                *a.offset(k1 as isize) + *a.offset(k2 as isize);
                            ak = *at.offset(j as isize) + ak;
                            *bt.offset(j as isize) =
                                *b.offset(k1 as isize) + *b.offset(k2 as isize);
                            bk = *bt.offset(j as isize) + bk;
                            j += 1;
                            *at.offset(j as isize) =
                                *a.offset(k1 as isize) - *a.offset(k2 as isize);
                            *bt.offset(j as isize) =
                                *b.offset(k1 as isize) - *b.offset(k2 as isize);
                            k1 = k1 + kspan;
                            if !(k1 < k2) {
                                break;
                            }
                        }
                        *a.offset(kk as isize) = ak;
                        *b.offset(kk as isize) = bk;
                        k1 = kk;
                        k2 = kk + kspnn;
                        j = 1 as libc::c_int;
                        loop {
                            k1 += kspan;
                            k2 -= kspan;
                            jj = j;
                            ak = aa;
                            bk = bb;
                            aj = 0.0f64;
                            bj = 0.0f64;
                            k = 1 as libc::c_int;
                            k = 2 as libc::c_int;
                            while k < jf {
                                ak += *at.offset(k as isize) * *ck.offset(jj as isize);
                                bk += *bt.offset(k as isize) * *ck.offset(jj as isize);
                                k += 1;
                                aj += *at.offset(k as isize) * *sk.offset(jj as isize);
                                bj += *bt.offset(k as isize) * *sk.offset(jj as isize);
                                jj += j;
                                if jj > jf {
                                    jj -= jf
                                }
                                k += 1
                            }
                            k = jf - j;
                            *a.offset(k1 as isize) = ak - bj;
                            *b.offset(k1 as isize) = bk + aj;
                            *a.offset(k2 as isize) = ak + bj;
                            *b.offset(k2 as isize) = bk - aj;
                            j += 1;
                            if !(j < k) {
                                break;
                            }
                        }
                        kk = kk + kspnn;
                        if kk <= nn {
                            continue;
                        }
                        kk = kk - nn;
                        if !(kk <= kspan) {
                            break;
                        }
                    }
                }
                /* multiply by rotation factor (except for factors of 2 and 4) */
                if i == m {
                    break;
                }
                kk = jc + 1 as libc::c_int;
                'c_6920: loop {
                    c2 = 1.0f64 - cd;
                    s1 = sd;
                    mm = imin2(kspan, klim);
                    loop
                    /* the following three statements compensate for truncation error.*/
                    /* if rounded arithmetic is used (nowadays always ?!), they may be deleted. */
                    /* L320: */
                    {
                        c1 = c2;
                        s2 = s1;
                        kk += kspan;
                        loop {
                            loop
                            /* L330: */
                            {
                                ak = *a.offset(kk as isize);
                                *a.offset(kk as isize) = c2 * ak - s2 * *b.offset(kk as isize);
                                *b.offset(kk as isize) = s2 * ak + c2 * *b.offset(kk as isize);
                                kk += kspnn;
                                if !(kk <= nt) {
                                    break;
                                }
                            }
                            ak = s1 * s2;
                            s2 = s1 * c2 + c1 * s2;
                            c2 = c1 * c2 - ak;
                            kk += -nt + kspan;
                            if !(kk <= kspnn) {
                                break;
                            }
                        }
                        kk += -kspnn + jc;
                        if kk <= mm {
                            /* L310: */
                            c2 = c1 - (cd * c1 + sd * s1);
                            s1 = s1 + (sd * c1 - cd * s1)
                        /* goto L320*/
                        } else if kk >= kspan {
                            kk = kk - kspan + jc + inc;
                            if kk <= jc + jc {
                                break;
                            } else {
                                break 'c_6920;
                            }
                        } else {
                            s1 = ((kk - 1 as libc::c_int) / jc) as libc::c_double * dr * rad;
                            c2 = s1.cos();
                            s1 = s1.sin();
                            mm = imin2(kspan, mm + klim)
                        }
                    }
                }
            } else {
                kspnn = kspan;
                kspan /= 4 as libc::c_int;
                loop {
                    c1 = 1.0f64;
                    s1 = 0 as libc::c_int as libc::c_double;
                    mm = imin2(kspan, klim);
                    loop {
                        k1 = kk + kspan;
                        k2 = k1 + kspan;
                        k3 = k2 + kspan;
                        akp = *a.offset(kk as isize) + *a.offset(k2 as isize);
                        akm = *a.offset(kk as isize) - *a.offset(k2 as isize);
                        ajp = *a.offset(k1 as isize) + *a.offset(k3 as isize);
                        ajm = *a.offset(k1 as isize) - *a.offset(k3 as isize);
                        *a.offset(kk as isize) = akp + ajp;
                        ajp = akp - ajp;
                        bkp = *b.offset(kk as isize) + *b.offset(k2 as isize);
                        bkm = *b.offset(kk as isize) - *b.offset(k2 as isize);
                        bjp = *b.offset(k1 as isize) + *b.offset(k3 as isize);
                        bjm = *b.offset(k1 as isize) - *b.offset(k3 as isize);
                        *b.offset(kk as isize) = bkp + bjp;
                        bjp = bkp - bjp;
                        if isn < 0 as libc::c_int {
                            akp = akm + bjm;
                            akm = akm - bjm;
                            bkp = bkm - ajm;
                            bkm = bkm + ajm;
                            if s1 != 0.0f64 {
                                current_block = 8794591204701163763;
                            } else {
                                current_block = 9725544835206080425;
                            }
                        } else {
                            akp = akm - bjm;
                            akm = akm + bjm;
                            bkp = bkm + ajm;
                            bkm = bkm - ajm;
                            if s1 == 0.0f64 {
                                current_block = 9725544835206080425;
                            } else {
                                current_block = 8794591204701163763;
                            }
                        }
                        match current_block {
                            9725544835206080425 => {
                                *a.offset(k1 as isize) = akp;
                                *b.offset(k1 as isize) = bkp;
                                *a.offset(k2 as isize) = ajp;
                                *b.offset(k2 as isize) = bjp;
                                *a.offset(k3 as isize) = akm;
                                *b.offset(k3 as isize) = bkm;
                                kk = k3 + kspan;
                                if kk <= nt {
                                    continue;
                                }
                            }
                            _ => {
                                *a.offset(k1 as isize) = akp * c1 - bkp * s1;
                                *b.offset(k1 as isize) = akp * s1 + bkp * c1;
                                *a.offset(k2 as isize) = ajp * c2 - bjp * s2;
                                *b.offset(k2 as isize) = ajp * s2 + bjp * c2;
                                *a.offset(k3 as isize) = akm * c3 - bkm * s3;
                                *b.offset(k3 as isize) = akm * s3 + bkm * c3;
                                kk = k3 + kspan;
                                if kk <= nt {
                                    continue;
                                }
                            }
                        }
                        kk = kk - nt + jc;
                        if kk <= mm {
                            c2 = c1 - (cd * c1 + sd * s1);
                            s1 = sd * c1 - cd * s1 + s1;
                            /* the following three statements compensate for truncation error. */
                            /* if rounded arithmetic is used (nowadays always ?!), substitute  c1=c2 */
                            c1 = c2
                        } else {
                            if !(kk < kspan) {
                                break;
                            }
                            s1 = ((kk - 1 as libc::c_int) / jc) as libc::c_double * dr * rad;
                            c1 = s1.cos();
                            s1 = s1.sin();
                            mm = imin2(kspan, mm + klim)
                        }
                        c2 = c1 * c1 - s1 * s1;
                        s2 = c1 * s1 * 2.0f64;
                        c3 = c2 * c1 - s2 * s1;
                        s3 = c2 * s1 + s2 * c1
                    }
                    kk = kk - kspan + inc;
                    if !(kk <= jc) {
                        break;
                    }
                }
                if kspan == jc {
                    break;
                }
            }
        } else {
            /* transform for factor of 2 (including rotation factor) */
            kspan /= 2 as libc::c_int;
            k1 = kspan + 2 as libc::c_int;
            loop {
                loop {
                    k2 = kk + kspan;
                    ak = *a.offset(k2 as isize);
                    bk = *b.offset(k2 as isize);
                    *a.offset(k2 as isize) = *a.offset(kk as isize) - ak;
                    *b.offset(k2 as isize) = *b.offset(kk as isize) - bk;
                    *a.offset(kk as isize) += ak;
                    *b.offset(kk as isize) += bk;
                    kk = k2 + kspan;
                    if !(kk <= nn) {
                        break;
                    }
                }
                kk -= nn;
                if !(kk <= jc) {
                    break;
                }
            }
            if kk > kspan {
                break;
            }
            loop {
                c1 = 1.0f64 - cd;
                s1 = sd;
                mm = imin2(k1 / 2 as libc::c_int, klim);
                loop {
                    loop {
                        k2 = kk + kspan;
                        ak = *a.offset(kk as isize) - *a.offset(k2 as isize);
                        bk = *b.offset(kk as isize) - *b.offset(k2 as isize);
                        *a.offset(kk as isize) += *a.offset(k2 as isize);
                        *b.offset(kk as isize) += *b.offset(k2 as isize);
                        *a.offset(k2 as isize) = c1 * ak - s1 * bk;
                        *b.offset(k2 as isize) = s1 * ak + c1 * bk;
                        kk = k2 + kspan;
                        if !(kk < nt) {
                            break;
                        }
                    }
                    k2 = kk - nt;
                    c1 = -c1;
                    kk = k1 - k2;
                    if kk > k2 {
                        continue;
                    }
                    kk += jc;
                    if kk <= mm {
                        ak = c1 - (cd * c1 + sd * s1);
                        s1 = sd * c1 - cd * s1 + s1;
                        /* the following three statements compensate for truncation error. */
                        /* if rounded arithmetic is used (nowadays always ?!), substitute  c1=ak */
                        c1 = ak
                    } else {
                        if kk >= k2 {
                            break;
                        }
                        s1 = ((kk - 1 as libc::c_int) / jc) as libc::c_double * dr * rad;
                        c1 = s1.cos();
                        s1 = s1.sin();
                        mm = imin2(k1 / 2 as libc::c_int, mm + klim)
                    }
                }
                k1 = k1 + inc + inc;
                kk = (k1 - kspan) / 2 as libc::c_int + jc;
                if !(kk <= jc + jc) {
                    break;
                }
            }
        }
    }
    /*------------------------------------------------------------*/
    /* permute the results to normal order---done in two stages */
    /* permutation for square factors of n */
    *np.offset(1 as libc::c_int as isize) = ks;
    if !(kt_0 == 0 as libc::c_int) {
        k = kt_0 + kt_0 + 1 as libc::c_int;
        if m < k {
            k -= 1
        }
        *np.offset((k + 1 as libc::c_int) as isize) = jc;
        j = 1 as libc::c_int;
        while j < k {
            *np.offset((j + 1 as libc::c_int) as isize) =
                *np.offset(j as isize) / *nfac_0.offset(j as isize);
            *np.offset(k as isize) =
                *np.offset((k + 1 as libc::c_int) as isize) * *nfac_0.offset(j as isize);
            j += 1;
            k -= 1
        }
        k3 = *np.offset((k + 1 as libc::c_int) as isize);
        kspan = *np.offset(2 as libc::c_int as isize);
        kk = jc + 1 as libc::c_int;
        k2 = kspan + 1 as libc::c_int;
        j = 1 as libc::c_int;
        if n == ntot {
            'c_7742: loop
            /* permutation for single-variate transform (optional code) */
            {
                loop {
                    ak = *a.offset(kk as isize);
                    *a.offset(kk as isize) = *a.offset(k2 as isize);
                    *a.offset(k2 as isize) = ak;
                    bk = *b.offset(kk as isize);
                    *b.offset(kk as isize) = *b.offset(k2 as isize);
                    *b.offset(k2 as isize) = bk;
                    kk += inc;
                    k2 += kspan;
                    if !(k2 < ks) {
                        break;
                    }
                }
                'c_7743: loop {
                    loop {
                        k2 -= *np.offset(j as isize);
                        j += 1;
                        k2 += *np.offset((j + 1 as libc::c_int) as isize);
                        if !(k2 > *np.offset(j as isize)) {
                            break;
                        }
                    }
                    j = 1 as libc::c_int;
                    loop {
                        if kk < k2 {
                            break 'c_7743;
                        }
                        kk += inc;
                        k2 += kspan;
                        if !(k2 < ks) {
                            break;
                        }
                    }
                    if !(kk < ks) {
                        break 'c_7742;
                    }
                }
            }
            jc = k3
        } else {
            'c_7529: loop
            /* permutation for multivariate transform */
            {
                k = kk + jc;
                loop {
                    ak = *a.offset(kk as isize);
                    *a.offset(kk as isize) = *a.offset(k2 as isize);
                    *a.offset(k2 as isize) = ak;
                    bk = *b.offset(kk as isize);
                    *b.offset(kk as isize) = *b.offset(k2 as isize);
                    *b.offset(k2 as isize) = bk;
                    kk += inc;
                    k2 += inc;
                    if !(kk < k) {
                        break;
                    }
                }
                kk += ks - jc;
                k2 += ks - jc;
                if kk < nt {
                    continue;
                }
                k2 += -nt + kspan;
                kk += -nt + jc;
                if k2 < ks {
                    continue;
                }
                's_1717: loop {
                    loop {
                        k2 -= *np.offset(j as isize);
                        j += 1;
                        k2 += *np.offset((j + 1 as libc::c_int) as isize);
                        if !(k2 > *np.offset(j as isize)) {
                            break;
                        }
                    }
                    j = 1 as libc::c_int;
                    loop {
                        if kk < k2 {
                            break 's_1717;
                        }
                        kk += jc;
                        k2 += kspan;
                        if !(k2 < ks) {
                            break;
                        }
                    }
                    if !(kk < ks) {
                        break 'c_7529;
                    }
                }
            }
            jc = k3
        }
    }
    if 2 as libc::c_int * kt_0 + 1 as libc::c_int >= m {
        return;
    }
    kspnn = *np.offset((kt_0 + 1 as libc::c_int) as isize);
    /* permutation for square-free factors of n */
    /* Here, nfac[] is overwritten... -- now CUMULATIVE ("cumprod") factors */
    nn = m - kt_0;
    *nfac_0.offset((nn + 1 as libc::c_int) as isize) = 1 as libc::c_int;
    j = nn;
    while j > kt_0 {
        *nfac_0.offset(j as isize) *= *nfac_0.offset((j + 1 as libc::c_int) as isize);
        j -= 1
    }
    kt_0 += 1;
    nn = *nfac_0.offset(kt_0 as isize) - 1 as libc::c_int;
    jj = 0 as libc::c_int;
    j = 0 as libc::c_int;
    loop {
        k2 = *nfac_0.offset(kt_0 as isize);
        k = kt_0 + 1 as libc::c_int;
        kk = *nfac_0.offset(k as isize);
        j += 1;
        if !(j <= nn) {
            break;
        }
        loop {
            jj += kk;
            if !(jj >= k2) {
                break;
            }
            jj -= k2;
            k2 = kk;
            k += 1;
            kk = *nfac_0.offset(k as isize)
        }
        *np.offset(j as isize) = jj
    }
    /* determine the permutation cycles of length greater than 1 */
    j = 0 as libc::c_int;
    loop {
        loop {
            j += 1;
            kk = *np.offset(j as isize);
            if !(kk < 0 as libc::c_int) {
                break;
            }
        }
        if kk != j {
            loop {
                k = kk;
                kk = *np.offset(k as isize);
                *np.offset(k as isize) = -kk;
                if !(kk != j) {
                    break;
                }
            }
            k3 = kk
        } else {
            *np.offset(j as isize) = -j;
            if !(j != nn) {
                break;
            }
        }
    }
    maxf_0 *= inc;
    loop {
        j = k3 + 1 as libc::c_int;
        nt = nt - kspnn;
        i = nt - inc + 1 as libc::c_int;
        if !(nt >= 0 as libc::c_int) {
            break;
        }
        loop
        /* reorder a and b, following the permutation cycles */
        {
            loop {
                j -= 1;
                if !(*np.offset(j as isize) < 0 as libc::c_int) {
                    break;
                }
            }
            jj = jc;
            loop {
                kspan = imin2(jj, maxf_0);
                jj -= kspan;
                k = *np.offset(j as isize);
                kk = jc * k + i + jj;
                k1 = kk + kspan;
                k2 = 1 as libc::c_int;
                while k1 != kk {
                    *at.offset(k2 as isize) = *a.offset(k1 as isize);
                    *bt.offset(k2 as isize) = *b.offset(k1 as isize);
                    k1 -= inc;
                    k2 += 1
                }
                loop {
                    k1 = kk + kspan;
                    k2 = k1 - jc * (k + *np.offset(k as isize));
                    k = -*np.offset(k as isize);
                    loop {
                        *a.offset(k1 as isize) = *a.offset(k2 as isize);
                        *b.offset(k1 as isize) = *b.offset(k2 as isize);
                        k1 -= inc;
                        k2 -= inc;
                        if !(k1 != kk) {
                            break;
                        }
                    }
                    kk = k2;
                    if !(k != j) {
                        break;
                    }
                }
                k1 = kk + kspan;
                k2 = 1 as libc::c_int;
                while k1 > kk {
                    *a.offset(k1 as isize) = *at.offset(k2 as isize);
                    *b.offset(k1 as isize) = *bt.offset(k2 as isize);
                    k1 -= inc;
                    k2 += 1
                }
                if !(jj != 0 as libc::c_int) {
                    break;
                }
            }
            if !(j != 1 as libc::c_int) {
                break;
            }
        }
    }
}
/* fftmx */
static mut old_n: libc::c_int = 0 as libc::c_int;
static mut nfac: [libc::c_int; 20] = [0; 20];
static mut m_fac: libc::c_int = 0;
static mut kt: libc::c_int = 0;
static mut maxf: libc::c_int = 0;
static mut maxp: libc::c_int = 0;
/* At the end of factorization,
 * nfac[] contains the factors,
 * m_fac contains the number of factors and
 * kt contains the number of square factors  */
/* non-API, but used by package RandomFields */
#[no_mangle]
pub unsafe extern "C" fn fft_factor(
    mut n: libc::c_int,
    mut pmaxf: *mut libc::c_int,
    mut pmaxp: *mut libc::c_int,
) {
    /* fft_factor - factorization check and determination of memory
     *  requirements for the fft.
     *
     * On return, *pmaxf will give the maximum factor size
     * and  *pmaxp will give the amount of integer scratch storage required.
     *
     * If *pmaxf == 0, there was an error, the error type is indicated by *pmaxp:
     *
     *  If *pmaxp == 0  There was an illegal zero parameter among nseg, n, and nspn.
     *  If *pmaxp == 1  There we more than 15 factors to ntot.  */
    let mut j: libc::c_int = 0;
    let mut jj: libc::c_int = 0;
    let mut k: libc::c_int = 0;
    let mut sqrtk: libc::c_int = 0;
    let mut kchanged: libc::c_int = 0;
    /* check series length */
    if n <= 0 as libc::c_int {
        old_n = 0 as libc::c_int;
        *pmaxf = 0 as libc::c_int;
        *pmaxp = 0 as libc::c_int;
        return;
    } else {
        old_n = n
    }
    /* determine the factors of n */
    m_fac = 0 as libc::c_int; /* k := remaining unfactored factor of n */
    k = n;
    if k == 1 as libc::c_int {
        return;
    }
    /* extract square factors first ------------------ */
    /* extract 4^2 = 16 separately
     * ==> at most one remaining factor 2^2 = 4, done below */
    while k % 16 as libc::c_int == 0 as libc::c_int {
        let fresh0 = m_fac;
        m_fac = m_fac + 1;
        nfac[fresh0 as usize] = 4 as libc::c_int;
        k /= 16 as libc::c_int
    }
    /* extract 3^2, 5^2, ... */
    kchanged = 0 as libc::c_int;
    sqrtk = sqrt(k as libc::c_double) as libc::c_int;
    j = 3 as libc::c_int;
    while j <= sqrtk {
        jj = j * j;
        while k % jj == 0 as libc::c_int {
            let fresh1 = m_fac;
            m_fac = m_fac + 1;
            nfac[fresh1 as usize] = j;
            k /= jj;
            kchanged = 1 as libc::c_int
        }
        if kchanged != 0 {
            kchanged = 0 as libc::c_int;
            sqrtk = sqrt(k as libc::c_double) as libc::c_int
        }
        j += 2 as libc::c_int
    }
    if k <= 4 as libc::c_int {
        kt = m_fac;
        nfac[m_fac as usize] = k;
        if k != 1 as libc::c_int {
            m_fac += 1
        }
    } else {
        if k % 4 as libc::c_int == 0 as libc::c_int {
            let fresh2 = m_fac;
            m_fac = m_fac + 1;
            nfac[fresh2 as usize] = 2 as libc::c_int;
            k /= 4 as libc::c_int
        }
        /* all square factors out now, but k >= 5 still */
        kt = m_fac;
        maxp = imax2(kt + kt + 2 as libc::c_int, k - 1 as libc::c_int);
        j = 2 as libc::c_int;
        loop {
            if k % j == 0 as libc::c_int {
                let fresh3 = m_fac;
                m_fac = m_fac + 1;
                nfac[fresh3 as usize] = j;
                k /= j
            }
            if j > 2147483647 as libc::c_int - 2 as libc::c_int {
                break;
            }
            j = (j + 1 as libc::c_int) / 2 as libc::c_int * 2 as libc::c_int + 1 as libc::c_int;
            if !(j <= k) {
                break;
            }
        }
    }
    if m_fac <= kt + 1 as libc::c_int {
        maxp = m_fac + kt + 1 as libc::c_int
    }
    if m_fac + kt > 20 as libc::c_int {
        /* error - too many factors */
        old_n = 0 as libc::c_int;
        *pmaxf = 0 as libc::c_int;
        *pmaxp = 0 as libc::c_int;
        return;
    } else {
        if kt != 0 as libc::c_int {
            j = kt;
            while j != 0 as libc::c_int {
                j -= 1;
                let fresh4 = m_fac;
                m_fac = m_fac + 1;
                nfac[fresh4 as usize] = nfac[j as usize]
            }
        }
        maxf = nfac[(m_fac - kt - 1 as libc::c_int) as usize];
        /* The last squared factor is not necessarily the largest PR#1429 */
        if kt > 0 as libc::c_int {
            maxf = imax2(nfac[(kt - 1 as libc::c_int) as usize], maxf)
        }
        if kt > 1 as libc::c_int {
            maxf = imax2(nfac[(kt - 2 as libc::c_int) as usize], maxf)
        }
        if kt > 2 as libc::c_int {
            maxf = imax2(nfac[(kt - 3 as libc::c_int) as usize], maxf)
        }
    }
    *pmaxf = maxf;
    *pmaxp = maxp;
}
#[no_mangle]
pub unsafe extern "C" fn fft_work(
    mut a: *mut libc::c_double,
    mut b: *mut libc::c_double,
    mut nseg: libc::c_int,
    mut n: libc::c_int,
    mut nspn: libc::c_int,
    mut isn: libc::c_int,
    mut work: *mut libc::c_double,
    mut iwork: *mut libc::c_int,
) -> Rboolean {
    let mut nf: libc::c_int = 0;
    let mut nspan: libc::c_int = 0;
    let mut ntot: libc::c_int = 0;
    /* check that factorization was successful */
    if old_n == 0 as libc::c_int {
        return FALSE;
    }
    /* check that the parameters match those of the factorization call */
    if n != old_n || nseg <= 0 as libc::c_int || nspn <= 0 as libc::c_int || isn == 0 as libc::c_int
    {
        return FALSE;
    }
    /* perform the transform */
    nf = n;
    nspan = nf * nspn;
    ntot = nspan * nseg;
    fftmx(
        a,
        b,
        ntot,
        nf,
        nspan,
        isn,
        m_fac,
        kt,
        &mut *work.offset(0 as libc::c_int as isize),
        &mut *work.offset(maxf as isize),
        &mut *work
            .offset((2 as libc::c_int as libc::c_ulong).wrapping_mul(maxf as size_t) as isize),
        &mut *work
            .offset((3 as libc::c_int as libc::c_ulong).wrapping_mul(maxf as size_t) as isize),
        iwork,
        nfac.as_mut_ptr(),
    );
    return TRUE;
}