// quadfit.c: plot quadratic smoothness as a function of smoothing window
// Vaughan Pratt
// February 25, 2011
//
// Reads samples from stdin
// Plots 1/(1 - r2) of least squares fit of a quadratic to smoothed stdin,
// as a function of running average window width from 55 to 75 years
// Usage: quadfit < datafile > outputfile
// Datafile should contain white-spaced numbers (any number per line).
// Outputfile can be input to your favourite plotting program, e.g.
// Excel (Insert > Chart), gnuplot, etc.

#include <math.h>
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>

#define LW	(55*12)			 // Lower width for running avg
#define UW	(75*12)			 // Upper width for running avg
#define sq(x)	((x)*(x))		 // squaring macro

// Average of x[-hw], x[-hw+1], ..., x[hw-1], x[hw]
double
avg(double x[], int hw)
{
    int i;
    double ac;
    for (i = -hw, ac = 0; i <= hw; i++)
    	ac += x[i];
    return ac/(2*hw + 1);
}

// Set dst[i] = average of x[i] for i in [-n,n]
// n: half-width of x[] and dst[], hw: half-width of averaging window
void
runavg(double dst[], double x[], int n, int hw)
{
    int i;
    for (i = -n; i <= n; i++)
    	dst[i] = avg(&x[i], hw);
    return;
}

// Return 1/(1 - r2) of least-squares quadratic fit to x[-n], ..., x[n]
double
lsfitquad(double x[], int n)
{
    int i;
    double a, b, c;				// quadratic coefficients
    double xdi0, xdi1, xdi2;			// x dot [i^{0,1,2}]
    double d4;					// determinant
    double si0, si2, si4;			// sij = (-n)^j + ... + n^j
    double xsq, rsq;				// x.x, r.r
    si0 = 2*n + 1;				// (-n)^0 + ... + n^0
    si2 = ((2*n + 3)*n + 1)*n/3;		// (-n)^2 + ... + n^2
    si4 = (((6*n + 15)*n + 10.)*n*n - 1)*n/15;	// (-n)^4 + ... + n^4
    // (Don't need sij for odd j because it's zero)
    d4 = sq(si2) - si0*si4;
    for (i = -n, xdi2 = xdi1 = xdi0 = 0; i <= n; i++) {
	xdi0 += x[i];
	xdi1 += x[i]*i;
	xdi2 += x[i]*sq(i);
    }
    a = (si2*xdi0 - si0*xdi2)/d4;
    b = xdi1/si2;
    c = (si2*xdi2 - si4*xdi0)/d4;
    for (i = -n, xsq = rsq = 0; i <= n; i++) {
    	xsq += sq(x[i] - xdi0/si0);		// central x[i]
	rsq += sq(x[i] - ((a*i + b)*i + c));
    }
    return xsq/rsq;				// 1/(1 - r2)
}

int
main()
{
    int m;				// midpoint of input data
    int hw;				// half-width of smoothing window
    double rsrc[10000], rsmo[10000];	// raw arrays, room to spare
    double *src, *smo;    		// centers thereof
    for (m = 0; scanf("%lf", &rsrc[m]) == 1; m++);// read data from stdin
    m = (m - 1)/2;			// midpoint of input
    src = rsrc + m, smo = rsmo + m;	// center src, smo
    for (hw = LW/2; hw <= UW/2; hw++) {	// plot fits
	runavg(smo, src, m - hw, hw);	// smooth source
    	printf("%12f %12f\n", hw/6., lsfitquad(smo, m - hw));
    }
    return 0;    
}