#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))
   using control variates MC with control
     phi(x) = 1.0 - cos(x) + 0.5*cos(x)*cos(x).
   Although in principle, integration of the 
   first order term should vanish, this term is
   important for the control. Constant 1.25 is
   the integral over the control:
       1.25 = (1/pi) int(phi,x=0..pi)

   Abramowitz and Stegun give I0(1) = 1.266066

                                    wpp 8/1/2004
*/

   int iN,N,i,j;
   double x,cu,cmu,fmcv,asval;
   double ggl();
   static double seed=331.0;
   double Ex[M],Eav,VEx;
 
   asval = 1.266065878;
   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++){
         fmcv = 0.0;
         for(i=0;i<N;i++){
            x     = ggl(&seed);
            cu    = cos(pi*x);
            cmu   = exp(-cu);
            fmcv += cmu - 1.0 + cu - 0.5*cu*cu;
         }
/* (1/pi) int(phi,0..pi) = 1.25, the integral of the control */
         Ex[j] = 1.25 + fmcv/((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)));
}