subroutine aptcomg (n, lg, ng, tol, ncomb, comb, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTCOMG
c
c     call aptcomg (n, lg, ng, tol, ncomb, comb, nerr)
c
c     Version:  aptcomg  Updated    2007 February 12 14:00.
c               aptcomg  Originated 2007 February 2 13:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c     Source:   In directory ~edwards/work/apt/src.
c
c     Purpose:  to find C(n, lg(1), lg(2), ..., lg(ng)), the number of ways to
c               select ng groups of items of lengths lg(1), lg(2), ..., lg(ng),
c               from a total of n items, without considering the order of items
c               in each group, or the order of the groups, and return the
c               integer and floating point values of C, ncomb and comb,
c               respectively.
c               The sum lgtot of the group lengths lg must not exceed n.
c               If lgtot is less than n, the result is the same as if one
c               additional group length lg(ng+1) = n - lgtot is added.
c               
c
c     Note:     The formula is C = ncomb = comb = n! / (PL * PR), where
c               PL = lg(1)! * lg(2)! * ... * lg(ng)! * (n - lgtot)!, and
c               PR = nrep(1)! * nrep(2)! * ... * nrep(n)!, and
c               nrep(i), (i = 1, n) is the number of groups with length i.
c               
c
c     Input:    n, lg, ng, tol.
c
c     Output:   ncomb, comb, nerr.
c
c     Glossary:
c
c     comb      Output   The floating point value of C.
c                          If n < 0, returned as -1.e99, with nerr = 1.
c                          WARNING:  may overflow largest machine floating point
c                          number, if n is too large.
c
c     n         Input    Total number of items.
c
c     nerr      Output   Indicates an input or result error, if not zero.
c                          -1 if ncomb exceeds the largest machine integer.
c                           1 if n is not positive.
c                           2 if ng is not positive.
c                           3 if any lg is not positive.
c                           4 if lgtot > n.
c
c     ncomb     Output   The integer value of C.
c                          If n < 0, returned as -999999, with nerr = 1.
c                          If larger then the largest machine integer,
c                          returned as -999999, with nerr = -1, but comb will
c                          still be accurate, unless comb overflows the largest
c                          machine floating point number.
c
c     lg        Input    The lengths of the ng groups of items to be selected.
c
c     ng        Input    The number of groups of items to be selected.
c
c     tol       Input    The precision with which ncomb and comb must agree,
c                          to insure that ncomb does not exceed the largest
c                          machine integer.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

      implicit none

c.... Arguments.

      integer  lg(1)                  ! Items in each group selected from n.
      integer  n                      ! Total number of items.
      integer  nerr                   ! Error indicator.
      integer  ncomb                  ! Integer value of C.
      integer  ng                     ! Total number of groups.

      real     comb                   ! Real value of C.
      real     tol                    ! Reqired precision of comb, ncomb match.

c.... Local variables.

      integer  i                      ! Index in lg.
      integer  lgtot                  ! Sum of lg values.
      integer  nfact                  ! Integer factorial of an integer.
      integer  nn                     ! Possible group length.
      integer  nrem                   ! Difference n - lgtot. 
      integer  nrep                   ! Number of times a group length repeated.

      real     fact                   ! Floating point factorial of an integer.
      real     fdiff                  ! Difference between floating pt values.
      real     fncomb                 ! Real value of ncomb.

cbugc***DEBUG begins.
cbug 9001 format (/ 'aptcomg finding ways to separate ',i5,' items',
cbug     &  ' into',i5,'  groups of sizes' / 16i5)
cbug 9002 format ('  tol =',1pe22.14)
cbug      write ( 3, 9001) n, ng, (lg(i), i = 1, ng)
cbug      write ( 3, 9002) tol
cbugc***DEBUG ends.

c.... Initialize.

      nerr   = 0
      ncomb  = -999999
      fncomb = -999999.0
      comb   = -1.e99

c.... Test for errors, limits.

      if (n .le. 0) then
        nerr = 1
        go to 210
      endif

      if (ng .le. 0) then
        nerr = 2
        go to 210
      endif

      lgtot = 0
      do i = 1, ng
        if (lg(i) .le. 0) then
          nerr = 3
          go to 210
        endif
        lgtot = lgtot + lg(i)
      enddo

      if (lgtot .gt. n) then
        nerr = 4
        go to 210
      endif

c.... Find the number of ways to select ng groups of lengths lg(i), i = 1, ng,
c....   from n items.

c.... Find the numerator, the factorial of n.

      call aptfact (n, tol, ncomb, comb, nerr)

c.... Find and divide by the factorial of each group length.

      do i = 1, ng

        call aptfact (lg(i), tol, nfact, fact, nerr)

        ncomb = ncomb / nfact
        comb  = comb  / fact

      enddo

c.... Find and divide by the factorial of n - lgtot.

      nrem = n - lgtot

      if (nrem .gt. 1) then

        call aptfact (nrem, tol, nfact, fact, nerr)
  
        ncomb = ncomb / nfact
        comb  = comb  / fact

      endif

c.... Find and divide by the number of repeats of each group length.

      do nn = 1, n

        nrep = 0
        do i = 1, ng
          if (lg(i) .eq. nn) then
            nrep = nrep + 1
          endif
        enddo

        if (nrem .eq. nn) then
          nrep = nrep + 1
        endif

        if (nrep .gt. 0) then

          call aptfact (nrep, tol, nfact, fact, nerr)

          ncomb = ncomb / nfact
          comb  = comb  / fact
        endif

      enddo

c.... See if the integer value is correct.

      nerr   = 0
      ncomb  = comb + 0.5
      fncomb = ncomb
      fdiff  = fncomb - comb
      if (abs (fdiff) .gt. tol * comb) then
cbugc***DEBUG begins.
cbug 7770 format ('aptcomb DEBUG.' /
cbug     &  '   ncomb =',i22 /
cbug     &  '  fncomb =',1pe22.14 /
cbug     &  '    comb =',1pe22.14 )
cbug        write ( 3, 7770) ncomb, fncomb, comb
cbugc***DEBUG ends.
        nerr  = -1
        ncomb = -999999
      endif

  210 continue
cbugc***DEBUG begins.
cbug 9003 format (/ 'aptcomg found # of ways.  nerr = ',i3 /
cbug     &  'nc, fnc, chk =',i20,2x,1p2e22.14 /
cbug     &  'fdiff =',1pe22.14)
cbug      write ( 3, 9003) nerr, ncomb, comb, fncomb, fdiff
cbugc***DEBUG ends.

      return

      end

UCRL-WEB-209832