Gauss Central Difference Formula(Forward)

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,l,olda=1;  
8:       int i,j,k,n,mid=0;  
9:       printf("\n\t\tGauss Forward Difference Formula");  
10:       printf("\n=======================================================\n");  
11:       printf("Enter the number of terms ");  
12:       scanf("%d",&n);  
13:       printf("\n=======================================================\n");  
14:       printf("Enter the array ax \n");  
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;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++;  
71:       }while(ax[i]<x);  
72:       i--;  
73:       mid=i;  
74:       u=(x-ax[mid])/h;  
75:       sum=sum+y[mid]+(u*dy[mid][1]);  
76:       a=u;  
77:       for(i=2;i<n;i++)  
78:       {  
79:            a=u;  
80:            if(i%2==1)  
81:            {  
82:                 sum=sum + ((olda/factorial(i))*dy[mid][i]);  
83:            }  
84:            else  
85:            {  
86:                 for(l=1;l<=i/2;l++)  
87:                 {  
88:                      a=a*(u-l);  
89:                 }  
90:                 if(i>=3)  
91:                      a=a*(u+1);  
92:                 mid-=1;  
93:                 sum=sum + ((a/factorial(i))*dy[mid][i]);  
94:                 olda=a;  
95:            }  
96:       }  
97:       printf("\nValue at %f is %f ",x,sum);  
98:       getchar();  
99:  }  
100:  float factorial(int x)  
101:  {  
102:       int i,j=1;  
103:       for(i=1;i<=x;i++)  
104:            j=j*i;  
105:       return j;  
106:  }