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

#define VECTOR(first,last) (((double *) malloc(sizeof(double)*((last)-(first)+1)))-(first))

static void poly1(obs, deg, x, z, coef)
int obs, deg;
double *x, *z, *coef;
{
  int e, k, t;
  double alfa, beta, D, OLDSQSUM, SSQPOL, XPROD, ZPROD, olda;
  double *error, *ORPOL2, *ORPOL1, *POLC1, *POLC2;

  error=VECTOR(1,obs);
  ORPOL2=VECTOR(1,obs);
  ORPOL1=VECTOR(1,obs);
  POLC1=VECTOR(-1,deg);
  POLC2=VECTOR(-1,deg);

  for(k=1; k<=obs; k++)
  {
    ORPOL2[k] = 1;
    ORPOL1[k] = 0;
    error[k]  = z[k];
  }

  POLC1[-1] = 0;
  olda      = 0;
  alfa      = 0;
  OLDSQSUM = 1;
  beta = OLDSQSUM;

  for(t=0; t<=deg; t++)
  {
    XPROD  = 0;
    ZPROD  = 0;
    SSQPOL = 0;
#pragma ivdep
    for(k=1; k<=obs; k++)
    {
      error[k] = error[k] - olda*ORPOL1[k];
      D = ORPOL2[k]*beta + ORPOL1[k]*(x[k]+alfa);
      ORPOL2[k] = ORPOL1[k];
      ORPOL1[k] = D;
      ZPROD = ZPROD+error[k]*D;
      D = D*D;
      SSQPOL = SSQPOL+D;
      XPROD = XPROD+D*x[k];
    }
    coef[t] = ZPROD/SSQPOL;
    olda = coef[t];
    POLC1[t] = 1;
    POLC2[t] = 0;
    for(e=(t-1); e>=0; e--)
    {
      D        = POLC2[e]*beta;
      POLC2[e] = POLC1[e];
      POLC1[e] = POLC1[e]*alfa + D + POLC1[e-1];
      coef[e]  = coef[e] + coef[t]*POLC1[e];
    }
    beta = -SSQPOL/OLDSQSUM;
    OLDSQSUM = SSQPOL;
    alfa = -XPROD/OLDSQSUM;
  }
  free((char *) (POLC2-1));
  free((char *) (POLC1-1));
  free((char *) (ORPOL1+1));
  free((char *) (ORPOL2+1));
  free((char *) (error+1));
}

#define MAXN 10000
#define MAXD 3

void main(int argc, char **argv)
{
  double radius,finalRadius,thickness,length,finalLength;
  int rotation,n,deg,i,j;

  double *RADIUS,*LENGTH,*ROTATION,*COEF,maxerror,X,Y,Ycalc,ERROR;
  n=0;
  RADIUS=VECTOR(1,MAXN);
  LENGTH=VECTOR(1,MAXN);
  ROTATION=VECTOR(1,MAXN);
  COEF=VECTOR(0,MAXD);

  radius=11.3/2.0; /* cm */
  radius=5.0;
  radius=10.0;
  finalRadius=26.5/2.0;
  thickness=1.5e-3*2.54;
  thickness=1.0e-3;
  length=0.0;
  finalLength=2400.0*12.0*2.54;
  rotation=0;
  n++;
  RADIUS[n]=radius+rotation*thickness;
  LENGTH[n]=length;
  ROTATION[n]=rotation;
  while(length<finalLength)
  {
    rotation++;
    length+=2.0*M_PI*(radius+rotation*thickness);
//  printf("%d\t%.18lf\t%.18lf\n", rotation,radius+rotation*thickness,length);
    if(rotation%50==0)
    {
      n++;
      RADIUS[n]=radius+rotation*thickness;
      LENGTH[n]=length;
      ROTATION[n]=rotation;
    }
  }
  printf("length=F(radius):\n");
  for(deg=2;deg<=MAXD;deg++)
  {
    poly1(n, deg, RADIUS, LENGTH, COEF);
    for(i=0;i<=deg;i++) printf("%.10lg\t",COEF[i]);
    printf("\n");
    maxerror=0.0;
    for(i=1;i<=n;i++)
    {
      X=RADIUS[i];
      Y=LENGTH[i];
      Ycalc=0.0;
      for(j=deg;j>=0;j--) Ycalc=Ycalc*X+COEF[j];
      ERROR=fabs(Y-Ycalc);
      if(ERROR>maxerror) maxerror=ERROR;
    }
    printf("%d\t%.10lg\n",deg,maxerror);
  }
  printf("length=F(rotation):\n");
  for(deg=2;deg<=MAXD;deg++)
  {
    poly1(n, deg, ROTATION, LENGTH, COEF);
    for(i=0;i<=deg;i++) printf("%.10lg\t",COEF[i]);
    printf("\n");
    maxerror=0.0;
    for(i=1;i<=n;i++)
    {
      X=ROTATION[i];
      Y=LENGTH[i];
      Ycalc=0.0;
      for(j=deg;j>=0;j--) Ycalc=Ycalc*X+COEF[j];
      ERROR=fabs(Y-Ycalc);
      if(ERROR>maxerror) maxerror=ERROR;
    }
    printf("%d\t%.10lg\n",deg,maxerror);
  }
}