#include #include /* gauss1d.c : 1999 3 27 by T.Oguni */ double f(double x) { return 1./(x*x*x); } double gauss1d(a, b) double a, b; { double r[2]={.652145158, .3478548451}, s[2]={.3399810435, .8611363115}; double x2, x1, y, u=0.0, q; int k; for (k=0; k<2; k++){ x1 =(s[k]*(b-a)+a+b)/2.; x2 = (-s[k]*(b-a)+a+b)/2.; y=f(x1)+f(x2); u=u+r[k]*y; } q=(b-a)*u/2.; return q; } /* main */ main() { printf("f(x)dx = %f\n", gauss1d(1., 2.)); } /* solution = 0.37 */