examples/smpi/NAS/common/randi8_safe.f

   1       double precision function randlc(x, a)
   2
   3 c---------------------------------------------------------------------
   4 c
   5 c   This routine returns a uniform pseudorandom double precision number in the
   6 c   range (0, 1) by using the linear congruential generator
   7 c
   8 c   x_{k+1} = a x_k  (mod 2^46)
   9 c
  10 c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  11 c   before repeating.  The argument A is the same as 'a' in the above formula,
  12 c   and X is the same as x_0.  A and X must be odd double precision integers
  13 c   in the range (1, 2^46).  The returned value RANDLC is normalized to be
  14 c   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
  15 c   the new seed x_1, so that subsequent calls to RANDLC using the same
  16 c   arguments will generate a continuous sequence.
  17
  18       implicit none
  19       double precision x, a
  20       integer*8 Lx, La, a1, a2, x1, x2, xa
  21       double precision d2m46
  22       parameter(d2m46=0.5d0**46)
  23
  24       Lx = x
  25       La = A
  26       a1 = ibits(La, 23, 23)
  27       a2 = ibits(La, 0, 23)
  28       x1 = ibits(Lx, 23, 23)
  29       x2 = ibits(Lx, 0, 23)
  30       xa = ishft(ibits(a1*x2+a2*x1, 0, 23), 23) + a2*x2
  31       Lx   = ibits(xa,0, 46)
  32       x    = dble(Lx)
  33       randlc = d2m46*x
  34       return
  35       end
  36
  37
  38 c---------------------------------------------------------------------
  39 c---------------------------------------------------------------------
  40
  41
  42       SUBROUTINE VRANLC (N, X, A, Y)
  43       implicit none
  44       integer n, i
  45       double precision x, a, y(*)
  46       integer*8 Lx, La, a1, a2, x1, x2, xa
  47       double precision d2m46
  48       parameter(d2m46=0.5d0**46)
  49
  50       Lx = X
  51       La = A
  52       a1 = ibits(La, 23, 23)
  53       a2 = ibits(La, 0, 23)
  54       do i = 1, N
  55          x1 = ibits(Lx, 23, 23)
  56          x2 = ibits(Lx, 0, 23)
  57          xa = ishft(ibits(a1*x2+a2*x1, 0, 23), 23) + a2*x2
  58          Lx   = ibits(xa,0, 46)
  59          y(i) = d2m46*dble(Lx)
  60       end do
  61       x = dble(Lx)
  62       return
  63       end
  64