变步长的Simpson求积法

#include<iostream.h>
#include<fstream.h>
#include<math.h>
#include<stdlib.h>

ifstream infile("in.dat");
ofstream outfile("out.out ");

double a[3],b[3];
int i,m[3],k;

void indata()
{
for(i=0;i<3;i++)
infile>>a[i]>>b[i]>>m[i];
m[i]/=2;
}

double cal0(double x)
{
return (6*exp(-x*x));
}

double cal1(double x)
{
return (3*log(x));
}

double cal2(double x)
{
return (sin(x)/x);
}

void simpson(double a,double b,int m,double (*g)(double))
{
double h[3],q[3],p[3],s[3],t;

h[0]=(b-a)/(2*m);
t=0;
for(i=1;i<m+1;i++)
t+=(g(a+(2*i-1)*h[0])+g(a+2*i*h[0]));
p[0]=g(a)-g(b)+2*t;
t=0;
for(i=1;i<m+1;i++)
t+=g(a+(2*i-1)*h[0]);
s[0]=h[0]*(p[0]+2*t)/3;

h[1]=h[0]/2;
q[0]=0;
for(i=1;i<2*m+1;i++)
q[0]+=g(a+(2*i-1)*h[1]);
s[1]=h[1]*(p[0]+4*q[0])/3;
h[2]=h[1]/2;
q[1]=0;
for(i=1;i<pow(2,2)*m+1;i++)
q[1]+=g(a+(2*i-1)*h[2]);
p[1]=p[0]+2*q[0];
s[2]=h[2]*(p[1]+4*q[1])/3;
for(k=2;;k++)
{
h[(k+1)%3]=(b-a)/(pow(2,k+2)*m);
q[k%3]=0;
for(i=1;i<pow(2,k+1)*m+1;i++)
q[k%3]+=g(a+(2*i-1)*h[(k+1)%3]);
p[k%3]=p[(k-1)%3]+2*q[(k-1)%3];
s[(k+1)%3]=h[(k+1)%3]*(p[k%3]+4*q[k%3])/3;
if(s[(k+1)%3]-s[k%3]<1e-7)
{ t=s[(k+1)%3];
break;
}
}
outfile<<t<<endl;
outfile<<"积分区间数目为:"<<pow(2,k+1)*m<<endl;
}

void main()
{
indata();

outfile<<"I("<<1<<")=";
simpson(a[0],b[0],m[0],cal0);

outfile<<"I("<<2<<")=";
simpson(a[1],b[1],m[1],cal1);

outfile<<"I("<<3<<")=";
simpson(a[2],b[2],m[2],cal2);

infile.close();
outfile.close();
}

input==================

0 1 2
1 2 2
1 4 2

output=================

I(1)=4.48094
积分区间数目为:16
I(2)=1.15888
积分区间数目为:32
I(3)=0.81212
积分区间数目为:32