解题的过程和复化SIMPSON公式

#include "stdio.h"
#include "math.h"

double P4(double x)
{
\x09return (1.0 + x + x*x/2.0 + pow(x,3)/6.0 + pow(x,4)/24.0);
}

double G(double x)
{
\x09if ( 0<x && x<=1)
\x09{
\x09\x09return ( (exp(x)-P4(x))/sqrt(x) );
\x09}else
\x09{
\x09\x09return 0;
\x09}
}

double H()
{
\x09int i;
\x09double f1=0,f2=0;
\x09double m[]={     1,    1,1.0/2,1.0/6,1.0/24};
\x09double n[]={-1.0/2,1.0/2,3.0/2,5.0/2,7.0/2 };
\x09for(i=0;i<5;i++)
\x09{
\x09\x09m[i]=m[i]/(n[i]+1);
\x09\x09n[i]=n[i]+1;
\x09}
\x09for(i=0;i<5;i++)
\x09{
\x09\x09f1+=m[i]*pow(0,n[i]);
\x09\x09f2+=m[i]*pow(1,n[i]);
\x09}
\x09
\x09return (f2-f1);
}

double Simpson(double a, double b, double h, double (*f)(double))
{
\x09double xk;
\x09double sk = 0,sk2 = 0;
\x09double Sn = 0; 
\x09double n;
\x09int i;
\x09n = (b-a)/h;
\x09sk2 = f(a+h/2.0);
\x09for (i=1; i<n; i++)
\x09{
\x09\x09xk = a+h*i;
\x09\x09sk += f(xk);
\x09\x09sk2 += f(xk+h/2.0);
\x09}
\x09Sn = f(a) + 2.0*sk + f(b) + 4.0*sk2;
\x09Sn = Sn*h/6.0; 
\x09return Sn;
}

double Int()
{
\x09return (H()+Simpson(0, 1, 0.25, G));
}

void main()
{
\x09printf("H = %lf\n", H());
\x09printf("Int(Gx) = %lf\n", Simpson(0, 1, 0.25, G));
\x09printf("answer = %lf\n", Int());
}