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

! Uses the algorithm in
! Marsaglia, G. and Tsang, W.W. (2000)

! 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


DAGPUNAR PROG
-------------

FUNCTION random_gamma(s, b, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

!     N.B. This version is in `double precision' and includes scaling

!     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
!     CALLS EITHER random_gamma1 (S > 1.0)
!     OR random_exponential (S = 1.0)
!     OR random_gamma2 (S < 1.0).

!     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
!     B = Scale parameter

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

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

! Local parameters
REAL (dp), PARAMETER  :: one = 1.0_dp, zero = 0.0_dp

IF (s <= zero) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
  STOP
END IF

IF (s >= one) THEN
  fn_val = random_gamma1(s, first)
ELSE IF (s < one) THEN
  fn_val = random_gamma2(s, first)
END IF

! Now scale the random variable
fn_val = b * fn_val
RETURN

CONTAINS


FUNCTION random_gamma1(s, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO GAMMA**(S-1)*EXP(-GAMMA),
! BASED UPON BEST'S T DISTRIBUTION METHOD

!     S = SHAPE PARAMETER OF DISTRIBUTION
!          (1.0 < REAL)

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

!     Local variables
REAL (dp)             :: d, r, g, f, x
REAL (dp), SAVE       :: b, h
REAL (dp), PARAMETER  :: sixty4 = 64.0_dp, three = 3.0_dp, pt75 = 0.75_dp,  &
                         two = 2.0_dp, half = 0.5_dp

IF (s <= one) THEN
  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE'
  STOP
END IF

IF (first) THEN                        ! Initialization, if necessary
  b = s - one
  h = SQRT(three*s - pt75)
END IF

DO
  CALL RANDOM_NUMBER(r)
  g = r - r*r
  IF (g <= zero) CYCLE
  f = (r - half)*h/SQRT(g)
  x = b + f
  IF (x <= zero) CYCLE
  CALL RANDOM_NUMBER(r)
  d = sixty4*g*(r*g)**2
  IF (d <= zero) EXIT
  IF (d*x < x - two*f*f) EXIT
  IF (LOG(d) < two*(b*LOG(x/b) - f)) EXIT
END DO
fn_val = x

RETURN
END FUNCTION random_gamma1



FUNCTION random_gamma2(s, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
! GAMMA2**(S-1) * EXP(-GAMMA2),
! USING A SWITCHING METHOD.

!    S = SHAPE PARAMETER OF DISTRIBUTION
!          (REAL < 1.0)

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

!     Local variables
REAL (dp)             :: r, x, w
REAL (dp), SAVE       :: a, p, c, uf, vr, d
REAL (dp), PARAMETER  :: vsmall = EPSILON(one)

IF (s <= zero .OR. s >= one) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
  STOP
END IF

IF (first) THEN                        ! Initialization, if necessary
  a = one - s
  p = a/(a + s*EXP(-a))
  IF (s < vsmall) THEN
    WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
    STOP
  END IF
  c = one/s
  uf = p*(vsmall/a)**s
  vr = one - vsmall
  d = a*LOG(a)
END IF

DO
  CALL RANDOM_NUMBER(r)
  IF (r >= vr) THEN
    CYCLE
  ELSE IF (r > p) THEN
    x = a - LOG((one - r)/(one - p))
    w = a*LOG(x)-d
  ELSE IF (r > uf) THEN
    x = a*(r/p)**c
    w = x
  ELSE
    fn_val = zero
    RETURN
  END IF

  CALL RANDOM_NUMBER(r)
  IF (one-r <= w .AND. r > zero) THEN
    IF (r*(w + one) >= one) CYCLE
    IF (-LOG(r) <= w) CYCLE
  END IF
  EXIT
END DO

fn_val = x
RETURN

END FUNCTION random_gamma2

END FUNCTION random_gamma



PROGRAM demo_gamma
! Demo program to generate a small smaple of random gamma's

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp)  :: shape, scale, rg(100), avge, mean
LOGICAL    :: first
INTEGER    :: i, order(100)

INTERFACE
  FUNCTION random_gamma(s, b, first) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: s, b
    LOGICAL, INTENT(IN)    :: first
    REAL (dp)              :: fn_val
  END FUNCTION random_gamma

  RECURSIVE SUBROUTINE quick_sort(list, order)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), DIMENSION (:), INTENT(IN OUT)  :: list
    INTEGER, DIMENSION (:), INTENT(OUT)      :: order
  END SUBROUTINE quick_sort
END INTERFACE

WRITE(*, '(a)', ADVANCE='NO') ' Enter the shape & scale parameter values: '
READ(*, *) shape, scale

first = .TRUE.
DO i = 1, 100
  rg(i) = random_gamma(shape, scale, first)
  first = .FALSE.
END DO

CALL quick_sort(rg, order)
WRITE(*, '(" ", 10f7.3)') rg

mean = shape * scale
avge = SUM( rg ) / 100._dp
WRITE(*, '(a, g12.5, a, g12.5)')  &
         ' Population mean = ', mean, '  Sample mean = ', avge

STOP
END PROGRAM demo_gamma



RECURSIVE SUBROUTINE quick_sort(list, order)

! Quick sort routine from:
! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to
! Fortran 90", McGraw-Hill  ISBN 0-07-000248-7, pages 149-150.
! Modified by Alan Miller to include an associated integer array which gives
! the positions of the elements in the original order.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), DIMENSION (:), INTENT(IN OUT)  :: list
INTEGER, DIMENSION (:), INTENT(OUT)       :: order

! Local variable
INTEGER :: i

DO i = 1, SIZE(list)
  order(i) = i
END DO

CALL quick_sort_1(1, SIZE(list))
RETURN


CONTAINS


RECURSIVE SUBROUTINE quick_sort_1(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL (dp)           :: reference, temp
INTEGER, PARAMETER  :: max_simple_sort_size = 6

IF (right_end < left_end + max_simple_sort_size) THEN
  ! Use interchange sort for small lists
  CALL interchange_sort(left_end, right_end)

ELSE
  ! Use partition ("quick") sort
  reference = list((left_end + right_end)/2)
  i = left_end - 1
  j = right_end + 1

  DO
    ! Scan list from left end until element >= reference is found
    DO
      i = i + 1
      IF (list(i) >= reference) EXIT
    END DO
    ! Scan list from right end until element <= reference is found
    DO
      j = j - 1
      IF (list(j) <= reference) EXIT
    END DO


    IF (i < j) THEN
      ! Swap two out-of-order elements
      temp = list(i)
      list(i) = list(j)
      list(j) = temp
      itemp = order(i)
      order(i) = order(j)
      order(j) = itemp
    ELSE IF (i == j) THEN
      i = i + 1
      EXIT
    ELSE
      EXIT
    END IF
  END DO

  IF (left_end < j) CALL quick_sort_1(left_end, j)
  IF (i < right_end) CALL quick_sort_1(i, right_end)
END IF

RETURN
END SUBROUTINE quick_sort_1


SUBROUTINE interchange_sort(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL (dp)           :: temp

DO i = left_end, right_end - 1
  DO j = i+1, right_end
    IF (list(i) > list(j)) THEN
      temp = list(i)
      list(i) = list(j)
      list(j) = temp
      itemp = order(i)
      order(i) = order(j)
      order(j) = itemp
    END IF
  END DO
END DO

RETURN
END SUBROUTINE interchange_sort

END SUBROUTINE quick_sort