#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void gauss(double a[12][12],double b[12],double x[12]);
void main()
{
	double x[13]={0.00,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00};
	double fx[13]={0.00,5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00};
	double fx1[12],fx2[11],fx3[10],px1[191],px2[191],px3[191];
	double a;
	int i,k,j,m;
	for(i=0;i<12;i++)    //一階差商
		fx1[i]=(fx[i+1]-fx[i])/(x[i+1]-x[i]);
    for(i=0;i<11;i++)    //二階差商
		fx2[i]=(fx1[i+1]-fx1[i])/(x[i+2]-x[i]);
	for(i=0;i<10;i++)    //三階差商
		fx3[i]=(fx2[i+1]-fx2[i])/(x[i+3]-x[i]);
    
//分段線性插值
	for(k=0;k<=190;k++)
	{
		for(j=0;j<12;j++)
		{if(k>=x[j]&&k<x[j+1])
			m=j;
		}
        px1[k]=fx[m]+fx1[m]*(k-x[m]);
        printf("%d  %f \n",k,px1[k]);
	}
//分段二次多項式插值	 	
    for(k=0;k<=190;k++)
	{   if(k<(x[1]+x[2])/2)
	    px2[k]=fx[0]+fx1[0]*(k-x[0])+fx2[0]*(k-x[0])*(k-x[1]);
		else 
		{   if(k>=(x[10]+x[11])/2)
		    px2[k]=fx[10]+fx1[10]*(k-x[10])+fx2[10]*(k-x[10])*(k-x[11]);	
	        else
			{for(j=1;j<10;j++)
			{if(k>=(x[j]+x[j+1])/2&&k<(x[j+1]+x[j+2])/2)
			  m=j;
			}
            px2[k]=fx[m]+fx1[m]*(k-x[m])+fx2[m]*(k-x[m])*(k-x[m+1]);
			}
		}
        printf("%d  %f \n",k,px2[k]);
	}
  ////分段三次多項式插值
   for(k=0;k<=190;k++)
	{   if(k<x[2])
	    px3[k]=fx[0]+fx1[0]*(k-x[0])+fx2[0]*(k-x[0])*(k-x[1])+fx3[0]*(k-x[0])*(k-x[1])*(k-x[2]);
		else
		{  if(k>=x[10])
		   px3[k]=fx[9]+fx1[9]*(k-x[9])+fx2[9]*(k-x[9])*(k-x[10])+fx3[9]*(k-x[9])*(k-x[10])*(k-x[11]);	
	       else
		   {for(j=1;j<9;j++)
		   {if(k>=x[j+1]&&k<x[j+2])
			 m=j;
		   }
           px3[k]=fx[m]+fx1[m]*(k-x[m])+fx2[m]*(k-x[m])*(k-x[m+1])+fx3[m]*(k-x[m])*(k-x[m+1])*(k-x[m+2]);
		   }
		}
        printf("%d  %f \n",k,px3[k]);
	}
//三次樣條插值
   double h[12],alpha[12],beta[12],gama[12],M[13],N[12],A[12][12]={0},pxn[191];
	double f;
	for(i=0;i<12;i++)
		h[i]=x[i+1]-x[i];
	for(i=0;i<11;i++)
	{
		alpha[i]=h[i+1]/(h[i]+h[i+1]);
		gama[i]=1-alpha[i];
	}
	
	for(i=1;i<12;i++)
		beta[i-1]=6*((fx[i+1]-fx[i])/h[i]-(fx[i]-fx[i-1])/h[i-1])/(h[i]+h[i-1]);

	alpha[11]=h[0]/(h[0]+h[11]);
	gama[11]=1-alpha[11];
	beta[11]=6*((fx[1]-fx[0])/h[0]-(fx[12]-fx[11])/h[11])/(h[0]+h[11]);

	for(i=0;i<11;i++)
	{
		A[i][i]=2;
		A[i+1][i]=gama[i+1];
		A[i][i+1]=alpha[i];
	}
	A[11][11]=2;
	A[0][11]=gama[0];
	A[11][0]=alpha[11];

	gauss(A,beta,N);

	for(i=1;i<13;i++)
		M[i]=N[i-1];
	M[0]=M[12];

	
for(k=0;k<=190;k++)
{for(i=1;i<=12;i++)
		{if(k>=x[i-1]&&k<x[i])
			j=i;
		}
 pxn[k]=M[j-1]*(x[j]-k)*(x[j]-k)*(x[j]-k)/(6*h[j-1])+M[j]*(k-x[j-1])*(k-x[j-1])*(k-x[j-1])/(6*h[j-1])+(fx[j-1]/h[j-1]-M[j-1]*h[j-1]/6)*(x[j]-k)+(fx[j]/h[j-1]-M[j]*h[j-1]/6)*(k-x[j-1]);
  printf("%d  %f \n",k,pxn[k]);
} 
}
//gauss消元
#define e 0.000000001

void gauss(double a[12][12],double b[12],double x[12])
{
	int k,i,j,n=12;
	double sum;
	double m[12];

	/* 消元過程 */
	for(k=0;k<n-1;k++)
	{
		if(a[k][k]==0)
			break;
		else
		{
			for(i=k+1;i<n;i++)
			{
				m[i]=a[i][k]/a[k][k];
				for(j=k;j<n;j++)
					a[i][j]=a[i][j]-m[i]*a[k][j];
				b[i]=b[i]-m[i]*b[k];
			}
		}
	}
	/* 回帶過程 */
	x[n-1]=b[n-1]/a[n-1][n-1];
	for(i=n-2;i>=0;i--)
	{   
		for(sum=0.0,j=i+1;j<n;j++) 
			sum=sum+a[i][j]*x[j];
		x[i]=b[i]-sum;
		x[i]=x[i]/a[i][i];
	}
}