ran1.f

网络带宽测试工具

      Function ran1(idum)                    Result( ran_out )
      Use         numerics
      Implicit    None
      Integer  :: idum
! ----------------------------------------------------------------------
! --- ran1 returns a uniform deviate in (0,1).

! --- The algorithm is taken from Press & Teukolsky et.al. and
!     based on the linear congruential method  with choices for
!     M, IA, and IC that are given by D. Knuth in "Semi-numerical
!     algorithms.

! --- Input-parameters:
!     Integer - idum. When idum < 0 the sequence of random values
!                     When idum >= 0, ran1 returns the next value
!                     in the sequence. When ran1 is called for
!                     the first time it is also initialised.

! --- Output-parameters:
!     Integer              - idum. Next value of seed produced by ran1.
!     Real (Kind = Double) - ran_out: the random number returned.
! ----------------------------------------------------------------------

      Real(l_)           :: r(97), ran_out
      Integer            :: iff, x1, x2, x3, j
!
! --- Definitions of the three linear congruences used in generating
!     the random number.
!
      Integer, Parameter :: m1 = 259200, m2 = 134456, m3 = 243000,
     &                      a1 =   7141, a2 =   8121, a3 =   4561,
     &                      c1 =  54773, c2 =  28411, c3 =  51349
      Real*8, Parameter  :: one = 1.0, rm1 = one/m1, rm2 = one/m2 
!
      Save             iff, r, x1, x2, x3
      Data             iff/0/
! ----------------------------------------------------------------------
 
! --- (Re)initialise if required.
 
      If ( idum < 0 .OR. iff == 0 ) Then
         iff = 1

! --- Seed first generator.
	 
	 x1 = Mod( c1 - idum,  m1 ) 
	 x1 = Mod( a1*x1 + c1, m1 )

! --- Use it to seed the second generator.

         x2 = Mod( x1, m2 )
         x1 = Mod( a1*x1 + c1, m1 )

! --- Use generator 1 again to seed generator 3.

         x3 = Mod( a1*x1, m3 )

! --- Now fill array with random values, using gen. 2 for the high
!     order bits and gen. 1 for the low order bits.

         Do j = 1,97
            x1 = Mod( a1*x1 + c1, m1 )
	    x2 = Mod( a2*x2 + c2, m2 )
            r(j) = ( Dble( x1 ) + Dble( x2 ) * rm2 ) * rm1
         End Do
         idum = 1
      End If

! --- This section is only reached when no (re)initialisation takes
!     place. A new random number is generated to fill the place of
!     a randomly picked element from array 'r' (the selection of the
!     index is done by gen. 3).

      x1 = Mod( a1*x1 + c1, m1 )
      x2 = Mod( a2*x2 + c2, m2 )
      x3 = Mod( a3*x3 + c3, m3 )
      j = 1 + ( 97 + x3 )/m3
      ran_out = r(j)
      r(j) = ( Real( x1, l_ ) + Real( x2, l_ ) * rm2 ) * rm1
! ----------------------------------------------------------------------
      End Function ran1