MODULE randdp ! Replaces COMMON /randpp/ and defines dp IMPLICIT NONE INTEGER, PARAMETER, PUBLIC :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), SAVE, PUBLIC :: poly(101), other, offset INTEGER, SAVE, PUBLIC :: index CONTAINS ! Nick Maclaren's double precision random number generator. ! This version, which is compatible with Lahey's ELF90 compiler, ! is by Alan Miller ( alan @ vic.cmis.csiro.au, www.vic.cmis.csiro.au/~alan ) ! Latest revision - 18 December 1997 ! Copyright (C) 1992 N.M. Maclaren ! Copyright (C) 1992 The University of Cambridge ! This software may be reproduced and used freely, provided that all ! users of it agree that the copyright holders are not liable for any ! damage or injury caused by use of this software and that this ! condition is passed onto all subsequent recipients of the software, ! whether modified or not. SUBROUTINE sdprnd (iseed) INTEGER, INTENT(IN) :: iseed ! Local variables REAL (dp) :: x REAL (dp), PARAMETER :: xmod = 1000009711.0_dp, ymod = 33554432.0_dp INTEGER :: ix, iy, iz, i LOGICAL, SAVE :: inital = .TRUE. ! ISEED should be set to an integer between 0 and 9999 inclusive; ! a value of 0 will initialise the generator only if it has not ! already been done. IF (inital .OR. iseed /= 0) THEN inital = .false. ELSE RETURN END IF ! INDEX must be initialised to an integer between 1 and 101 inclusive, ! POLY(1...N) to integers between 0 and 1000009710 inclusive (not all 0), ! and OTHER to a non-negative proper fraction with denominator 33554432. ! It uses the Wichmann-Hill generator to do this. ix = MOD(ABS(iseed), 10000) + 1 iy = 2*ix + 1 iz = 3*ix + 1 DO i = -10,101 IF (i >= 1) poly(i) = AINT(xmod*x) ix = MOD(171*ix, 30269) iy = MOD(172*iy, 30307) iz = MOD(170*iz, 30323) x = MOD(DBLE(ix)/30269.0_dp + DBLE(iy)/30307.0_dp + DBLE(iz)/30323.0_dp, 1.0_dp) END DO other = AINT(ymod*x)/ymod offset = 1.0_dp/ymod INDEX = 1 RETURN END SUBROUTINE sdprnd FUNCTION dprand() RESULT(fn_val) REAL (dp) :: fn_val ! Local variables REAL (dp) :: x, y ! N.B. ymod has been removed from the previous DATA statement; it caused a ! fatal error as it is not used. REAL (dp), PARAMETER :: xmod = 1000009711.0_dp, xmod2 = 2000019422.0_dp, & xmod4 = 4000038844.0_dp, tiny = 1.0E-17_dp, & zero = 0.0_dp, one = 1.0_dp INTEGER :: n LOGICAL, SAVE :: inital = .TRUE. ! This returns a uniform (0,1) random number, with extremely good ! uniformity properties. It assumes that REAL (dp) provides ! at least 33 bits of accuracy, and uses a power of two base. IF (inital) THEN CALL sdprnd (0) inital = .false. END IF ! See [Knuth] for why this implements the algorithm described in the paper. ! Note that this code is tuned for machines with fast REAL (dp), but ! slow multiply and divide; many, many other options are possible. n = INDEX - 64 IF (n <= 0) n = n + 101 x = poly(INDEX) + poly(INDEX) x = xmod4 - poly(n) - poly(n) - x - x - poly(INDEX) IF (x < zero) THEN IF (x < -xmod) x = x + xmod2 IF (x < zero) x = x + xmod ELSE IF (x >= xmod2) THEN x = x - xmod2 IF (x >= xmod) x = x - xmod END IF IF (x >= xmod) x = x - xmod END IF poly(INDEX) = x INDEX = INDEX + 1 IF (INDEX > 101) INDEX = INDEX - 101 ! Add in the second generator modulo 1, and force to be non-zero. ! The restricted ranges largely cancel themselves out. DO y = 37.0_dp*other + offset other = y - AINT(y) IF (other /= zero) EXIT END DO x = x/xmod + other IF (x >= one) x = x - one fn_val = x + tiny RETURN END FUNCTION dprand END MODULE randdp