#include <conio.h>        /* 此頭函數請不要刪除 */
#include <math.h>

float funI(x,y)                     /* input integrated function I=f(x,y) */
   float x,y;
  {float I;
   I=(6+2.5*x*y+3*x*x*y*y+2*x*x*x)*exp(-2.5*x*x-3*y*y);
   return(I);
  }

float solution_ab(funI,a,b,c,N)    /* child function,solute integration for y. */
        float funI();
        float a,b,c;
        int N;
  {int R;
   int M;
   float x[2001];
   float H;  
   float result_ab;
   float f_ac;
   float f_bc;
   float sum_2M,sum_4M;
  
   M=N/2;
   H=(b-a)/N;
   f_ac=(*funI)(a,c);            /* solute f(a,c). */
   f_bc=(*funI)(b,c);            /* solute f(b,c). */
   for(R=0;R<=N;R++)
     x[R]=a+R*H;
   sum_2M=0;
   for(R=1;R<=M-1;R++)                 /* for R=[1,M-1] ,solute sum of f(x[2*R],c). */
     sum_2M+=(*funI)(x[2*R],c);
   sum_4M=0;
   for(R=1;R<=M;R++)                   /* for R=[1,M] ,solute sum of f(x[2*R-1],c). */
     sum_4M+=(*funI)(x[2*R-1],c);
   result_ab=0;
   result_ab=f_ac+f_bc+2*sum_2M+4*sum_4M;   /* solute integration for y. return(result_ab) */

   return(result_ab);
  }

int solution_cd(solution_ab,a,b,c,d,N,n)   /* child function,solute integration for x. */
       float solution_ab();
       float a,b,c,d;
       int N,n;
  {int m;
   int k; 
   float f_abc;
   float f_abd;
   float f_aby2k;
   float f_aby2k1;
   float h;
   float y[2001];
   float sumf_aby2k;
   float sumf_aby2k1;
   float result_xy;   
   
   m=n/2;
   h=(d-c)/n;

   f_abc=0.0;
   f_abd=0.0;
   f_abc=(*solution_ab)(funI,a,b,c,N);    /* solute f(a,b,c) for x. */
   f_abd=(*solution_ab)(funI,a,b,d,N);    /* solute f(a,b,d) for x. */
   for(k=0;k<=n;k++)
     y[k]=c+k*h;
   sumf_aby2k=0;
   for(k=1;k<=m-1;k++)             /* for k=[1,m-1] ,solute sum of solution_ab(y[2*k],N). */
     {f_aby2k=(*solution_ab)(funI,a,b,y[2*k],N);
      sumf_aby2k+=f_aby2k;
     }
   sumf_aby2k1=0;
   for(k=1;k<=m;k++)               /* for k=[1,m] ,solute sum of solution_ab(y[2*k-1],N). */
     {f_aby2k1=(*solution_ab)(funI,a,b,y[2*k-1],N);
      sumf_aby2k1+=f_aby2k1;
     }
   result_xy=0.0;
   result_xy=f_abc+f_abd+2*sumf_aby2k+4*sumf_aby2k1;
                                   /* solute integration for x without count for factor of */
                                   /* H*h/9 . return(result_xy) */
   return(result_xy);
  }

main()
{int N,n;
 float a,b,c,d;
 float h,H;
 float result_tm;
 float result;
 float solution_ab();
 
 a=0;
 b=3;
 c=0;
 d=2;
 N=n=60;               /* N,n are oven int number,both must be smaller */ 
 h=(d-c)/n;             /* than 107,larger than 1.  */
 H=(b-a)/N;
 result_tm=solution_cd(solution_ab,a,b,c,d,N,n);
 result=result_tm*H*h/9;                       /* solute given integration. */  
 printf("RESULT= %8.6f",result);
 
 getch();        /* 此語句請不要刪除*/
}