程序代码：

#include <functional>

double SimpsonIntegration1d( const std::function<double(double)>& f, double xa, double xb, unsigned n )
{
    double sum = f(xa) + f(xb);
    for( unsigned i=1; i<2*n; ++i )
        sum += (i%2*2+2) * f((xb-xa)/(2*n)*i+xa);
    return (xb-xa)/(6*n) * sum;
}

double SimpsonIntegration2d( const std::function<double(double,double)>& f, double xa, double xb, double ya, double yb, unsigned n )
{
    auto g = [=](double y) {
        auto h = [=](double x){ return f(x,y); };
        return SimpsonIntegration1d(h,xa,xb,n);
    };
    return SimpsonIntegration1d(g,ya,yb,n);
}

#include <stdio.h>

int main( void )
{
    auto f = [](double x, double y){ return x*x+y*y; };
    printf( "%f\n", SimpsonIntegration2d(f,-1,+1,-1,+1,1000) );
}

程序代码：

#include <stdio.h>

double SimpsonIntegration1d( double (*f)(double x), double xa, double xb, unsigned n )
{
    double sum = f(xa) + f(xb);
    for( unsigned i=1; i<2*n; ++i )
        sum += (i%2*2+2) * f((xb-xa)/(2*n)*i+xa);
    return (xb-xa)/(6*n) * sum;
}

double SimpsonIntegration2d_( double (*f)(double x,double y), double y, double xa, double xb, unsigned n )
{
    double sum = f(xa,y) + f(xb,y);
    for( unsigned i=1; i<2*n; ++i )
        sum += (i%2*2+2) * f((xb-xa)/(2*n)*i+xa,y);
    return (xb-xa)/(6*n) * sum;
}

double SimpsonIntegration2d( double (*f)(double x,double y), double xa, double xb, double ya, double yb, unsigned n )
{
    double sum = SimpsonIntegration2d_(f,ya,xa,xb,n) + SimpsonIntegration2d_(f,yb,xa,xb,n);
    for( unsigned i=1; i<2*n; ++i )
        sum += (i%2*2+2) * SimpsonIntegration2d_(f,(yb-ya)/(2*n)*i+ya,xa,xb,n);
    return (yb-ya)/(6*n) * sum;
}

#include <stdio.h>

double f( double x, double y )
{
    return x*x + y*y;
}

int main( void )
{
    printf( "%f\n", SimpsonIntegration2d(f,-1,+1,-1,+1,1000) );
}