c solves problem 6.10 of Tibshirani and Efron by both methods
c fast method: enumerate only distinct samples
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL X(8), Y(8), B(8)
      INTEGER JP(8)
      LOGICAL FINISH
      DATA X /1.2, 3.5, 4.7, 7.3, 8.6, 12.4, 13.8, 18.1/
      CALL TIMER(ISTART)
      NB = 8
      X1 = 0.0D0
      X2 = 0.0D0
      DO K = 1, 8
       JP(K) = 0
      ENDDO
      FINISH = .FALSE.
      DO WHILE (.NOT.FINISH)
       CALL GENERATE_PARTITION(JP,NB,FINISH)
       KI = 0
       DO K = 1, 8
        IF (JP(K).GT.0) THEN
         DO I = 1, JP(K)
          KI = KI + 1
          B(KI) = X(K)
         ENDDO
        ENDIF
       ENDDO
       CALL QCKSRT(NB,B)
       TMEAN = 0.0D0
       DO K = 3, 6
        TMEAN = TMEAN + B(K)
       ENDDO
       TMEAN = TMEAN / 4.0D0
       PROBJ = PROBABILITY610(JP)
       X1 = X1 + PROBJ * TMEAN
       X2 = X2 + PROBJ * TMEAN ** 2
      ENDDO
      SE = DSQRT(X2-X1**2)
      CALL TIMER(IEND)
      SECS = FLOAT(IEND-ISTART) / 100.0D0
      WRITE (6,'(''   FAST METHOD'')')
      WRITE (6,'(''   ideal  se          = '',  F15.6)') SE
      WRITE (6,'(''   ideal mean         = '',  F15.6)') X1
      WRITE (6,'(''   time (in seconds)  = '',  F15.6)') SECS
      STOP
      END
c
      SUBROUTINE GENERATE_PARTITION(J,N,FINISH)
c output variables:
c generate all non-negative integer solutions to the equation
c  j_1+j_2+...+j_n = n
c Algorithm:
c Step 1. Start with j == (n,0,...,0) = (j_1, ..., j_n)
c Step 2. Find the largest index, i, such that j_i + j_(i+1) + ... + j_n < n.
c Step 3. Increment j_i and then set everything to its left to 0 except for
c        j_1 which is set to the largest possible value.
c Step 4. Repeat Steps 2 & 3 until (0,...,n) is reached. This occurs when
c        there is no value of j_i from Step 2 which can be found.
c

      LOGICAL FINISH
      INTEGER J(N)
      FINISH = .TRUE.
      NSUM = 0
      DO I = 1, N
       NSUM = NSUM + J(I)
      ENDDO
      IF (NSUM.EQ.0) THEN
       J(1) = N
      ELSE
       I = N
       JSUM = 0
       DO WHILE (JSUM.LT.N)
        JSUMPREV = JSUM
        JSUM = JSUM + J(I)
        I = I - 1
       ENDDO
       I = I + 2
       IF (I.EQ.N+1) RETURN
       J(I) = J(I) + 1
       IF (I.GT.1) THEN
        DO K = 2, I - 1
         J(K) = 0
        ENDDO
       ENDIF
       J(1) = N - JSUMPREV - 1
      ENDIF
      FINISH = .FALSE.
      RETURN
      END
c
      DOUBLE PRECISION FUNCTION PROBABILITY610(J)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER J(8)
      X = DLOGAM(DFLOAT(9))
      DO I = 1, 8
       X = X - DLOGAM(DFLOAT(J(I)+1))
      ENDDO
      PROBABILITY610 = DEXP(X-8.0D0*DLOG(8.0D0))
      RETURN
      END
c version: Feb 23 1992
      SUBROUTINE QCKSRT(N, ARR)
      PARAMETER (M = 7, NSTACK = 50, FM = 7875., FA = 211., FC = 1663., 
     1 FMI = 1.2698413E-4)
      DIMENSION ARR(N), ISTACK(NSTACK)
      JSTACK = 0
      L = 1
      IR = N
      FX = 0.
      DO WHILE( .TRUE. )
        IF( IR - L .LT. M ) THEN
          DO J = L + 1, IR
            A = ARR(J)
            DO I = J - 1, 1, -1
              IF( ARR(I) .LE. A ) GOTO 12
              ARR(I + 1) = ARR(I)
            ENDDO
            I = 0
12          ARR(I + 1) = A
          ENDDO
          IF( JSTACK .EQ. 0 ) RETURN
          IR = ISTACK(JSTACK)
          L = ISTACK(JSTACK - 1)
          JSTACK = JSTACK - 2
        ELSE
          I = L
          J = IR
          FX = MOD(FX*FA + FC, FM)
          IQ = L + (IR - L + 1)*(FX*FMI)
          A = ARR(IQ)
          ARR(IQ) = ARR(L)
          DO WHILE( .TRUE. )
            DO WHILE( J .GT. 0 )
              IF( A .GE. ARR(J) ) EXIT
              J = J - 1
            ENDDO
            IF( J .LE. I ) GOTO 31
            ARR(I) = ARR(J)
            I = I + 1
            DO WHILE( I .LE. N )
              IF( A .LE. ARR(I) ) EXIT
              I = I + 1
            ENDDO
            IF( J .LE. I ) EXIT
            ARR(J) = ARR(I)
            J = J - 1
          ENDDO
          ARR(J) = A
          I = J
          GOTO 30
31        ARR(I) = A
30        JSTACK = JSTACK + 2
          IF( JSTACK .GT. NSTACK ) PAUSE 'NSTACK must be made larger.'
          IF( IR - I .GE. I - L ) THEN
            ISTACK(JSTACK) = IR
            ISTACK(JSTACK - 1) = I + 1
            IR = I - 1
          ELSE
            ISTACK(JSTACK) = I - 1
            ISTACK(JSTACK - 1) = L
            L = I + 1
          ENDIF
        ENDIF
      ENDDO
      END
c DLOGAM.                                                                       
c LOG OF GAMMA FUNCTION.                                                        
c                                                                               
      DOUBLE PRECISION FUNCTION DLOGAM(X)
c                                                                               
c    THIS FUNCTION EVALUATES THE NATURAL LOGARITHM OF GAMMA(X)                  
c    FOR ALL X>0, ACCURATE TO 10 DECIMAL PLACES. STIRLING'S                     
c    FORMULA IS USED FOR THE CENTRAL POLYNOMIAL PART.                           
c                                                                               
c REFERENCE: PIKE AND HILL, CACM ALGORITHM 291,                                 
c     COMM. ACM 9 ,SEPT 1966, 684                                               
c                                                                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA ZERO, HALF, ONE, SEVEN /0.0D0, 0.5D0, 1.0D0, 7.0D0/
      DATA A1 /0.918938533204673D0/, A2 / -0.000595238095238D0/, A3 /
     +  0.000793650793651D0/, A4 / -0.00277777777777778D0/, A5 /
     +  0.0833333333333333D0/
      XX = X
      F = ZERO
      DLOGAM = ZERO
      IF (X.LE.ZERO) RETURN
      IF (X.LT.SEVEN) THEN
       F = ONE
       Z = X - ONE
   10  Z = Z + ONE
       IF (Z.LT.SEVEN) THEN
        XX = Z
        F = F * Z
        GOTO 10
       ENDIF
       XX = XX + ONE
       F = -DLOG(F)
      ENDIF
      Z = ONE / (XX*XX)
      DLOGAM = F + (XX-HALF) * DLOG(XX) - XX + A1 + (((A2*Z+A3)*Z+A4)*Z+
     +  A5) / XX
      RETURN
      END