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

long * seq;
long * apr;
int   size;

double * coeff;
int    degree;

long ** sdtb,** adtb;

double poly_calc(long ind)
{
  long x;
  int i;
  double res=0;

  for(x=1,i=0;i<degree;i++,x*=ind)
    res+=x*coeff[i];
  
  return res;
}

long der(long * dt,int m,int p)
{
  long ** dtbl=(dt==seq)?sdtb:adtb;

  if(m>=degree)
    return 0;

  if((p<=0)||(p>=size-m))
    return 0;

  if(m==0)
    return dt[p];

  return dtbl[m-1][p];
}

double der_avg(long * dt,int mag)
{
  long long dsum=0;
  int i;
  
  if(mag==size)
    return 0;

  for(i=mag;i<size;i++)
    dsum+=der(dt,mag,i);
  
  return dsum/(size-mag);
}

double der_diff(int mag)
{
  long long dsum=0;
  int i;
  
  if(mag==size)
    return 0;

  for(i=mag;i<size;i++)
    dsum+=der(seq,mag,i)-der(apr,mag,i);
  
  return dsum/(size-mag);
  
}

void approximate()
{
  int i;

  apr[0]=seq[0];

  for(i=1;i<size;i++)
    apr[i]=rint(poly_calc(degree));
}

void derive()
{
  int i,j;

  for(i=1;i<size-1;i++)
    adtb[0][i]=apr[i]-apr[i-1];
  for(j=1;j<degree;j++)
    for(i=1;i<size-1-j;i++)
      adtb[j][i]=der(apr,j-1,i);
}

int main(int argc,char ** argv)
{
  int i,j;

  if(argc!=3)
    {
      fprintf(stderr,"Usage: %s <number of samples> <polynomial degree>\n",argv[0]);
      return -1;
    }
  size=atoi(argv[1]);
  degree=atoi(argv[2]);

  puts("Allocating memory ...");

  if(!(size&&degree))
    {
      fprintf(stderr,"Usage: %s <number of samples> <polynomial degree>\n",argv[0]);
      return -1;
    }

  if((seq=malloc(size*sizeof(long)))==NULL)
    {
      fprintf(stderr,"Error allocating %d bytes for samples array\n",size*sizeof(long));
      return -1;
    }
  if((apr=malloc(size*sizeof(long)))==NULL)
    {
      fprintf(stderr,"Error allocating %d bytes for approximation array\n",size*sizeof(long));
      return -1;
    }
  if((coeff=malloc(degree*sizeof(double)))==NULL)
    {
      fprintf(stderr,"Error allocating %d bytes for polynomial coefficients\n",degree*sizeof(double));
      return -1;
    }

  if((sdtb=malloc(degree*sizeof(long *)))==NULL)
    {
      fprintf(stderr,"Error allocating %d bytes for data derivative table index\n",degree*sizeof(long*));
      return -1;
    }
  if((adtb=malloc(degree*sizeof(long *)))==NULL)
    {
      fprintf(stderr,"Error allocating %d bytes for approximation derivative table index\n",degree*sizeof(long*));
      return -1;
    }


  for(i=0;i<degree;i++)
    if((sdtb[i]=malloc((degree-i-1)*sizeof(long)))==NULL)
      {
	fprintf(stderr,"Error allocating %d bytes for data derivative table\n",(degree-i-1)*sizeof(long));
	return -1;
      }

  for(i=0;i<degree;i++)
    if((adtb[i]=malloc((degree-i-1)*sizeof(long)))==NULL)
      {
	fprintf(stderr,"Error allocating %d bytes for approximation derivative table\n",(degree-i-1)*sizeof(long));
	return -1;
      }
    

  puts("Initializing data structures ...");
  memset(apr,0,sizeof(long)*size);
  memset(coeff,0,sizeof(double)*degree);

  puts("Reading data ...");

  for(i=0;i<size;i++)
    if(!scanf("%ld",&seq[i]))
    {
      fprintf(stderr,"Error reading data element %d !!!\n",i);
      return -2;
    }
  printf("Got %d data elements\n",i);

  puts("Creating data derivative table ...");
  for(i=1;i<size-1;i++)
    sdtb[0][i]=seq[i]-seq[i-1];
  for(j=1;j<degree;j++)
    for(i=1;i<size-1-j;i++)
      sdtb[j][i]=der(seq,j-1,i);



  puts("Setting initial polynomial coefficients ...");
  for(i=0;i<degree;i++)
    coeff[i]=der_avg(seq,i);
  puts("Entering homing approximation loop ...");
  printf("Approximation: [");

  while(0!=der_diff(0))
    {
      approximate();
      derive();
      for(i=0;i<degree;i++)
	coeff[i]-=der_diff(i);
      putc('*',stdout);
    }
  printf("] %d\n",i);
  for(i=0;i<degree;i++)
    printf("K[%d] = %f\n",i,coeff[i]);
  return 0;
}