SUBROUTINE random_Zipf(a, r, setup)
! Generate a random deviate (r) from the Zipf (or Zeta) distribution.
! If setup = .TRUE. (or local variable b < 0), two local variables (b and
! const) are set/reset for future calls, otherwise it is assumed that
! the value of a is unchanged since the last call.

! The value of a must be > 1.0 (NOT CHECKED)

! Algorithm from page 551 of:
! Devroye, Luc (1986) `Non-uniform random variate generation',
! Springer-Verlag: Berlin.   ISBN 3-540-96305-7 (also 0-387-96305-7)

! Programmer: Alan Miller (amiller @ users,bigpond.net.au)
! Latest revision - 3 February 1998

IMPLICIT NONE
REAL, INTENT(IN)        :: a
INTEGER, INTENT(OUT)    :: r
LOGICAL, INTENT(IN OUT) :: setup

! Local variables
REAL, SAVE      :: b = -1.0, const
REAL            :: t, u, v
REAL, PARAMETER :: one = 1.0

IF (b < 0.0 .OR. setup) THEN
  b = 2.0**(a - one)
  const = -one/(a - one)
  setup = .FALSE.
END IF

DO
  CALL RANDOM_NUMBER(u)
  CALL RANDOM_NUMBER(v)
  r = FLOOR(u**const)
  t = (one + one/r)**(a - one)
  IF (v*r*(t - one)/(b - one) <= t/b) EXIT
END DO

RETURN
END SUBROUTINE random_Zipf



PROGRAM t_random_Zipf
IMPLICIT NONE
INTEGER              :: r, i, nfreq, count, i1
INTEGER, ALLOCATABLE :: freq(:)
REAL, ALLOCATABLE    :: expctd(:)
REAL                 :: a, total
REAL, PARAMETER      :: one = 1.0
LOGICAL              :: setup

INTERFACE
  SUBROUTINE random_Zipf(a, r, setup)
    IMPLICIT NONE
    REAL, INTENT(IN)        :: a
    INTEGER, INTENT(OUT)    :: r
    LOGICAL, INTENT(IN OUT) :: setup
  END SUBROUTINE random_Zipf
END INTERFACE

WRITE(*, *) 'Enter a (> 1.0): '
READ(*, *) a
setup = .TRUE.
! 1/20 = 5%   Store individual frequencies up to last 5%
nfreq = NINT( (20.0/(a - one))**(one/(a-one)) )
ALLOCATE( freq(nfreq), expctd(nfreq) )
freq = 0

! Generate 1000 random deviates from the Zipf distribution

DO i = 1, 1000
  CALL random_Zipf(a, r, setup)
  r = MIN(r, nfreq)
  freq(r) = freq(r) + 1
END DO

! Calculate approximate expected frequencies
total = 0.0
DO i = 1, nfreq-1
  expctd(i) = REAL(i)**(-a)
  total = total + expctd(i)
END DO
expctd(nfreq) = one / ( (a-one)*(nfreq - 0.5)**(a-one) )
total = total + expctd(nfreq)
expctd = expctd * (1000. / total)

WRITE(*, *) '    Range      Obs.freq.  Expctd.freq.'
i = 1
DO
  count = 0
  total = 0
  i1 = i
  DO
    count = count + freq(i)
    total = total + expctd(i)
    IF (total > 25. .OR. i >= nfreq) EXIT
    i = i + 1
  END DO
  IF (i == nfreq) i = 999999
  WRITE(*, '(i6, " - ", i6, i7, "    ", f9.2)') i1, i, count, total
  IF (i == 999999) EXIT
  i = i + 1
END DO

STOP
END PROGRAM t_random_Zipf