#include <stdio.h>
#include <math.h>
#define M 100
#define pi 3.14159265358979323
main()
{
/* 
   Computes integral
     I0(1) = (1/pi) int(z=0..pi) dz exp(-cos(z))
   by raw MC. Abramowitz and Stegun give I0(1) = 1.266066
 
                             wpp 8/1/2004
*/

   int iN,N,i,j;
   double u,asval=1.266065878;
   double ggl();
   static double seed=331.0;
   double Ex[M],Eav,VEx;
 
   printf(" A and S tables:  I0(1) = %12.6e\n",asval);
   printf(" sample     variance        MC I0 val \n");
   printf(" ------   ------------    ------------ \n");
   for(iN=1;iN<6;iN++){
      N = pow(10,iN);
      for(j=0;j<M;j++){
         Ex[j] = 0.0;
         for(i=0;i<N;i++){
            u      = cos(pi*ggl(&seed));
            Ex[j] += exp(-u);
         }
         Ex[j] /= (double)N;
      }
      Eav = 0.0;
      VEx = 0.0;
      for(j=0;j<M;j++){
         Eav += Ex[j];
         VEx += Ex[j]*Ex[j];
      }
      Eav /= (double) M;
      VEx /= (double) M;
      VEx -= Eav*Eav;

      printf("%6d %15.6e %15.6e\n",N,VEx,Eav);
   }
}
#include <math.h>
double ggl(double *ds)
{
/* generate u(0,1) distributed random numbers.
   Seed ds must be saved between calls. ggl is
   essentially the same as the IMSL routine RNUM.

   W. Petersen and M. Troyer, 24 Oct. 2002, ETHZ */

   double t,d2=0.2147483647e10;
   t   = (double) *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = (double) t;
   return((double)((t-1.0e0)/(d2-1.0e0)));
}