double simpson(unsigned int n,
        double a,double b,
        double (*foo)(double)) {
    unsigned int k;
    double x,s_odd=0,s_even=0,h=(b-a)/(n+n);
    for (k=1,x=a+h+h; k<n; ++k,x=x+h+h)
        s_even+=foo(x);
    for (k=0,x=a+h; k<n; ++k,x=x+h+h)
        s_odd+=foo(x);
    return h*(foo(a)+foo(b)+2*s_even+4*s_odd)/3;
}