| c--------------------------------------------------------------------- |
| double precision function randlc (x, a) |
| c--------------------------------------------------------------------- |
| |
| c--------------------------------------------------------------------- |
| c |
| c This routine returns a uniform pseudorandom double precision number in the |
| c range (0, 1) by using the linear congruential generator |
| c |
| c x_{k+1} = a x_k (mod 2^46) |
| c |
| c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers |
| c before repeating. The argument A is the same as 'a' in the above formula, |
| c and X is the same as x_0. A and X must be odd double precision integers |
| c in the range (1, 2^46). The returned value RANDLC is normalized to be |
| c between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain |
| c the new seed x_1, so that subsequent calls to RANDLC using the same |
| c arguments will generate a continuous sequence. |
| c |
| c This routine should produce the same results on any computer with at least |
| c 48 mantissa bits in double precision floating point data. On 64 bit |
| c systems, double precision should be disabled. |
| c |
| c David H. Bailey October 26, 1990 |
| c |
| c--------------------------------------------------------------------- |
| |
| implicit none |
| |
| double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z |
| parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23, |
| > t46 = t23 ** 2) |
| |
| c--------------------------------------------------------------------- |
| c Break A into two parts such that A = 2^23 * A1 + A2. |
| c--------------------------------------------------------------------- |
| t1 = r23 * a |
| a1 = int (t1) |
| a2 = a - t23 * a1 |
| |
| c--------------------------------------------------------------------- |
| c Break X into two parts such that X = 2^23 * X1 + X2, compute |
| c Z = A1 * X2 + A2 * X1 (mod 2^23), and then |
| c X = 2^23 * Z + A2 * X2 (mod 2^46). |
| c--------------------------------------------------------------------- |
| t1 = r23 * x |
| x1 = int (t1) |
| x2 = x - t23 * x1 |
| |
| |
| t1 = a1 * x2 + a2 * x1 |
| t2 = int (r23 * t1) |
| z = t1 - t23 * t2 |
| t3 = t23 * z + a2 * x2 |
| t4 = int (r46 * t3) |
| x = t3 - t46 * t4 |
| randlc = r46 * x |
| return |
| end |
| |
| |
| |
| c--------------------------------------------------------------------- |
| c--------------------------------------------------------------------- |
| |
| subroutine vranlc (n, x, a, y) |
| |
| c--------------------------------------------------------------------- |
| c--------------------------------------------------------------------- |
| |
| c--------------------------------------------------------------------- |
| c This routine generates N uniform pseudorandom double precision numbers in |
| c the range (0, 1) by using the linear congruential generator |
| c |
| c x_{k+1} = a x_k (mod 2^46) |
| c |
| c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers |
| c before repeating. The argument A is the same as 'a' in the above formula, |
| c and X is the same as x_0. A and X must be odd double precision integers |
| c in the range (1, 2^46). The N results are placed in Y and are normalized |
| c to be between 0 and 1. X is updated to contain the new seed, so that |
| c subsequent calls to RANDLC using the same arguments will generate a |
| c continuous sequence. |
| c |
| c This routine generates the output sequence in batches of length NV, for |
| c convenience on vector computers. This routine should produce the same |
| c results on any computer with at least 48 mantissa bits in double precision |
| c floating point data. On Cray systems, double precision should be disabled. |
| c |
| c David H. Bailey August 30, 1990 |
| c--------------------------------------------------------------------- |
| |
| integer n |
| double precision x, a, y(*) |
| |
| double precision r23, r46, t23, t46 |
| integer nv |
| parameter (r23 = 2.d0 ** (-23), r46 = r23 * r23, t23 = 2.d0 ** 23, |
| > t46 = t23 * t23, nv = 64) |
| double precision xv(nv), t1, t2, t3, t4, an, a1, a2, x1, x2, yy |
| integer n1, i, j |
| external randlc |
| double precision randlc |
| |
| c--------------------------------------------------------------------- |
| c Compute the first NV elements of the sequence using RANDLC. |
| c--------------------------------------------------------------------- |
| t1 = x |
| n1 = min (n, nv) |
| |
| do i = 1, n1 |
| xv(i) = t46 * randlc (t1, a) |
| enddo |
| |
| c--------------------------------------------------------------------- |
| c It is not necessary to compute AN, A1 or A2 unless N is greater than NV. |
| c--------------------------------------------------------------------- |
| if (n .gt. nv) then |
| |
| c--------------------------------------------------------------------- |
| c Compute AN = AA ^ NV (mod 2^46) using successive calls to RANDLC. |
| c--------------------------------------------------------------------- |
| t1 = a |
| t2 = r46 * a |
| |
| do i = 1, nv - 1 |
| t2 = randlc (t1, a) |
| enddo |
| |
| an = t46 * t2 |
| |
| c--------------------------------------------------------------------- |
| c Break AN into two parts such that AN = 2^23 * A1 + A2. |
| c--------------------------------------------------------------------- |
| t1 = r23 * an |
| a1 = aint (t1) |
| a2 = an - t23 * a1 |
| endif |
| |
| c--------------------------------------------------------------------- |
| c Compute N pseudorandom results in batches of size NV. |
| c--------------------------------------------------------------------- |
| do j = 0, n - 1, nv |
| n1 = min (nv, n - j) |
| |
| c--------------------------------------------------------------------- |
| c Compute up to NV results based on the current seed vector XV. |
| c--------------------------------------------------------------------- |
| do i = 1, n1 |
| y(i+j) = r46 * xv(i) |
| enddo |
| |
| c--------------------------------------------------------------------- |
| c If this is the last pass through the 140 loop, it is not necessary to |
| c update the XV vector. |
| c--------------------------------------------------------------------- |
| if (j + n1 .eq. n) goto 150 |
| |
| c--------------------------------------------------------------------- |
| c Update the XV vector by multiplying each element by AN (mod 2^46). |
| c--------------------------------------------------------------------- |
| do i = 1, nv |
| t1 = r23 * xv(i) |
| x1 = aint (t1) |
| x2 = xv(i) - t23 * x1 |
| t1 = a1 * x2 + a2 * x1 |
| t2 = aint (r23 * t1) |
| yy = t1 - t23 * t2 |
| t3 = t23 * yy + a2 * x2 |
| t4 = aint (r46 * t3) |
| xv(i) = t3 - t46 * t4 |
| enddo |
| |
| enddo |
| |
| c--------------------------------------------------------------------- |
| c Save the last seed in X so that subsequent calls to VRANLC will generate |
| c a continuous sequence. |
| c--------------------------------------------------------------------- |
| 150 x = xv(n1) |
| |
| return |
| end |
| |
| c----- end of program ------------------------------------------------ |
| |
