CBelow is the f77 source code for the random number generator
Cwe all know and love as rannyu.  The gererator is seeded by
Cthe f77 call
C
C      call setrn(iseed)
C
Cwhere iseed is dimensioned intger iseed(4).  Tor retrieve the
Ccurrent seed (useful for continuing a series of runs), the call
Cis
C      call savern(iseed)
C
CThe random number is returned as a function
C
C     x = rannyu(0)
C
Cwhich returns a floating point number on the interval (0,1).
CLet me know if you have any problems.
C
C
      function rannyu(x)
      common /rnyucm/ m1,m2,m3,m4,l1,l2,l3,l4
c 
c     multiplier is represented as m1*t**3 + m2*t**2 + m3*t + m4
c     where t=2**12.  the seed is similarly stored.
c 
c     rannyu generates pseudo-random numbers by
c 
c     r <== mod(r*11**13,2**48)
c 
c     and returns a number uniform on (0,1).
c 
      i1=l1*m4+l2*m3+l3*m2+l4*m1
      i2=l2*m4+l3*m3+l4*m2
      i3=l3*m4+l4*m3
      i4=l4*m4
      l4=mod(i4,2**12)
      i3=i3+i4/2**12
      l3=mod(i3,2**12)
      i2=i2+i3/2**12
      l2=mod(i2,2**12)
      l1=mod(i1+i2/2**12,2**12)
      rannyu=2.0**(-12)*(float(l1)+
     +       2.0**(-12)*(float(l2)+
     +       2.0**(-12)*(float(l3)+
     +       2.0**(-12)*(float(l4)))))
      return
      end
      subroutine ranyuv(a,n,j)
c
c   this routine will return a vector of random numbers
c
      common /rnyucm/ m1,m2,m3,m4,l1,l2,l3,l4
      dimension a(1)
      do 99 i=1,n,j
      i1=l1*m4+l2*m3+l3*m2+l4*m1
      i2=l2*m4+l3*m3+l4*m2
      i3=l3*m4+l4*m3
      i4=l4*m4
      l4=mod(i4,2**12)
      i3=i3+i4/2**12
      l3=mod(i3,2**12)
      i2=i2+i3/2**12
      l2=mod(i2,2**12)
      l1=mod(i1+i2/2**12,2**12)
   99 a(i) = 2.0**(-12)*(float(l1)+
     +       2.0**(-12)*(float(l2)+
     +       2.0**(-12)*(float(l3)+
     +       2.0**(-12)*(float(l4)))))
      return
      end
      subroutine setrn(iseed)
c
c  seed the generator. care is taken to ensure
c  seed is odd. default seed is 1
c
      common /rnyucm/ m(4),l(4)
      integer iseed(4)
      data m / 502,1521,4071,2107/
      data l /   0,   0,   0,   1/
         do 10 i=1,4
         l(i)=iseed(i)
   10    continue
      l(4)=2*(l(4)/2)+1
      return
      end
      subroutine savern(iseed)
c
c  routine to return current value of iseed
c
      common /rnyucm/ m(4),l(4)
      integer iseed(4)
         do 10 i=1,4
         iseed(i)=l(i)
   10    continue
      return
      end