subroutine aptbpac (nbase, idiga, ndiga, idigb, ndigb,
& idigaf, idigbf, idigr, idigv, idigvf, idigw,
& idigx, idigy, idigz, idigcw, idigpw, ndigm,
& idigp, ndigp, idigc, ndigc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBPAC
c
c call aptbpac (nbase, idiga, ndiga, idigb, ndigb,
c & idigaf, idigbf, idigr, idigv, idigvf, idigw,
c & idigx, idigy, idigz, idigcw, idigpw, ndigm,
c & idigp, ndigp, idigc, ndigc, nerr)
c
c Version: aptbpac Updated 2006 May 15 18:20.
c aptbpac Originated 2005 May 15 18:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the number of permutations "p" and combinations "c" of
c a set of "a" things taken "b" at a time, where "a", "b", "c"
c and "p" are big integers with lengths ndiga, ndigb, ndigc and
c ndigp, respectively, stored as one base nbase digit per machine
c word, from most significant to least significant.
c
c Integer arrays idigaf, idigbf, idigv, idigvf,
c idigw, idigx, idigy and idigz, of length ndigm or more, are
c needed as working arrays for intermediate results.
c
c Flag nerr indicates any input error.
c
c See aptbtod, aptdtob, aptbadd, aptbsub, aptbmul, aptbdiv,
c aptbrev, aptbrtn, aptbfac.
c
c Input: nbase, idiga, ndiga, idigb and ndigb (data), idigaf, idigbf,
c idigr, idigv, idigvx, idigw, idigx, idigy and idigz
c (working space), and ndigm (max array size).
c
c Output: idigc, ndigw, idigp, ndigp, nerr.
c
c Calls: aptbdiv, aptbfac, aptbtod, aptbmul
c
c Glossary:
c
c idiga Input The number of items "a" from which "b" items may be
c selected. Big integer "a" must be positive, and is
c stored as an array of ndiga non-negative base nbase
c digits, each representing a single "digit" in the
c representation of "a", from most to least
c significant. If nbase > 10, each "digit" may require
c 2 or more integers. For example, for decimal integer
c 4821 and nbase = 16 (hexadecimal),
c idiga(n) = (5, 13, 2, 1), with ndiga = 4, or
c 4821 (dec) = 5 * 1 + 13 * 16 + 2 * 256 + 1 * 4096
c
c idigb Input See idiga. Big integer "b" is the number of items
c which may be selected from "a" items, and must be
c in the range from 1 to "a".
c
c idigc Output See idiga. Big integer "c" is the number of
c combinations "c" of "a" things taken "b" at a time,
c and is equal to the factorial of "a" divided by the
c product of the factorials of "b" and "a" - "b".
c
c idigaf Input Temporary array for the factorial of "a".
c Size ndigm or more.
c
c idigbf Input Temporary array for the factorial of "b".
c Size ndigm or more.
c
c idigr Input Temporary array for the remainders of divisions.
c Size ndigm or more.
c
c idigv Input Temporary working array for "a" - "b".
c Size ndigm or more.
c
c idigvf Input Temporary working array for factorial of "a" - "b".
c Size ndigm or more.
c
c idigw-z Input Temporary working arrays.
c Size ndigm or more.
c
c nbase Input The number base for which the digits of integer arrays
c are the coefficients of the powers of nbase,
c in order from most to least significant.
c
c ndiga Input The length of the integer array "a".
c Must be positive.
c
c ndigb Input The length of the integer array "b".
c Must be positive.
c
c ndigc Output The length of the integer array "c".
c Must be large enough to hold the digits of the
c factorial of "a". For example, if "a" = 450,
c ndigc must be at least 1001.
c
c ndigm Input The maximum number of words allowed in integer arrays
c idiga, idigb, idigcw, idigpw, idigw, idigx, idigy and
c idigz. Memory space for idigw, idigx, idigy, idigz,
c idigcw and idigpw must be big enough to hold all the
c digits of the factorial of "a".
c
c ndigp Output The length of the integer array "p".
c Must be large enough to hold the digits of the
c factorial of "a".
c
c nerr Output Indicates an input or a result error, if not zero.
c 1 if nbase is less than 2.
c 2 if ndiga is not positive
c 3 if any digits idiga in "a" are negative.
c 4 if any digits idigb in "b" are negative.
c 5 if "b" is bigger than "a".
c 6 if ndigm is too small.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
integer idiga(1) ! Big integer "a".
integer idigb(1) ! Big integer "b".
integer idigc(1) ! Big integer "c", number of combinations.
integer idigcw(1) ! Big integer "c", number of combinations.
integer idigaf(1) ! Working array "af" = "a" factorial.
integer idigbf(1) ! Working array "bf" = "b" factorial.
integer idigp(1) ! Big integer "p", number of permutations.
integer idigpw(1) ! Big integer "p", number of permutations.
integer idigr(1) ! Working array "r", temporary remainder.
integer idigv(1) ! Working array "v" = "a" - "b".
integer idigvf(1) ! Working array "vf" = "a"-"b" factorial.
integer idigw(1) ! Working array.
integer idigx(1) ! Working array.
integer idigy(1) ! Working array.
integer idigz(1) ! Working array.
c.... Local variables.
integer n ! Index in a big number array.
cbugc***DEBUG begins.
cbug 9900 format (/)
cbug 9901 format (/ 'aptbpac finding permutations and combinations',
cbug & ' of the base nbase =',i7 /
cbug & ' digit array idiga of length',i7 /
cbug & ' to get base nbase digit array idigb.' /
cbug & ' ndigm =',i7,'.')
cbug 9902 format (' a n =',i7,', idiga =',i7,'.')
cbug 9903 format (' b n =',i7,', idigb =',i7,'.')
cbug write ( 3, 9901) nbase, ndiga, ndigm
cbug write ( 3, 9900)
cbug write ( 3, 9902) (n, idiga(n), n = 1, ndiga)
cbug write ( 3, 9900)
cbug write ( 3, 9903) (n, idigb(n), n = 1, ndigb)
cbugc***DEBUG ends.
c.... Test for input errors.
nerr= 0
if (nbase .lt. 2) then
nerr = 1
go to 210
endif
if ((ndiga .lt. 1) .or.
& ((ndiga .eq. 1) .and. (idiga(1) .eq. 0))) then
nerr = 2
go to 210
endif
if (ndiga .ge. 1) then
do n = 1, ndiga
if (idiga(n) .lt. 0) then
nerr = 3
go to 210
endif
enddo
endif
if (ndigb .ge. 1) then
do n = 1, ndigb
if (idigb(n) .lt. 0) then
nerr = 4
go to 210
endif
enddo
endif
c.... Find "v" = "a" - "b".
ndigv = 1
do n = 1, ndigm
idigv(n) = 0
enddo
call aptbsub (nbase, idiga, ndiga, idigb, ndigb,
& idigx, ndigm, idigv, ndigv, isign, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9979 format ('nbug=',i3,' "v" n =',i7,', idigv =',i7,'.')
cbugcbug 9980 format ( / )
cbugcbug nbug = 1
cbugcbug if (ndigv .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9979) (nbug, n, idigv(n), n = 1, ndigv)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... See if "b" exceeds "a".
if (nerra .ne. 0) then
nerr = 5
go to 210
endif
c.... See if "b" = "a".
if ((ndigv .le. 0) .or.
& ((ndigv .eq. 1) .and. (idigv(1) .eq. 0))) then
ndigcw = 1
idigcw = 1
ndigpw = 1
idigpw = 1
go to 210
endif
c.... Initialize.
ndigaf = 1
ndigbf = 1
ndigpw = 1
ndigcw = 1
ndigr = 1
ndigvf = 1
ndigw = 1
ndigx = 1
ndigy = 1
ndigz = 1
do n = 1, ndigm
idigaf(n) = 0
idigbf(n) = 0
idigpw(n) = 0
idigcw(n) = 0
idigr(n) = 0
idigvf(n) = 0
idigw(n) = 0
idigx(n) = 0
idigy(n) = 0
idigz(n) = 0
enddo
c.... Find "af", the factorial of "a".
call aptbfac (nbase, idiga, ndiga, idigw, idigx, idigy,
& idigz, ndigm, idigaf, ndigaf, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9982 format ('nbug=',i3,' "af" n =',i7,', idigaf =',i7,'.')
cbugcbug nbug = 2
cbugcbug if (ndigaf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9982) (nbug, n, idigaf(n), n = 1, ndigaf)
cbugcbug endif
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 6
go to 210
endif
c.... Find "bf", the factorial of "b".
call aptbfac (nbase, idigb, ndigb, idigw, idigx, idigy,
& idigz, ndigm, idigbf, ndigbf, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9983 format ('nbug=',i3,' "bf" n =',i7,', idigbf =',i7,'.')
cbugcbug nbug = 3
cbugcbug if (ndigbf .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9983) (nbug, n, idigbf(n), n = 1, ndigbf)
cbugcbug endif
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 6
go to 210
endif
c.... Find "vf", the factorial of "v".
call aptbfac (nbase, idigv, ndigv, idigw, idigx, idigy,
& idigz, ndigm, idigvf, ndigvf, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9984 format ('nbug=',i3,' "vf" n =',i7,', idigvf =',i7,'.')
cbugcbug nbug = 4
cbugcbug if (ndigvf .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9984) (nbug, n, idigvf(n), n = 1, ndigvf)
cbugcbug endif
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 6
go to 210
endif
c.... Find "pw", the number of permutations of "a" things taken "b" at a time.
call aptbdiv (nbase, idigaf, ndigaf, idigvf, ndigvf,
& idigv, idigw, idigx, idigy, idigz, ndigm,
& idigpw, ndigpw, idigr, ndigr, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9977 format ('nbug=',i3,' "r" n =',i7,', idigr =',i7,'.')
cbugcbug 9981 format ('nbug=',i3,' "p" n =',i7,', idigpw =',i7,'.')
cbugcbug nbug = 5
cbugcbug if (ndigr .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9977) (nbug, n, idigr(n), n = 1, ndigr)
cbugcbug endif
cbugcbug if (ndigpw .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9981) (nbug, n, idigpw(n), n = 1, ndigpw)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Find "cw", the number of combinations of "a" things taken "b" at a time.
call aptbdiv (nbase, idigpw, ndigpw, idigbf, ndigbf,
& idigv, idigw, idigx, idigy, idigz, ndigm,
& idigcw, ndigcw, idigr, ndigr, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9978 format ('nbug=',i3,' "c" n =',i7,', idigcw =',i7,'.')
cbugcbug nbug = 6
cbugcbug if (ndigr .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9977) (nbug, n, idigr(n), n = 1, ndigr)
cbugcbug endif
cbugcbug if (ndigcw .ne. 0) then
cbugcbug write ( 3, 9980)
cbugcbug write ( 3, 9978) (nbug, n, idigcw(n), n = 1, ndigcw)
cbugcbug endif
cbugcbugc***DEBUG ends.
210 continue
c.... Put "cw" into "c", "pw" into "p".
ndigc = ndigcw
do n = 1, ndigc
idigc(n) = idigcw(n)
enddo
ndigp = ndigpw
do n = 1, ndigp
idigp(n) = idigpw(n)
enddo
c.... Zero out working arrays.
ndigaf = 1
ndigbf = 1
ndigcw = 1
ndigpw = 1
ndigr = 1
ndigv = 1
ndigvf = 1
ndigw = 1
ndigx = 1
ndigy = 1
ndigz = 1
do n = 1, ndigm
idigaf(n) = 0
idigbf(n) = 0
idigcw(n) = 0
idigpw(n) = 0
idigr(n) = 0
idigv(n) = 0
idigvf(n) = 0
idigw(n) = 0
idigx(n) = 0
idigy(n) = 0
idigz(n) = 0
enddo
if (nerr .ne. 0) then
ndigcw = 0
ndigpw = 0
do n = 1, ndigm
idigcw(n) = -999999
idigpw(n) = -999999
enddo
endif
cbugc***DEBUG begins.
cbug 9911 format (/ 'aptbpac results: nerr=',i2,', ndigc =',i7,
cbug & ', ndigp =',i7,'.')
cbug 9913 format (' "c" n =',i7,', idigc =',i7,'.')
cbug 9914 format (' "p" n =',i7,', idigp =',i7,'.')
cbug write ( 3, 9911) nerr, ndigc, ndigp
cbug write ( 3, 9900)
cbug write ( 3, 9913) (n, idigc(n), n = 1, ndigc)
cbug write ( 3, 9900)
cbug write ( 3, 9914) (n, idigp(n), n = 1, ndigp)
cbugc***DEBUG ends.
return
c.... End of subroutine aptbpac. (+1 line.)
end
UCRL-WEB-209832