1 | c subroutines ... for pseudo-random uniform generations
   2 | c                 in the unit interval !
   3 | c
   4 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   5 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   6 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   7 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   8 | c.......................................................................
   9 | c
  10 | c     generateur de nombres aleatoires uniformes entre 0 et 1
  11 | c     theoriquement compatible avec toutes les ecrans de la creation 
  12 | c
  13 | c....................................................................... 
  14 |       real function urand(iy)
  15 |       integer iy
  16 | c.......................................................................
  17 | c     urand is a uniform random number generator based on theory and 
  18 | c     suggestions given in d.e. knuth (1969), vol 2.
  19 | c     the integer iy should be initialized to an arbitrary integer prior
  20 | c     to the first call to urand. the calling program should not alter
  21 | c     the value of iy between subsequent calls to urand. values of urand
  22 | c     will be returned in the interval (0,1).
  23 | c.......................................................................
  24 |       integer ia,ic,itwo,m2,m,mic
  25 |       double precision halfm
  26 |       real s
  27 |       double precision datan,dsqrt
  28 |       data m2/0/,itwo/2/ 
  29 |       save
  30 |       if (m2.ne.0) go to 20
  31 | c
  32 | c     if first entry, compute machine integer word length
  33 | c
  34 |       m=1
  35 | 10    m2=m
  36 |       m=itwo*m2
  37 |       if (m.gt.m2) go to 10
  38 |       halfm=m2
  39 | c
  40 | c     compute multiplier and increment for linear congruential method
  41 | c
  42 |       ia=8*idint(halfm*datan(1.d0)/8.d0)+5
  43 |       ic=2*idint(halfm*(0.5d0-dsqrt(3.d0)/6.d0))+1
  44 |       mic=(m2-ic)+m2
  45 | c
  46 | c     s is the scale factor for converting to floating point
  47 | c
  48 |       s=0.5/halfm
  49 | c
  50 | c     compute next random number
  51 | c
  52 | 20    iy=iy*ia
  53 | c
  54 | c     the following statement is for computers which do not allow
  55 | c     integer overflow on addition 
  56 | c
  57 |       if (iy.gt.mic) iy=(iy-m2)-m2
  58 | c
  59 |       iy=iy+ic
  60 | c
  61 | c     the following statement is for computers where the word length
  62 | c     for addition is greater than for multiplication
  63 | c
  64 |       if (iy/2.gt.m2) iy=(iy-m2)-m2
  65 | c
  66 | c     the following statement is for computers where integer overflow
  67 | c     affects the sign bit
  68 | c
  69 |       if (iy.lt.0) iy=(iy+m2)+m2
  70 |       urand=float(iy)*s
  71 | c.......................................................................
  72 |       return
  73 |       end