//==================================================== file = moments.c =====
//=  Program to compute the first M moments for a series X of size N        =
//=   - Computes both raw (about the origin) and central moments            =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Input from input file "in.dat" to stdin (see example below)       =
//=        * Comments are bounded by "&" characters at the beginning and    =
//=          end of the comment block                                       =
//=    2) Must manually set #define M for number of moments to compute      =
//=    3) Output is to stdout                                               =
//=-------------------------------------------------------------------------=
//= Example "in.dat" file:                                                  =
//=                                                                         =
//=    & Sample series of data which can be integers or reals.              =
//=      There are 11 values in this file. &                                =
//=    50                                                                   =
//=    42                                                                   =
//=    48                                                                   =
//=    61                                                                   =
//=    60                                                                   =
//=    53                                                                   =
//=    39                                                                   =
//=    54                                                                   =
//=    42                                                                   =
//=    59                                                                   =
//=    53                                                                   =
//=-------------------------------------------------------------------------=
//= Example output (for above "in.dat" and M = 5):                          =
//=                                                                         =
//=   ----------------------------------------------- moments.c -----       =
//=     Total of 11 values (sum = 561.000000)                               =
//=       Raw moment 1 = 51.000000                                          =
//=       Raw moment  2 = 2653.545455                                       =
//=       Raw moment  3 = 140605.909091                                     =
//=       Raw moment  4 = 7572951.727273                                    =
//=       Raw moment  5 = 413780658.272727                                  =
//=   ---------------------------------------------------------------       =
//=       Central moment  1 = 0.000000                                      =
//=       Central moment  2 = 52.545455                                     =
//=       Central moment  3 = -84.545455                                    =
//=       Central moment  4 = 4973.636364                                   =
//=       Central moment  5 = -15913.636364                                 =
//=   ---------------------------------------------------------------       =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 moments.c                                                 =
//=-------------------------------------------------------------------------=
//=  Execute: moments < in.dat                                              =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (11/20/02) - Genesis (from summary1.c)                    =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf() and feof()
#include <stdlib.h>                // Needed for exit() and atof()
#include <string.h>                // Needed for strcmp()
#include <math.h>                  // Needed for pow()

//----- Defines -------------------------------------------------------------
#define MAX_SIZE 2000000           // Maximum size of time series data array
#define M              5           // Number of moments to compute

//----- Globals -------------------------------------------------------------
double     X[MAX_SIZE];            // Time series read from "in.dat"
int        N;                      // Number of values in "in.dat"

//----- Function prototypes -------------------------------------------------
void   load_X_array(void);         // Load X array

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  double   sum;                    // Sum of values
  double   mean;                   // Mean of all values
  double   raw_mom[M+1];           // Raw moments
  double   central_mom[M+1];       // Central moments
  int      i, j;                   // Loop counters

  // Load the series X
  printf("----------------------------------------------- moments.c -----\n");
  load_X_array();

  // Initialize all moments to zero
  for (i=1; i<=M; i++)
    raw_mom[i] = central_mom[i] = 0;

  // Loop to compute mean (needed for all central moment calculations) and sum
  mean = 0;
  for (i=0; i<N; i++)
  {
    sum = sum + X[i];
    mean = mean + X[i] / N;
  }

  // Loop to compute moments
  for (i=0; i<N; i++)
  {
    for (j=1; j<=M; j++)
      raw_mom[j] = raw_mom[j] + (pow(X[i], (double) j) / N);
    for (j=1; j<=M; j++)
      central_mom[j] = central_mom[j] + (pow((X[i] - mean), (double) j) / N);
  }

  // Output results
  printf("  Total of %ld values (sum = %f) \n", N, sum);
  for (i=1; i<=M; i++)
    printf("    Raw moment %2d = %f \n", i, raw_mom[i]);
  printf("---------------------------------------------------------------\n");
  for (i=1; i<=M; i++)
    printf("    Central moment %2d = %f \n", i, central_mom[i]);
  printf("---------------------------------------------------------------\n");
}

//===========================================================================
//=  Function to load X array from stdin and determine N                    =
//===========================================================================
void load_X_array(void)
{
  char      temp_string[1024];     // Temporary string variable

  // Read all values into X
  N = 0;
  while(1)
  {
    scanf("%s", temp_string);
    if (feof(stdin)) goto end;

    // This handles a comment bounded by "&" symbols
    while (strcmp(temp_string, "&") == 0)
    {
      do
      {
        scanf("%s", temp_string);
        if (feof(stdin)) goto end;
      } while (strcmp(temp_string, "&") != 0);
      scanf("%s", temp_string);
      if (feof(stdin)) goto end;
    }

    // Enter value in array and increment array index
    X[N] = atof(temp_string);
    N++;

    // Check if MAX_SIZE data values exceeded
    if (N >= MAX_SIZE)
    {
      printf("*** ERROR - greater than %ld data values \n", MAX_SIZE);
      exit(1);
    }
  }

  // End-of-file escape
  end:

  return;
}