C
      REAL FUNCTION rnd(seed)
*     混合法を併用した線形合同法の乱数発生プログラム
*     "Numerical Recipes" の第二版の ran1 を基礎にしている
      IMPLICIT NONE
      INTEGER a,m,q,p,n,ndiv,j,k,seed
      REAL rm,rmax
*     2^31 - 1
*     m = a*q + p
      PARAMETER (a = 16807, m = 2147483647, rm = 1.0/m)
      PARAMETER (q = 127773, p = 2836, n = 32, ndiv = 1 + (m-1)/n)
      PARAMETER (rmax = 1.0 - 1.2e-7)
      INTEGER r(n),r0,r1
      SAVE r,r0,r1
      DATA r/n*0/
*     乱数表の初期化
      IF(seed .NE. 0) THEN
        IF(seed .LT. 0) THEN
          r1 = 584287
        ELSE
          r1 = seed
        ENDIF
        DO j = n+8,1,-1
          k = r1/q
          r1 = a*(r1-k*q) - p*k
          IF (r1 .LT. 0) r1 = r1 + m
          IF (j .LE. n) r(j) = r1
        ENDDO
        r0 = r(1)
      END IF
*     初期化を行わないときの開始点
*     オーバーフローすることなく r1 = mod(a*r1,m) を計算する
      k = r1/q
      r1 = a*(r1 - k*q) - p*k
      IF (r1 .LT. 0) r1 = r1 + m
      j = 1 + r0/ndiv
      r0 = r(j)
      r(j) = r1
      rnd = min(rm*r0,rmax)
      END