*******************************************************************
      PROGRAM XRANDU3D
*
*     This program is written to test
*     the RANDU random number generator
*     in 3D (See Chapter 11:Monte Carlo).
*
      double precision randu
      double precision random_low, random_med, random_high
      open(unit=10,file='output_randu3.dat',status='new')
      random_low=randu()
      random_med=randu()
      do 10 i=1,2500
      random_high=randu()
      write(10,*), random_low, random_med, random_high
      random_low=random_med
      random_med=random_high
 10   continue
*
*     Now view the 3D plot in 'output_randu3.dat'
*     using Matlab.
*
      end

*******************************************************************
      double precision function randu ()
*
* This function is a modification of ranu() given
* in Chapter 11: Monte Carlo.
* Based on good implementation of SURAND. See
* Park & Miller, reference given in Chapter 11.
* Note that the parameters q,r have been modified
* according to randu's a,m.
*
      integer a,m,q,r,seed
      double precision rm
      parameter (a=65539, m=2147483647, q=32766, r=32773, rm=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
      randu=seed*rm
      return
      end
*******************************************************************