*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
PROGRAM entropy ! compute the Brillouin and DeMoivre-Shannon indices
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
IMPLICIT NONE
INTEGER s ; PARAMETER (s = 144 000)
INTEGER N(s), i, j, sact
REAL f, f2, f3, brill, brill2, shannon, RGAUSS
DOUBLE PRECISION DGAMLN, d
! DATA N / 24, 20, 7, 3, 3, 1 /
! sact = 6
CALL SRAND(TIME())
DO j = 1, 100
* choose random size array
sact = INT(RAND(0) * s) + 1
* get a random input array
DO i = 1, sact
f = RAND(0) * 180.
N(i) = ABS(INT(RGAUSS(f, f / 2)))
END DO ! next i
! CYCLE
f = brill(N,sact)
f2 = brill2(N,sact)
f3 = shannon(N,sact)
PRINT *, 'result: ', f, ' vs ', f2, ' vs ', f3
END DO
END
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
* Brillouin's index by Sterling's approximation
* H = (1 / N) @ LOG(N! / PRODUCT(n[i]!))
* and n! = (2 @ pi @ n)^(1/2) @ (n/e)^n
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
FUNCTION brill(N, s)
IMPLICIT NONE
INTEGER N, s
REAL brill
DIMENSION N(s)
* local variables
INTEGER i, strim
REAL sum, snphln, ln2pi, f
PARAMETER (ln2pi = 1.8378771)
* begin
brill = -999.
strim = s
* find SUM n[i] and SUM (n[i]+1/2) @ LOG(n[i])
sum = 0 ; snphln = 0
DO i = 1, s
f = FLOAT(N(i))
IF (f < 0.5) THEN ! skip empty groups
strim = strim - 1
CYCLE
END IF
sum = sum + f
snphln = snphln + (f + 0.5) * LOG(f)
END DO
IF (strim .EQ. 0) RETURN ! return negative if no data
brill = (1 / sum) * ( (sum+0.5) * LOG(sum)
& - snphln - (FLOAT(strim-1) / 2) * ln2pi )
END ! of brill
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
* Brillouin's Index by the log-gamma function
* H = (1 / N) @ LOG(N! / PRODUCT(n[i]!))
* gamma(n+1) = n!
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
FUNCTION brill2(N, s)
IMPLICIT NONE
INTEGER N, s
REAL brill2
DIMENSION N(s)
* local variables
INTEGER i
DOUBLE PRECISION d, DGAMLN, sum, dif
* find SUM n[i]
sum = 0 ; dif = 0
DO i = 1, s
d = DBLE(N(i))
sum = sum + d
dif = dif + DGAMLN(d + 1)
END DO
brill2 = ( DGAMLN(sum + 1) - dif ) / sum
END ! of brill2
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
* the usual DeMoivre-Shannon entropy
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
FUNCTION shannon(N, s)
IMPLICIT NONE
INTEGER N, s
REAL shannon
DIMENSION N(s)
* local variables
INTEGER i
REAL sum, snln, f
* begin
sum = 0 ; snln = 0
DO i = 1, s
IF (N(i) .EQ. 0) CYCLE
f = FLOAT(N(i))
sum = sum + f
snln = snln + f * LOG(f)
END DO
shannon = LOG(sum) - snln / sum
END ! of shannon
C***BEGIN PROLOGUE RGAUSS
C***PURPOSE Generate a normally distributed (Gaussian) random number.
C***LIBRARY SLATEC (FNLIB)
C***CATEGORY L6A14
C***TYPE SINGLE PRECISION (RGAUSS-S)
C***KEYWORDS FNLIB, GAUSSIAN, NORMAL, RANDOM NUMBER, SPECIAL FUNCTIONS
C***AUTHOR Fullerton, W., (LANL)
C***DESCRIPTION
C
C Generate a normally distributed random number, i.e., generate random
C numbers with a Gaussian distribution. These random numbers are not
C exceptionally good -- especially in the tails of the distribution,
C but this implementation is simple and suitable for most applications.
C See R. W. Hamming, Numerical Methods for Scientists and Engineers,
C McGraw-Hill, 1962, pages 34 and 389.
C
C Input Arguments --
C XMEAN the mean of the Guassian distribution.
C SD the standard deviation of the Guassian function
C EXP (-1/2 * (X-XMEAN)**2 / SD**2)
C
C***REFERENCES (NONE)
C***ROUTINES CALLED RAND
C***REVISION HISTORY (YYMMDD)
C 770401 DATE WRITTEN
C 861211 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 910819 Added EXTERNAL statement for RAND due to problem on IBM
C RS 6000. (WRB)
C***END PROLOGUE RGAUSS
FUNCTION RGAUSS (XMEAN, SD)
C***FIRST EXECUTABLE STATEMENT RGAUSS
RGAUSS = -6.0
DO 10 I=1,12
RGAUSS = RGAUSS + RAND(0)
10 CONTINUE
C
RGAUSS = XMEAN + SD*RGAUSS
C
RETURN
END