積分計算

なんとなく必要になりそうなのでなんとなく移植しました。
http://d.hatena.ne.jp/n-trino/20091124#p1
元がhaskellなので関数型言語っぽい実装になりました。

//台形
struct Trapezoid{
	template< class FUNCTION >
	double operator()( const FUNCTION& funct, double a, double b )const{
		return (funct( a ) + funct( b )) * ( b - a ) / 2.0;
	}
};
//シンプソン
struct Simpson{
	template< class FUNCTION >
	double operator()( const FUNCTION& funct, double a, double b )const{
		return (funct( a ) + 4.0 * funct( (a + b)/2.0 ) + funct( b )) * (b - a) / 6.0;
	}
};

//数値積分
template< class FUNCTION = double(*)(double), class INTERPOLATION = Simpson >
class NumericalIntegral
{
public:
	//コンストラクタで積分したい関数および補間方法を入れる。
	NumericalIntegral( const FUNCTION& funct = FUNCTION(), const INTERPOLATION& inter = INTERPOLATION() ):
	  m_Funct(funct),m_Inter(inter){}
	virtual ~NumericalIntegral(){}
	
	//計算開始。aからbの区間を積分。2^divだけ分割する。
	template< unsigned int DIV >
	double calc( double a, double b )const{
		const double m = (a + b) / 2;
		return calc< DIV-1 >( a, m ) + calc< DIV-1 >( m, b );
	}
	template<>
	double calc< 0 >( double a, double b )const{
		return m_Inter( m_Funct, a, b );
	}
private:
	const FUNCTION m_Funct;
	const INTERPOLATION m_Inter;
};

使い方は特に説明するまでもないと思います。シンプソンがあれば台形が必要な理由はないと思いますがなんとなく作っちゃいました。
実際の計算の再帰部分ですがテンプレートでなく普通の引数にしようか迷ったのですがなんとなくテンプレートにしちゃいました。しかし実際には積分区間の長さによって調整する可能性があるので普通の引数の方が良かったかもしれません。

さて、実際に使うかどうかは不明です。