/*
 *  R : A Computer Language for Statistical Data Analysis

 *  Copyright (C) 1999-2016  The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/.
 */

#include <R.h>

static void
burg(int n, double*x, int pmax, double *coefs, double *var1, double *var2)
{
    double d, phii, *u, *v, *u0, sum;

    u = (double *) R_alloc(n, sizeof(double));
    v = (double *) R_alloc(n, sizeof(double));
    u0 = (double *) R_alloc(n, sizeof(double));

    for(int i = 0; i < pmax*pmax; i++) coefs[i] = 0.0;
    sum = 0.0;
    for(int t = 0; t < n; t++) {
	u[t] = v[t] = x[n - 1 - t];
	sum += x[t] * x[t];
    }
    var1[0] = var2[0] = sum/n;
    for(int p = 1; p <= pmax; p++) { /* do AR(p) */
	sum = 0.0;
	d = 0;
	for(int t = p; t < n; t++) {
	    sum += v[t]*u[t-1];
	    d += v[t]*v[t] + u[t-1]*u[t-1];
	}
	phii = 2*sum/d;
	coefs[pmax*(p-1) + (p-1)] = phii;
	if(p > 1)
	    for(int j = 1; j < p; j++)
		coefs[p-1 + pmax*(j-1)] =
		    coefs[p-2 + pmax*(j-1)] - phii* coefs[p-2 + pmax*(p-j-1)];
	/* update u and v */
	for(int t = 0; t < n; t++)
	    u0[t] = u[t];
	for(int t = p; t < n; t++) {
	    u[t] = u0[t-1] - phii * v[t];
	    v[t] = v[t] - phii * u0[t-1];
	}
	var1[p] = var1[p-1] * (1 - phii * phii);
	d = 0.0;
	for(int t = p; t < n; t++) d += v[t]*v[t] + u[t]*u[t];
	var2[p] = d/(2.0*(n-p));
    }
}

#include <Rinternals.h>

SEXP Burg(SEXP x, SEXP order)
{
    x = PROTECT(coerceVector(x, REALSXP));
    int n = LENGTH(x), pmax = asInteger(order);
    SEXP coefs = PROTECT(allocVector(REALSXP, pmax * pmax)),
	var1 = PROTECT(allocVector(REALSXP, pmax + 1)),
	var2 = PROTECT(allocVector(REALSXP, pmax + 1));
    burg(n, REAL(x), pmax, REAL(coefs), REAL(var1), REAL(var2));
    SEXP ans = PROTECT(allocVector(VECSXP, 3));
    SET_VECTOR_ELT(ans, 0, coefs);
    SET_VECTOR_ELT(ans, 1, var1);
    SET_VECTOR_ELT(ans, 2, var2);
    UNPROTECT(5);
    return ans;
}