FUNCTION ran_gamma(s, first) RESULT(fn_val)
USE random
IMPLICIT NONE

! Uses the algorithm in
! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
! gamma variables', Trans. om Math. Software (TOMS).

! Generates a random gamma deviate for shape parameter s >= 1.

REAL, INTENT(IN)    :: s
LOGICAL, INTENT(IN) :: first
REAL                :: fn_val

! Local variables
REAL, SAVE  :: c, d
REAL        :: u, v, x

IF (s < 1.0) THEN
   WRITE(*, *) 'Shape parameter must be >= 1'
   STOP
END IF

IF (first) THEN
  d = s - 1./3.
  c = 1.0/SQRT(9.0*d)
END IF

! Start of main loop
DO

! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.

  DO
    x = random_normal()
    v = (1.0 + c*x)**3
    IF (v > 0.0) EXIT
  END DO

! Generate uniform variable U

  CALL RANDOM_NUMBER(u)
  IF (u < 1.0 - 0.0331*x**4) THEN
    fn_val = d*v
    EXIT
  ELSE IF (LOG(u) < 0.5*x**2 + d*(1.0 - v + LOG(v))) THEN
    fn_val = d*v
    EXIT
  END IF
END DO

RETURN
END FUNCTION ran_gamma



PROGRAM test_rgamma
USE random
IMPLICIT NONE

INTERFACE
  FUNCTION ran_gamma(s, first) RESULT(fn_val)
    USE random
    IMPLICIT NONE
    REAL, INTENT(IN)    :: s
    LOGICAL, INTENT(IN) :: first
    REAL                :: fn_val
  END FUNCTION ran_gamma
END INTERFACE

REAL, PARAMETER  :: shape(8) = (/ 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 20.0, 50.0 /)
INTEGER          :: i, j
REAL             :: aver, dev, sumsq, start, finish, x

DO i = 1, 8
  WRITE(*, *)
  WRITE(*, '(a, f7.1)') ' Shape parameter = ', shape(i)

! First time the new algorithm

  CALL CPU_TIME(start)
  x = ran_gamma(shape(i), .TRUE.)
  aver = x
  sumsq = 0.0
  DO j = 2, 1000000
    x = ran_gamma(shape(i), .FALSE.)
    dev = x - aver
    aver = aver + dev/j
    sumsq = sumsq + dev*(x-aver)
  END DO
  CALL CPU_TIME(finish)
  sumsq = SQRT(sumsq/999999.)
  WRITE(*, *) 'Using Marsaglia & Tsang algorithm'
  WRITE(*, '(a, f7.2, a, 2f7.3)')  &
        ' Time = ', finish - start, 'sec.   Mean & st.devn. = ', aver, sumsq

! Repeat using algorithm in module random

  CALL CPU_TIME(start)
  x = random_gamma(shape(i), .TRUE.)
  aver = x
  sumsq = 0.0
  DO j = 2, 1000000
    x = random_gamma(shape(i), .FALSE.)
    dev = x - aver
    aver = aver + dev/j
    sumsq = sumsq + dev*(x-aver)
  END DO
  CALL CPU_TIME(finish)
  sumsq = SQRT(sumsq/999999.)
  WRITE(*, *) 'Using algorithm in module random'
  WRITE(*, '(a, f7.2, a, 2f7.3)')  &
        ' Time = ', finish - start, 'sec.   Mean & st.devn. = ', aver, sumsq
END DO

STOP
END PROGRAM test_rgamma