subroutine aptbfac (nbase, idiga, ndiga, idigf, idigg, idigbw,
& idiggw, ndigm, idigb, ndigb, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBFAC
c
c call aptbfac (nbase, idiga, ndiga, idigf, idigg, idigbw,
c idiggw, ndigm, idigb, ndigb, nerr)
c
c Version: aptbfac Updated 2006 May 19 14:20.
c aptbfac Originated 2005 June 7 15:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the factorial function of the big integer "a", stored as
c the array of non-negative base nbase digits idiga, of length
c ndiga, and store the result as the array of non-negative base
c nbase digits idigb, of length ndigb. Integer arrays idigf
c and idigg, of length ndigm, are needed for intermediate
c results.
c
c The digits of idiga and idigb are the base nbase
c digits representing the non-negative decimal integers ideca
c and idecb, in order from most significant to least
c significant, using the equations:
c ideca = sum (idiga(n) * nbase^(N-n), n = 1, N = ndiga),
c idecb = sum (idigb(n) * nbase^(N-n), n = 1, N = ndigb),
c idecb = ideca!.
c
c Working arrays idigf, idigg, idiggw and idigbw are needed for
c intermediate results, and must each have its own memory space,
c of length ndigm or more.
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, idigf, idigg, idigbw, idiggw, ndigm.
c
c Output: idigb, ndigb, nerr.
c
c Calls: aptbmul, aptbsub
c
c Glossary:
c
c idiga Input A big integer "a", stored as an array of ndiga
c non-negative base nbase digits, each representing a
c single "digit" in the base nbase representation of
c "a", from most to least significant. If nbase > 10,
c each "digit" may require 2 or more integers.
c For example, for decimal integer 4821 and
c nbase = 16 (hexadecimal), idiga(n) = (5, 13, 2, 1),
c with ndiga = 4, or
c 4821 (dec) = 5 * 1 + 13 * 16 + 2 * 256 + 1 * 4096
c Big integer "a" must not be negative.
c
c idigb Output See idiga. Big integer "b" is the factorial function
c of "a".
c
c idigf-h Input See idiga. Arrays idigf, idigg, idiggw and idigbw are
c needed for temporary results, and must each have its
c own memory space, at least ndigm.
c
c nbase Input The number base for which the digits of integer arrays
c idiga, idigf, idigg and idigb are the
c coefficients of the powers of nbase, in order from
c most to least significant.
c
c ndiga Input The length of the integer array idiga.
c
c ndigb Output The length of the integer array idigb.
c Must be at big enough for the factorial of "a".
c
c ndigm Input The maximum number of words allowed in integer arrays
c idiga, idigf, idigg, idigbw and idigb. Memory space
c for idigf must be at least ndiga, and for idigg
c and idigb at least ndigm, and ndigm must be big
c enough to contain all of the digits of the factorial
c 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 negative.
c 3 if any digits idiga are negative.
c 4 if ndigm is too small.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
integer idiga(1) ! Big integer "a".
integer idigb(1) ! Factorial of "a".
integer idigbw(1) ! Temporary "b", the factorial of "a".
integer idigf(1) ! Current factor.
integer idigg(1) ! Current result, "f" - 1.
integer idiggw(1) ! Working space for aptbmul.
c.... Local variables.
integer n ! Index in big number array.
cbugc***DEBUG begins.
cbug 9900 format (/)
cbug 9901 format (/ 'aptbfac finding factorial 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 ('nbug=',i3,' a n =',i7,', idiga =',i7,'.')
cbug write ( 3, 9901) nbase, ndiga, ndigm
cbug nbug = 0
cbug if (ndiga .ne. 0) then
cbug write ( 3, 9900)
cbug write ( 3, 9902) (nbug, n, idiga(n), n = 1, ndiga)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
c.... Test for input errors.
nerr= 0
if (nbase .lt. 2) then
nerr = 1
go to 210
endif
if (ndiga .lt. 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
c.... Initialize.
ndigbw = 1
ndigf = 1
ndiggw = 1
ndigg = 1
do n = 1, ndigm
idigbw(n) = 0
idigf(n) = 0
idiggw(n) = 0
idigg(n) = 0
enddo
cbugcbugc***DEBUG begins.
cbugcbug 9974 format ('nbug=',i3,' bw n =',i7,', idigbw =',i7,'.')
cbugcbug 9976 format ('nbug=',i3,' f n =',i7,', idigf =',i7,'.')
cbugcbug 9977 format ('nbug=',i3,' gw n =',i7,', idiggw =',i7,'.')
cbugcbug 9978 format ('nbug=',i3,' g n =',i7,', idigg =',i7,'.')
cbugcbug nbug = 1
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Find "b", the factorial of "a".
c.... Test for "a" = 0. The factorial of zero is 1.
if ((ndiga .eq. 0) .or.
& ((ndiga .eq. 1) .and. (idiga(1) .eq. 0))) then
ndigbw = 1
idigbw(1) = 1
cbugcbugc***DEBUG begins.
cbugcbug nbug = 2
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbugc***DEBUG ends.
go to 210
endif
cbugcbugc***DEBUG begins.
cbugcbug nbug = 3
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Put "a" (reversed) into "bw" and "f".
ndigbw = ndiga
ndigf = ndiga
do n = 1, ndiga
na = ndiga + 1 - n
idigbw(n) = idiga(na)
idigf(n) = idiga(na)
enddo
cbugcbugc***DEBUG begins.
cbugcbug nbug = 4
cbugcbug if (ndiga .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9902) (nbug, n, idiga(n), n = 1, ndiga)
cbugcbug endif
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Return here for each factor.
nfact = 1
110 nfact = nfact + 1
cbugcbugc***DEBUG begins.
cbugcbug 9777 format (/ 'nfact=',i5)
cbugcbug write ( 3, 9777) nfact
cbugcbugc***DEBUG ends.
c.... Change "f" to most significant digit first.
call aptbrev (idigf, ndigf, nerrb)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 5
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Subtract 1 from factor "f" to get factor "g".
call aptbsub (nbase, idigf, ndigf, 1, 1, idiggw, ndigm,
& idigg, ndigg, isign, nerra)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 6
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Change "f", "g" to least significant digit first.
call aptbrev (idigf, ndigf, nerrb)
call aptbrev (idigg, ndigg, nerrb)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 7
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Stop when "g" is 1.
if ((ndigg .eq. 1) .and. (idigg(1) .eq. 0)) then
go to 210
endif
c.... Put new factor "g" into old factor "f".
ndigf = ndigg
do n = 1, ndigg
idigf(n) = idigg(n)
enddo
cbugcbugc***DEBUG begins.
cbugcbug nbug = 8
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Change "bw" and "f" to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigf, ndigf, nerrb)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 9
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Find new product "g" = "bw" * "f".
call aptbmul (nbase, idigbw, ndigbw, idigf, ndigf, idiggw, ndigm,
& idigg, ndigg, nerra)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 10
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Change "bw", "f" and "g" (if OK) to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigf, ndigf, nerrb)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 11
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Test for too many digits in "g".
if (nerra .eq. 0) then
call aptbrev (idigg, ndigg, nerrb)
cbugcbugc***DEBUG begins.
cbugcbug nbug = 12
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
else
nerr = 4
go to 210
endif
cbugcbugc***DEBUG begins.
cbugcbug if (nerra .ne. 0) then
cbugcbug write ( 3, 9911) nerra, ndigbw
cbugcbug endif
cbugcbug nbug = 13
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigf .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9976) (nbug, n, idigf(n), n = 1, ndigf)
cbugcbug endif
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Put product "g" into "bw".
ndigbw = ndigg
do n = 1, ndigg
idigbw(n) = idigg(n)
enddo
cbugcbugc***DEBUG begins.
cbugcbug nbug = 14
cbugcbug if (ndigg .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9978) (nbug, n, idigg(n), n = 1, ndigg)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Repeat until done.
go to 110
210 continue
c.... Put final product "bw" into factorial "b".
if (nerr .eq. 0) then
ndigb = ndigbw
do n = 1, ndigb
nb = ndigb + 1 - n
idigb(nb) = idigbw(n)
enddo
endif
cbugcbugc***DEBUG begins.
cbugcbug nbug = 15
cbugcbug if (ndigbw .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9974) (nbug, n, idigbw(n), n = 1, ndigbw)
cbugcbug endif
cbugcbug if (ndigb .ne. 0) then
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9913) (nbug, n, idigb(n), n = 1, ndigb)
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Zero out working arrays.
ndigbw = 1
ndigf = 1
ndigg = 1
ndiggw = 1
do n = 1, ndigm
idigbw(n) = 0
idigf(n) = 0
idigg(n) = 0
idiggw(n) = 0
enddo
if (nerr .ne. 0) then
ndigb = 0
do n = 1, ndigm
idigb(n) = -999999
enddo
endif
cbugc***DEBUG begins.
cbug 9911 format (/ 'aptbfac results: nerr=',i2,', ndigb =',i7,'.')
cbug 9913 format (' a! n =',i7,', idigb =',i7,'.')
cbug write ( 3, 9911) nerr, ndigb
cbug if (ndigb .ne. 0) then
cbug write ( 3, 9900)
cbug write ( 3, 9913) (nbug, n, idigb(n), n = 1, ndigb)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptbfac. (+1 line.)
end
UCRL-WEB-209832