normal distr. simpsons

[nikki1337]

I am looking to modify this piece of code to use simpson's rule

  double integral(double a, double b, double(*f)(double))
{
	double sum = 0;
	double err = 1e10; 
	double c = (b-a)/n;

	for (double k = 1.; k < err-1; k+=1.)
		sum += f(a + k*c);

	return c * ((f(a) + f(b)) / 2 + sum);
}
const double pi = 3.14159265359;

class NormalDistrib
{
public:
    NormalDistrib(double _mu, double _sigma) : mu(_mu), sigma(_sigma) {}
    inline double pdf(double x);
    double cdf(double x);
    double mu;
    double sigma;
};

inline double NormalDistrib::pdf(double x)
{ 
    return exp( -1 * (x - mu) * (x - mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * pi));
}

double NormalDistribution::cdf(double x)
{
    const double y = mu - 10 * sigma;
    double err = 1e10;
    double sum = 0;
    double c = (x - y) / err;
    for (double k = 1; k < err-1; k++)
        {
        sum += pdf( y + k*c);
        }
    return c * ((pdf(x) + pdf(y))/2 + sum);
}