/*
 ============================================================================
 Name        : korrelacio.c
 Author      : Stippinger Marcell
 Version     : 1.0
 Copyright   : BME-TTK FI SzamSzim gyak segedanyag
 Description : Autocorrelation and momentum analysis in C, Ansi-style
 ============================================================================
 */


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

const int 	k=1000,		/* 1+the shift of the farthest data correlation calculated */
			m=20;		/* The index of the highest momentum calculated            */

/* We consider the original and the time-shifted data series as two
 * genuinely independent data series. That is we truncate the data
 * to the overlapping region:
 *    1   2   3   4   ...  N-1  N			(first)
 *            1   2   ...  N-3 N-2 N-1  N	(second)
 *           \_________ _________/
 *                     V
 *           These are the two independent data series we evaluate
 *
 * Therefore the average and the variance of the two truncated
 * series will differ, we deal ith this in the code.
 */
double  *correlation,  /* Result of the correlation calculation   */
		*sumhead,		/* Sum of the first i input (head)        */
		*sum2head,		/* Sum of the square of the first i input */
		*sumtail,		/* Sum of the last i input (tail)         */
		*sum2tail,		/* Sum of the square of the last i input  */
		sumall=0,		/* Sum of the all the input               */
		sum2all=0,		/* Sum of the square of all the inplut    */
		*shiftproduct,	/* sum_{i=k}^{N} X(i)X(N-i) type sums     */
		*values,		/* For computing purposes we have to store
						 * the last k values, we overwrite these
						 * cyclically
						 */
		*momentum;		/* The ith momentum * N                   */

int n=0;  /* Number of input */

void myread();
void mywrite();

int main()
{
	/* Memory allocation */
	correlation	= calloc(k,sizeof(double));
	sumhead	= calloc(k+1,sizeof(double));
	sum2head	= calloc(k+1,sizeof(double));
	sumtail		= calloc(k+1,sizeof(double));
	sum2tail	= calloc(k+1,sizeof(double));
	shiftproduct= calloc(k,sizeof(double));
	values		= calloc(k,sizeof(double));
	momentum	= calloc(m+1,sizeof(double));

	/* Check success */
	if ((correlation==NULL)||(sumhead==NULL)||(sum2head==NULL)||
			(sumtail==NULL)||(sum2tail==NULL)||(shiftproduct==NULL)||
			(values==NULL)||(momentum==NULL))
	{
		printf("Memory allocation failed"); exit(1);
	}

	/* Call the data analysis */
	myread();
	/* Call the output module */
	mywrite();

	/* Memory freeing */
	correlation	= calloc(k,sizeof(double));
	sumhead	= calloc(k+1,sizeof(double));
	sum2head	= calloc(k+1,sizeof(double));
	sumtail		= calloc(k+1,sizeof(double));
	sum2tail	= calloc(k+1,sizeof(double));
	shiftproduct= calloc(k,sizeof(double));
	values		= calloc(k,sizeof(double));
	momentum	= calloc(m+1,sizeof(double));

	/* Return code */
	return 0;
}

/* Read and analyze data */
void myread()
{
	double x,	/* Input variable */
		   meanfirst, meansecond, varfirst, varsecond; /* Temporary variables */
	int i,		/* Time shift index */
		j;		/* Momentum index */

	/* Initialize cumulative sum of first input */
	sumhead[0]=0.;
	sum2head[0]=0.;

	/* Read the input as long as there is any */
	while(!feof(stdin))
	{
		/* Check whether we could read exactly one entry
		 * (some systems indicate error with negative values)
		 */
		if(scanf("%lf",&x)==1)
		{
			/* Increment input count and update sums */
			n++;
			sumall+=x;
			sum2all+=x*x;
			/* Calculate the cumulative sum for the first (head) entries */
			if (n<=k) {
				sumhead[n]=sumhead[n-1]+x;
				sum2head[n]=sum2head[n-1]+x*x;
			}
			/* Cyclically store values, this is needed for the tail sums */
			values[n%k]=x;
			/* Update momenta */
			for (j=0; j<=m; j++) {
				momentum[j]+=pow(x,j);
			}
			/* Update correlation data */
			for (i=0; i<k; i++) {
				shiftproduct[i]+=values[n%k]*values[(n+k-i)%k];
			}
		}
	}

	/* Initialize cumulative sum of last input */
	sumtail[0]=0.;
	sum2tail[0]=0.;

	/* Calculate the cumulative sum for the last (tail) entries */
	for (i=1; (i<=k) && (i<=n); i++) {
		sumtail[i]=sumtail[i-1]+values[(n+1-i)%k];
		sum2tail[i]=sum2tail[i-1]+values[(n+1-i)%k]*values[(n+1-i)%k];
	}

	for (i=0; (i<k) && (i<n); i++) {
		/* Calculating the correlation by definition
		 * Please note that D2(X)=E((X-E(X))^2)=E(X^2-2*X*E(X)+E(X)^2)=E(X^2)-E(X)^2
		 */
		meanfirst = (sumall-sumtail[i])/(double)(n-i);	/* Mean of the first data series for shift i */
		meansecond = (sumall-sumhead[i])/(double)(n-i);	/* Mean of the second data series for shift i */
		varfirst = (sum2all-sum2tail[i])/(double)(n-i)
				- meanfirst*meanfirst;	/* Variance of the first data series for shift i */
		varsecond = (sum2all-sum2head[i])/(double)(n-i)
				- meansecond*meansecond;/* variance of second data series for shift i */

		correlation[i] =
				(shiftproduct[i]/(double)(n-i) - meanfirst*meansecond) /
					(sqrt(varfirst*varsecond));		/* Correlation by definition */
	}
}

/* Writeout */
void mywrite()
{
	int i,		/* Time shift index */
		j;		/* Momentum index */
	/* Print normalized momenta */
	printf("%s","# Momenta j=0..m\n");
	for(j=0; (j<=m); j++) {
		printf("%d\t%+8.5lf\n",j,momentum[j]/(double)n);
	}
	/* Print correlation */
	printf("%s","\n\n# Autocorrelation i=0..k-1\n");
	for(i=0; (i<k) && (i<n); i++) {
		printf("%d\t%+8.5lf\n",i,correlation[i]);
	}
}