(list of modules)(lisf of routines)
    1   module knuth
    2   
    3   ! This module implements the Fortran 90 version of D. E. Knuth
    4   ! prng.
    5   
    6   implicit none
    7   
    8   private
    9   
   10   public :: rRan, iRan, iRanRange, Randomize, ModuleName
   11   
   12   integer,   parameter :: dpr = selected_real_kind(15,200)
   13   integer,   parameter :: kk = 100
   14   integer,   parameter :: ll = 37
   15   integer,   parameter :: mm = 2**30
   16   integer,   parameter :: tt = 70
   17   integer,   parameter :: kkk = kk + kk - 1
   18   integer,   parameter :: nn = 1024 ! generate RNs nn at a time
   19   
   20   character(*), parameter :: ModuleName = 'knuth'
   21   
   22   real(dpr), parameter :: ulp = 1.d0 / ( 2.d0**52 )
   23   
   24   real(dpr), dimension(kk) :: ranx
   25   integer,   dimension(kk) :: iranx
   26   
   27   real(dpr), dimension(nn) :: internal_rns
   28   integer,   dimension(nn) :: internal_irns
   29   integer :: ipos, rpos
   30   
   31   contains
   32   
   33   !--------------------------------------------------------------------------------
   34   subroutine Randomize( seed, rank )
   35   !--------------------------------------------------------------------------------
   36     integer,           intent(inout) :: seed
   37     integer, optional, intent(in)    :: rank
   38     integer, dimension(8) :: v
   39   
   40     if ( seed == 0 ) then
   41        call date_and_time ( values = v )
   42        seed = 1000 * v(8) + 100 * v(7) + 10 * v(6) + v(8)*v(7)
   43     end if
   44   
   45     if ( present( rank ) ) seed = 1000 * seed / (3*rank + 1)
   46   
   47     call RNFSTR( seed )
   48     call rnstrt( seed )
   49     ipos = nn+1
   50     rpos = nn+1
   51   
   52   end subroutine Randomize
   53   
   54   
   55   
   56   !--------------------------------------------------------------------------------
   57   function rRan( )
   58   !--------------------------------------------------------------------------------
   59   
   60     real(dpr) :: rRan
   61   
   62     if ( rpos > nn ) then
   63        call rnfarr( internal_rns, nn )
   64        rpos = 1
   65     end if  
   66     rRan = internal_rns( rpos )
   67     rpos = rpos + 1
   68   end function rRan
   69   
   70   
   71   !--------------------------------------------------------------------------------
   72   function iRan( )
   73   !--------------------------------------------------------------------------------
   74     integer :: iRan
   75     if ( ipos > nn ) then
   76        call rnarry( internal_irns, nn )
   77        ipos = 1
   78     end if
   79     iRan = internal_irns( ipos )
   80     ipos = ipos + 1
   81   end function iRan
   82   
   83   
   84   !--------------------------------------------------------------------------------
   85   function iRanRange( low, high )
   86   !--------------------------------------------------------------------------------
   87     integer, intent(in) :: low, high
   88     integer :: iRanRange
   89   
   90     iRanRange = low + floor( (high-low+1._dpr) * rRan() )
   91   
   92   end function iRanRange
   93   
   94   
   95   
   96   !--------------------------------------------------------------------------------
   97   subroutine rnfarr( aa, n )
   98   !--------------------------------------------------------------------------------
   99   !       FORTRAN 77 version of "ranf_array"
  100   !       from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
  101   !       including the MODIFICATIONS made in the 9th printing (2002)
  102   !       ********* see the book for explanations and caveats! *********
  103   
  104     real(dpr), dimension(:), intent(inout) :: aa
  105     integer,                 intent(in)    :: n
  106   
  107     integer   :: j
  108     real(dpr) :: y
  109   
  110     do j = 1, kk
  111        aa(j) = ranx(j)
  112     end do
  113     ! do 2
  114     do j = kk + 1, n
  115        y = aa(j-kk) + aa(j-ll)
  116        aa(j) = y - idint(y)
  117     end do
  118     ! do 3
  119     do j = 1, ll
  120        y = aa(n+j-kk) + aa(n+j-ll)
  121        ranx(j) = y - idint(y)
  122     end do
  123     ! do 4
  124     do j = ll+1, kk
  125        y = aa(n+j-kk) + ranx(j-ll)
  126        ranx(j) = y - idint(y)
  127     end do
  128   
  129   end subroutine rnfarr
  130   
  131   
  132   !--------------------------------------------------------------------------------
  133   subroutine rnfstr( seed )
  134   !--------------------------------------------------------------------------------
  135   
  136     integer, intent(in) :: seed
  137   
  138     integer   :: s, sseed, j, t
  139     real(dpr) :: ss, v
  140     real(dpr), dimension(kkk) :: u
  141     
  142     if ( seed < 0 ) then
  143        sseed = mm - 1 - mod(-1-seed,mm)
  144     else
  145        sseed = mod(seed,mm)
  146     end if
  147     ss = 2._dpr * ulp * (sseed+2)
  148     do j = 1, kk
  149        u(j) = ss
  150        ss = ss + ss
  151        if ( ss >= 1._dpr ) ss = ss - 1._dpr + 2._dpr * ulp
  152     end do
  153     u(2) = u(2) + ulp
  154     s = sseed
  155     t = tt - 1
  156     10 continue
  157     do j = kk, 2, -1
  158        u(j+j-1) = u(j)
  159        u(j+j-2) = 0
  160     end do
  161     do j = kkk, kk+1, -1
  162        v = u(j-(kk-ll)) + u(j)
  163        u(j-(kk-ll)) = v - idint(v)
  164        v = u(j-kk) + u(j)
  165        u(j-kk) = v - idint(v)
  166     end do
  167     if ( mod(s,2) == 1 ) then
  168        do j = kk, 1, -1
  169           u(j+1) = u(j)
  170        end do
  171        u(1) = u(kk+1)
  172        v = u(ll+1) + u(kk+1)
  173        u(ll+1) = v - idint(v)
  174     end if
  175   
  176     if ( s /= 0 ) then
  177        s = s / 2
  178     else
  179        t = t - 1
  180     end if
  181   
  182     if ( t > 0 ) goto 10
  183   
  184     do j = 1, ll
  185        ranx(j+kk-ll) = u(j)
  186     end do
  187   
  188     do j = ll+1, kk
  189        ranx(j-ll) = u(j)
  190     end do
  191     
  192     do j = 1, 10
  193        call rnfarr(u,kkk)
  194     end do
  195   
  196   end subroutine rnfstr
  197   
  198   
  199   !--------------------------------------------------------------------------------
  200   subroutine rnarry( aa, n )
  201   !--------------------------------------------------------------------------------
  202   !       FORTRAN 77 version of "ran_array"
  203   !       from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
  204   !       including the MODIFICATIONS made in the 9th printing (2002)
  205   !       ********* see the book for explanations and caveats! *********
  206   
  207     integer, dimension(:), intent(inout) :: aa
  208     integer,               intent(in)    :: n
  209   
  210     integer :: j
  211   
  212     do j = 1, kk
  213        aa(j) = iranx(j)
  214     end do
  215     do j = kk+1, n
  216        aa(j) = aa(j-kk) - aa(j-ll)
  217        if ( aa(j) < 0 ) aa(j) = aa(j) + mm
  218     end do
  219     do j = 1, ll
  220        iranx(j) = aa(n+j-kk) - aa(n+j-ll)
  221        if ( iranx(j) < 0 ) iranx(j) = iranx(j) + mm
  222     end do
  223     do j = ll+1, kk
  224        iranx(j) = aa(n+j-kk) - iranx(j-ll)
  225        if ( iranx(j) < 0 ) iranx(j) = iranx(j) + mm
  226     end do
  227   
  228   end subroutine rnarry
  229   
  230   
  231   !--------------------------------------------------------------------------------
  232   subroutine rnstrt( seed )
  233   !--------------------------------------------------------------------------------
  234   
  235     integer, intent(inout) :: seed
  236   
  237     integer, dimension(kkk) :: x
  238     integer                 :: j, sseed, ss, t
  239   
  240     IF (SEED .LT. 0) THEN
  241        SSEED=MM-1-MOD(-1-SEED,MM)
  242     ELSE
  243        SSEED=MOD(SEED,MM)
  244     END IF
  245     SS=SSEED-MOD(SSEED,2)+2
  246     DO J=1,KK
  247        X(J)=SS
  248        SS=SS+SS
  249        IF (SS .GE. MM) SS=SS-MM+2
  250     end do
  251     X(2)=X(2)+1
  252     SS=SSEED
  253     T=TT-1
  254   10 continue
  255     DO J=KK,2,-1
  256        X(J+J-1)=X(J)
  257        X(J+J-2)=0
  258     end do
  259     DO J=KKK,KK+1,-1
  260        X(J-(KK-LL))=X(J-(KK-LL))-X(J)
  261        IF (X(J-(KK-LL)) .LT. 0) X(J-(KK-LL))=X(J-(KK-LL))+MM
  262        X(J-KK)=X(J-KK)-X(J)
  263        IF (X(J-KK) .LT. 0) X(J-KK)=X(J-KK)+MM
  264     end do
  265     IF (MOD(SS,2) .EQ. 1) THEN
  266        DO J=KK,1,-1
  267           X(J+1)=X(J)
  268        end do
  269        X(1)=X(KK+1)
  270        X(LL+1)=X(LL+1)-X(KK+1)
  271        IF (X(LL+1) .LT. 0) X(LL+1)=X(LL+1)+MM
  272     END IF
  273     IF (SS .NE. 0) THEN
  274        SS=SS/2
  275     ELSE
  276        T=T-1
  277     END IF
  278     IF (T .GT. 0) GO TO 10
  279     DO J=1,LL
  280        iranx(J+KK-LL)=X(J)
  281     end do
  282     DO J=LL+1,KK
  283        iranx(J-LL)=X(J)
  284     end do
  285     DO J=1,10
  286        CALL RNARRY(X,KKK)
  287     end do
  288     
  289   end subroutine rnstrt
  290   
  291   
  292   
  293   end module knuth
  294   
  295   
  296   !!$PROGRAM MAIN
  297   !!$!      a rudimentary test program:
  298   !!$  use knuth
  299   !!$  IMPLICIT none
  300   !!$  real(dpr), dimension(2009) :: A
  301   !!$  integer :: i
  302   !!$  CALL RNFSTR(310952)
  303   !!$  DO I=1,2010
  304   !!$     CALL RNFARR(A,1009)
  305   !!$  end do
  306   !!$  PRINT '(F22.20,a,F22.20)',A(1), ' should be ', 0.36410514377569680455_dpr
  307   !!$!                   the number should be 0.36410514377569680455
  308   !!$  CALL RNFSTR(310952)
  309   !!$  DO I=1,1010
  310   !!$          CALL RNFARR(A,2009)
  311   !!$  end do
  312   !!$  PRINT '(F22.20,a)',A(1), ' (should be the same)'
  313   !!$!                                 again, 0.36410514377569680455
  314   !!$END PROGRAM MAIN
List of routines:
Randomize
rnarry
rnfarr
rnfstr
rnstrt
iRan
iRanRange
rRan
(list of modules)