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