/* Integral(Simpson's Method) from a to b        */
/* This is a simple version of intsimp in GAUSS. */
new; cls;
print simpsonint(-100,0,1000000,&f);




proc simpsonint(a,b,n,&f);
   local h,s1,s2,s,i,f:proc;
   h=(b-a)/(2*n);
   s1=0;
   i=1;
   do while i<=n;
      s1=s1+f(a+(2*i-1)*h);
      i=i+1;
   endo;
   s2=0;
   i=1;
   do while i<=n-1;
      s2=s2+f(a+2*i*h);
      i=i+1;
   endo;
   s=h/3*(f(a)+4*s1+2*s2+f(b));
   retp(s);
endp;

fn f(x)=pdfn(x);