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