C                Code                                                 FORTRAN                Code
#include <math.h>                                 /* math library header file */  
#include<stdio.h>                                  /* io header file  ( eg ., printf) */  
#include <stdlib.h>                               /* system library header file  
                                                                    * (eg.,  malloc ) */  

#define max(a,b)     ((a)>(b)?(a):(b))     /* counterpart of fortran  
                                                                            * max function */  
int seed = 1;                                           /*global variables:seed cseed*/  
extern unsigned long long int cseed = 1; 

      
double ranu1 (void){  
 int a,m ;  
 double minv;  
 extern   int seed;                             /* indicate a  global variable  */  
                                                            /* would be referred outside the*/  
a=16807;                                         /* current scope. */  
 m=2147483647;  
 minv = (double) 1.0/m;  

 /* seed = a*seed%m; */ 
/* (%: mod operation .)  intuitive implentation, may cause   
  overflow. */ 
seed=((unsigned int )a*seed)% m ;   
/* circumvent the overflow problem. */ 
 return (double) seed*minv;  
}  
 

        double precision function ranu ()  
c Naive implementation of SURAND. Subroutine ranset should be called  
c (once) before the first function call for initialization.  
      integer a, m, seed  
      double precision minv  
      parameter (a=16807, m=2147483647, minv=1d0/m)  
      common /random/ seed  
      save /random/  
      data seed /1/  
      seed = mod ( a * seed, m)  
      ranu = seed * minv  
      return  
      end 
double ranu2 (void){  
  int   a=16807, m=2147483647, q=127773, r=2836;  
 double minv = (double) 1.0/m;  
        extern  int seed;  

 seed = a*(seed%q)-r*(seed/q);  

     if (seed < 0) seed = seed + m;  
 return (double) seed*minv;  
}  
 

       double precision function ranu ()  
c Good implementation of SURAND. See Park & Miller, Comm.  
c ACM 31:1192, 1988.  
c Subroutine ranset should be called (once) before the first function call.  
      integer a, m, q, r, seed  
      double precision minv  
      parameter (a=16807, m=2147483647, q=127773, r=2836, minv=1d0/m)  
      common /random/ seed  
      save /random/  
      data seed /1/  
      seed = a * mod(seed, q) - r * (seed/q)  
      if (seed .le. 0) seed = seed + m  
      ranu = seed * minv  
      return  
      end 
float  matgen ( float **amat, int n )      /* preallocated 2D array */  
 {                                                                   /* pointer amat. */  
                                                                     /* norma is returned as  */  
                                                                     /*function     value */  
  
int a = 3125, m=65536;  
 float norma, halfm = 32768.0, quartm=16384.0;  
 int iseed = 1325;  
 int i,j;  

 norma = 0.0;  
 for (j=0;j<n;j++){                                   /* outer loop begin */  
  for (i=0;i<n;i++){                                 /* inner loop begin */  
   iseed = a*iseed%m;  
   amat [i][j] = (float) (iseed - halfm) / quartm ;  
   norma = max (amat[i][j], norma);  
  }                                                            /*  inner loop end */  
 }                                                              /* outer loop end */  
 return norma;  
}

       subroutine matgen(a,Lda,n,b,norma)  
      real amat(Lda,1), b(1), norma, halfm, quartm  
      integer a, m  
      parameter  (a = 3125, m = 65536)  
      parameter  (halfm = 32768.0, quartm = 16384.0)  
  
      iseed = 1325  
      norma = 0.0  
      do 30 j = 1, n  
         do 20 i = 1, n  
            iseed = mod (a * iseed, m)  
            amat(i,j) = (iseed - halfm) / quartm  
            norma  = max ( amat(i,j), norma)  
   20    continue  
   30 continue  
      return  
      end
void ranuv (int n, double *vec){  
 int   a=16807, m=2147483647, q=127773, r=2836;  
 double minv = (double) 1.0/m;  
 int i;  
        extern int seed;  

 if (n < 1) return;  
 for (i=0;i<n;i++){  
  seed= a*(seed%q)-r*(seed/q);  
  if (seed < 0 ) seed = seed + m;  
  vec[i]=(double) seed*minv;  
 }  
 return;  
}  
 

        subroutine ranuv (n, vec)  
c Generate a vector of n pseudorandom uniform variates  
      integer n, a, m, q, r, seed  
      double precision vec(n), minv  
      parameter (a=16807, m=2147483647, q=127773, r=2836, minv=1d0/m)  
      common /random/ seed  
      save /random/  
      if (n .lt. 1) return  
      do 10 i = 1, n  
         seed = a * mod(seed, q) - r * (seed/q)  
         if (seed .le. 0) seed = seed + m  
         vec(i) = seed * minv  
   10 continue  
      return  
      end
void rannv1 (int n, double *vec, double mean, double *var){  
 double temp ;  
 double p0=-.322232431088, p1=-1.0, p2=-.342242088547;  
 double p3=-.204231210245e-1, p4=-.453642210148e-4;  
 double q0=.99348462606e-1, q1=.588581570495, q2=.531103462366;  
 double q3=.10353775285, q4=.38560700634e-2;  
         int i;  

 if (n < 1) return;  
 ranuv (n,&vec[0]);  
 for (i=0;i<n;i++){  
  temp = vec[i];  
  if (temp > 0.5 ) vec[i] = 1.0 -vec[i];  
  vec[i] = sqrt(log(1.0/(vec[i]*vec[i])));  
  vec[i]= vec[i]+ ((((vec[i]*p4 + p3)*vec[i] + p2)*  
         vec[i] + p1) * vec[i] + p0)/  
        ((((vec[i] * q4 + q3) * vec[i] + q2) *  

          vec[i] + q1) * vec[i] + q0);  
  if (temp > 0.5 )  vec[i] = -vec[i];  
 }  

 for (i=0;i<n;i++){  
  vec[i] = sqrt(var[i]) * vec[i] + mean;  
 }  
 return;  

}  
 

      subroutine rannv1 (n, vec, mean, var)  
c A vector of n pseudorandom numbers is generated, each from a standard  
c normal distribution (mean zero, variance one), based on Odeh and Evans,  
c App. Stat. 23:96 (1974). For a nonzero mean MU and/or non unity variance,  
c set  vec(i) = mu + sqrt(sigma(i))*vec(i).  
c Subroutine ranset should be called once before the first subroutine call.  
      integer n  
      double precision vec(n),mean,var(n),  
     *      temp,p0,p1,p2,p3,p4,q0,q1,q2,q3,q4  
      parameter (p0=-.322232431088d0, p1=-1d0, p2=-.342242088547d0,  
     *      p3=-.204231210245d-1, p4=-.453642210148d-4,  
     *      q0=.99348462606d-1, q1=.588581570495d0, q2=.531103462366d0,  
     *      q3=.10353775285d0, q4=.38560700634d-2)  
      if (n .lt. 1) return  
      call ranuv(n, vec)  
      do 10 i = 1, n  
         temp = vec(i)  
         if (temp .gt. 0.5d0) vec(i) = 1d0 - vec(i)  
         vec(i) = sqrt(log(1d0/vec(i)**2))  
         vec(i) = vec(i) +  
     *         ((((vec(i) * p4 + p3) * vec(i) + p2) *  
     *         vec(i) + p1) * vec(i) + p0) /  
     *         ((((vec(i) * q4 + q3) * vec(i) + q2) *  
     *         vec(i) + q1) * vec(i) + q0)  
         if (temp .lt. 0.5d0) vec(i) = -vec(i)  
   10 continue  
      do 20 i = 1, n  
         vec(i) = sqrt(var(i)) * vec(i) + mean  
   20 continue  
      return  
      end
void s_crand( int iseed){                             /*  fortran rand  function return */  
extern  unsigned long long int cseed;          /*  a value between [0,1 ),  c  rand  */  
 cseed = iseed;                                              /*  function return  a  integer */  
 return;  
}  

double crand (void){                                   /* this crand  function  provides   */  
unsigned long long int    M,ia,ic;               /* identical random values between */  
              int    shfM,shfc ;                           /* [0,1) */  
M  = 2147483648;  
ia = 1103515245;  
ic = 12345;  
            cseed= (ia*cseed+ic)%M;  
            shfM=       M>>16;                         /* >>:  right shift ,  counterpart */  
            shfc=   cseed>>16;                          /* of fortran rshift function       */  
            return  (double) shfc/ (double)shfM ;  
  
}

  
void monte ( int nstep){  
 int i,nin ;  
        extern int seed;  
 double  x, y,tmp;  
 nin = 0;  
  
 s_crand(seed);  
 for (i=0;i<nstep;i++){  
  x = crand();  
  y = crand();  
  tmp = sqrt(x*x + y*y);  
  if (tmp < 1.0)  nin = nin + 1;  
 }  
 printf("%d  %lf\n",nstep, (double) (4.0 * nin)/nstep);  
 return;  
}  
  
  
 		       subroutine monte(nstep)  
      implicit none  
      integer nstep, i, nin, iseed  
      double precision  x, y, tmp, rand  

      nin = 0  
      iseed = 12345  
      call srand(iseed)  

      do 30 i = 1, nstep  
         x = rand()  
         y = rand()  
         tmp = sqrt(x*x + y*y)  
         if (tmp .lt. 1.d0) then  
            nin = nin + 1  
         endif  
   30 continue  

      print *, nstep, (4.d0 * nin)/nstep  
      return  
      end