#include <stdio.h>
#include <math.h>
#define mx 513
#define my 513

double ff[mx][my]; // A large size of array is taken as an extern

double f(double x, double y)
  {
   double fx,pi;
     pi=atan(1.0)*4.0;
     fx=sin(pi*x)*sin(pi*y);
     return fx;
  }


 double simp2d(double ff[mx][my], double dx, double dy, int nx, int ny) 
  {
    double ax[mx],ay[my],s;
    int i,j;
    for(i=0;i<=nx-2;i+=2){
		ax[i]=2.0; ax[i+1]=4.0;
		 }
    ax[0]=1.0;ax[nx]=1.0;
    for(j=0;j<=ny-2;j+=2){
		ay[j]=2.0; ay[j+1]=4.0; }
    ay[0]=1.0;ay[ny]=1.0;

    s=0.0; 
    for (i=0;i<=nx;i++) {
	  for (j=0; j<=ny;j++) {
		s=s+ax[i]*ay[j]*ff[i][j];
        } }
    s=s*dx*dy/9.0;
    return s;
 }
int main() 
  {
    double extern ff[mx][my];
    double pi,s,s1,ds,dx,dy,x,y;
    int i,j,k,nx,ny;
    printf("check1");
    pi=atan(1.0)*4.0;
    nx=8;ny=8;
    printf("check1"); s1=4.0/pi/pi;
    for (k=1;k<=6;k++) {
	  nx=2*nx;ny=2*ny;
	  dx=1.0/(double)nx;dy=1.0/(double)ny;
      for (i=0;i<=nx;i++) {
		x=(double)i*dx; 
		for (j=0;j<=ny;j++) {
		y=(double)j*dy;ff[i][j]=f(x,y);} }
  	s=simp2d(ff,dx,dy,nx,ny);	
    ds=fabs(s-s1);
      printf("nx= %4d, ny=%4d, simpson = %20.13e, error = %12.5e\n",nx,ny,s,ds);
     }
 }