Laplace-Everett's Formula

1:  #include<stdio.h>  
2:  #include<math.h>  
3:  float factorial(int);  
4:  void main()  
5:  {  
6:       float ax[10],y[10];  
7:       float dy[10][10],u,h,x,a,sum=0,sum1=0,sum2=0,l,w,o,p,s,m;  
8:       int i,j,k,n,fact;  
9:       printf("\n\t\tLaPlace-Everett Formula");  
10:       printf("\n=======================================================\n");  
11:       printf("Enter the number of terms\n ");  
12:       scanf("%d",&n);  
13:       printf("\n=======================================================\n");  
14:       printf("Enter the array ax ");  
15:       for(i=0;i<n;i++)  
16:       {  
17:            printf("Enter the value of ax[%d] = ",i);  
18:            scanf("%f",&ax[i]);  
19:       }  
20:       printf("\n=======================================================\n");  
21:       for(i=0;i<n;i++)  
22:       {  
23:            printf("Enter the value of y[%d] = ",i);  
24:            scanf("%f",&y[i]);  
25:       }  
26:       printf("\nEnter the interpolation point ");  
27:       scanf("%f",&x);  
28:       h=ax[1]-ax[0];  
29:       for(i=0;i<n;i++)  
30:       {  
31:            for(j=0;j<=n;j++)  
32:                 dy[i][j]=-1;  
33:       }  
34:       for(i=0;i<n;i++)  
35:       {  
36:            dy[i][0]=y[i];  
37:       }  
38:       k=1;  
39:       j=1;  
40:       while(j<n)  
41:       {  
42:            i=1;  
43:            while(i<=n-k)  
44:            {  
45:                 dy[i-1][j]=dy[i][j-1]-dy[i-1][j-1];  
46:                 i++;  
47:            }  
48:            j++;  
49:            k++;  
50:       }  
51:       i=0;  
52:       printf("\n================================================================================\n");  
53:       printf("\ty");  
54:       for(i=1;i<n;i++)  
55:            printf("\td^%dy",i);  
56:       printf("\n");  
57:       for(i=0;i<n;i++)  
58:       {  
59:            for(j=0;j<n-i;j++)  
60:            {  
61:                 if(dy[i][j]==-1.0)  
62:                      printf("\t_");  
63:                 else  
64:                      printf("\t%-4.4f",dy[i][j]);  
65:            }  
66:            printf("\n");  
67:       }  
68:       i=0;  
69:       do{  
70:            i=i+1;  
71:       }while(ax[i]<x);  
72:       i--;  
73:       u=(x-ax[i])/h;  
74:       w=1-u;  
75:       sum1=sum1+u*(dy[i+1][0]);  
76:       sum2=sum2+w*(dy[i][0]);  
77:       k=2;  
78:       o=1;  
79:       s=3;  
80:       for(j=i;j>=0;j--)  
81:       {  
82:            p=u;  
83:            for(m=1;m<=o;m++)  
84:            {  
85:                 p=p*(pow(u,2)-pow(m,2));  
86:            }  
87:            o++;  
88:            p=(p/factorial(s))*dy[j][k];  
89:            sum1=sum1+p;  
90:            k=k+2;  
91:            if(k>=n)  
92:                 break;  
93:       }  
94:       k=2;  
95:       o=1;  
96:       s=3;  
97:       for(j=i-1;j>=0;j--)  
98:       {  
99:            p=w;  
100:            for(m=1;m<=o;m++)  
101:            {  
102:                 p=p*(pow(w,2)-pow(m,2));  
103:            }  
104:            p=(p/factorial(s))*dy[j][k];  
105:            sum2=sum2+p;  
106:            o++;  
107:            k=k+2;  
108:            if(k>=n)  
109:                 break;  
110:       }  
111:       sum=sum1+sum2;  
112:       printf("\nValue at %f is %f ",x,sum);  
113:       getchar();  
114:  }  
115:  float factorial(int x)  
116:  {  
117:       int i,j=1;  
118:       for(i=1;i<=x;i++)  
119:            j=j*i;       
120:       return j;  
121:  }