subroutine aptbdiv (nbase, idiga, ndiga, idigb, ndigb,
& idigqw, idigrw, idigbw, idigy, idigz, ndigm,
& idigq, ndigq, idigr, ndigr, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBDIV
c
c call aptbdiv (nbase, idiga, ndiga, idigb, ndigb,
c & idigqw, idigrw, idigbw, idigy, idigz, ndigm,
c & idigq, ndigq, idigr, ndigr, nerr)
c
c Version: aptbdiv Updated 2006 June 6 13:20.
c aptbdiv Originated 2005 May 23 14:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To divide "a" by "b" to find the quotient "q" and the remainder
c "r". Each of these are big integers, stored as the arrays of
c non-negative base nbase digits idiga, idigb, idigq and idigr,
c of lengths ndiga, ndigb, ndigq, and ndigr, respectively,
c stored in normal order (most significant digit first).
c
c Integer arrays idigqw, idigrw, idigbw, idigy and idigz are
c needed for intermediate results, and each must have its own
c memory space, of length ndigm or more.
c
c Note: to find nbase^N times the reciprocal of "b", enter
c "a" as nbase^N, where N is much larger than ndigb.
c
c The digits of idiga, idigb, idigq and idigr are the
c base nbase digits representing the non-negative decimal integers
c ideca, idecb, idecq and idecr, in order from most significant
c to least 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 idecq = sum (idigq(n) * nbase^(N-n), n = 1, N = ndigq)
c idecr = sum (idigr(n) * nbase^(N-n), n = 1, N = ndigr)
c
c Flag nerr indicates any input error.
c
c See aptbtod, aptdtob, aptbadd, aptbsub, aptbmul, aptbpow,
c aptbrtn, aptbfac.
c
c Input: nbase, idiga, ndiga, idigb, ndigb, idigqw, idigrw, idigbw,
c idigy, idigz, ndigm.
c
c Output: idigq, ndigq, idigr, ndigr, nerr.
c
c Calls: aptbsub
c
c Glossary:
c
c idiga Input Dividend "a", a big integer stored as idiga, an array
c of ndiga non-negative base nbase digits, each
c representing a single "digit" in the base nbase
c representation of "a". If nbase > 10, each "digit"
c may require 2 or more integers. For example, for
c nbase = 16, and "a" = decimal integer 4821,
c idiga(n) = (5, 13, 2, 1), ndiga = 4. This means that
c 4821 (dec) = 5 * 1 + 13 * 16 + 2 * 256 + 1 * 4096.
c
c idigb Input Divisor "b", a big integer stored as idigb. The memory
c space for idigb must be at least as large as ndiga.
c See idiga.
c
c idigbw Input Temporary array for divisor idigb. Initially reversed
c "b". Must not be idiga, idigb, idigq or idigr, but
c must have its own memory space, of length ndigm.
c
c idigq Output Quotient "q", a big integer stored as idigq. Array
c idigq is the integer part of the quotient resulting
c from dividing big integer "a" by big integer "b".
c The memory space for idigq must be at least as large
c as ndiga. See idiga.
c
c idigqw Input Temporary array for idigq. Must not be idiga, idigb,
c idigq or idigr, but must have its own memory space,
c of length ndigm.
c
c idigr Output Remainder "r", a big integer stored as idigr.
c Array idigr is the remainder resulting from
c dividing big integer "a" by big integer "b".
c Array idigq must not be idigq, but must have its own
c memory space, at least as large as ndiga. See idiga.
c
c idigrw Input Temporary array for idigr. Initially reversed "a".
c Must not be idiga, idigb, idigq or idigr,
c but must have its own memory space, of length
c ndigm.
c
c idigy-z Input Working arrays of integers, of length ndigm.
c Must not be idiga, idigb, idigq or idigr,
c but must have its own memory space, of length
c ndigm or more.
c
c nbase Input The number base for which the digits of integer arrays
c idiga, idigb, idigq and idigr 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 Input The length of the integer array idigb.
c
c ndigm Input The maximum number of words allowed in integer arrays
c idigqw, idigrw, idigbw, idigy, idigz. Memory space
c for each of these must be at least ndigm, and ndigm
c must be at least as big as ndiga.
c
c ndigq Output The length of the integer array idigq, <= ndiga.
c
c ndigr Output The length of the integer array idigr, <= ndiga.
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 of idiga are negative.
c 4 if big integer "b" is not positive.
c 5 if any digits of idigb are negative.
c 6 if ndigm is less than ndiga.
c 7 if ndigm is less than ndigb.
c 8 if "a" and "b" are both zero.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
integer idiga(1) ! Array "a". Numerator.
integer idigb(1) ! Array "b". Denominator.
integer idigbw(1) ! Array "bw". Temporary for idigb.
integer idigq(1) ! Array "q". Quotient.
integer idigqw(1) ! Array "qw". Temporary for "q".
integer idigr(1) ! Array "r". Remainder.
integer idigrw(1) ! Array "rw". Temporary for "r".
integer idigy(1) ! Array "y". Temporary for aptbsub use.
integer idigz(1) ! Array "z". Temporary for aptbsub use.
cbugcbugc***DEBUG begins.
cbugcbug character*64 abug ! Debug comment.
cbugcbugc***DEBUG ends.
c.... Local variables.
cbugc***DEBUG begins.
cbug 9900 format (/)
cbug 9901 format (/ 'aptbdiv dividing the base nbase digit array',
cbug & ' idiga by idigb' /
cbug & ' to get base nbase digit array idigq and',
cbug & ' remainder idigr.' /
cbug & ' nbase =',i7,', ndiga =',i7,', ndigb =',i7,'.' /
cbug & ' ndigm =',i7 )
cbug write ( 3, 9900)
cbug write ( 3, 9901) nbase, ndiga, ndigb, ndigm
cbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug 9100 format (/ 'aptbdiv DEBUG. ',a64)
cbugcbug 9301 format (' "a" n ',i7,', idiga =',i7,'.')
cbugcbug 9302 format (' "b" n ',i7,', idigb =',i7,'.')
cbugcbug 9303 format (' "bw" n ',i7,', idigbw =',i7,'.')
cbugcbug 9304 format (' "q" n ',i7,', idigq =',i7,'.')
cbugcbug 9305 format (' "qw" n ',i7,', idigqw =',i7,'.')
cbugcbug 9306 format (' "r" n ',i7,', idigr =',i7,'.')
cbugcbug 9307 format (' "rw" n ',i7,', idigrw =',i7,'.')
cbugcbug 9308 format (' "y" n ',i7,', idigy =',i7,'.')
cbugcbug 9309 format (' "z" n ',i7,', idigz =',i7,'.')
cbugcbug abug = 'Entered aptbdiv to divide "a" by "b" to get "c", "r".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = '"a" = 20, "b" = 16, "c" = 1 and "r" = 4."'
cbugcbug write ( 3, 9100) abug
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9301) (n, idiga(n), n = 1, ndiga)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9302) (n, idigb(n), n = 1, ndigb)
cbugcbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug 9701 format ('kbug =',i3,' nit=',i3)
cbugcbug abug = 'Initializing kbug = 1, nit = 0'
cbugcbug write ( 3, 9100) abug
cbugcbug kbug = 1
cbugcbug nit = 0
cbugcbug write ( 3, 9701) kbug, nit
cbugcbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
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
if ((ndigb .lt. 1) .or.
& ((ndigb .eq. 1) .and. (idigb(1) .eq. 0))) then
nerr = 4
go to 210
endif
if (ndigb .ge. 1) then
do n = 1, ndigb
if (idigb(n) .lt. 0) then
nerr = 5
go to 210
endif
enddo
endif
if (ndiga .gt. ndigm) then
nerr = 6
go to 210
endif
if (ndigb .gt. ndigm) then
nerr = 7
go to 210
endif
c.... See if big integers "a" and "b" are both zero. If so, error.
if ((ndiga .eq. 1) .and.
& (ndigb .eq. 1) .and.
& (idiga(1) .eq. 0) .and.
& (idigb(1) .eq. 0)) then
nerr = 8
go to 210
endif
cbugcbugc***DEBUG begins.
cbugcbug abug = 'No errors found.'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
c.... Zero out temporary working arrays.
ndigbw = 1
ndigqw = 1
ndigrw = 1
ndigy = 1
ndigz = 1
do n = 1, ndigm
idigbw(n) = 0
idigqw(n) = 0
idigrw(n) = 0
enddo
cbugcbugc***DEBUG begins.
cbugcbug abug = 'Zeroed out "bw", "qw", "rw", "y" and "z".'
cbugcbug write ( 3, 9100) abug
cbugcbug kbug = 7
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9303) (n, idigbw(n), n = 1, ndigbw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9305) (n, idigqw(n), n = 1, ndigqw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9307) (n, idigrw(n), n = 1, ndigrw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9308) (n, idigy(n) , n = 1, ndigy )
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9309) (n, idigy(n) , n = 1, ndigz )
cbugcbugc***DEBUG ends.
c.... Compare the sizes of "a" and "b".
if (ndiga .gt. ndigb) then
cbugcbugc***DEBUG begins.
cbugcbug abug = '"a" is bigger than "b".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = 'This MUST happen when "a" = 20, "b" = 16'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
go to 150
elseif (ndigb .gt. ndiga) then
cbugcbugc***DEBUG begins.
cbugcbug abug = '"b" is bigger than "a".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = 'This can not happen when "a" = 20, "b" = 16'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
go to 110
else ! "a" and "b" have same length.
do n = 1, ndiga
if (idiga(n) .gt. idigb(n)) then
cbugcbugc***DEBUG begins.
cbugcbug abug = '"a" is bigger than "b".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = 'This MUST happen when "a" = 20, "b" = 16'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
go to 150
elseif (idigb(n) .gt. idiga(n)) then
cbugcbugc***DEBUG begins.
cbugcbug abug = '"b" is bigger than "a".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = 'This can not happen when "a" = 20, "b" = 16'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
go to 110
endif
enddo
endif
c.... Big integers "a" and "b" are identical. Quotient is 1, remainder is zero.
cbugcbugc***DEBUG begins.
cbugcbug abug = '"a" is equal to "b". Quotient = 1, residual = 0.'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
ndigqw = 1
idigqw(1) = 1
ndigrw = 1
idigrw(1) = 0
cbugcbugc***DEBUG begins.
cbugcbug kbug = 8
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9305) (n, idigqw(n), n = 1, ndigqw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9307) (n, idigrw(n), n = 1, ndigrw)
cbugcbugc***DEBUG ends.
go to 210
c.... Big integer divisor "b" is bigger than "a".
c.... Quotient "qw" (temporary for "q") is zero.
c..... Remainder "rw" (temporary for "r") is "a" (reversed).
110 continue
cbugcbugc***DEBUG begins.
cbugcbug abug = '"a" is smaller than "b". Quotient = 0, residual = "a".'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
ndigqw = 1
idigqw(1) = 0
ndigrw = ndiga
do n = 1, ndigrw
na = ndiga + 1 - n
idigrw(n) = idiga(na)
enddo
cbugcbugc***DEBUG begins.
cbugcbug kbug = 9
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9305) (n, idigqw(n), n = 1, ndigqw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9307) (n, idigrw(n), n = 1, ndigrw)
cbugcbugc***DEBUG ends.
go to 210
c.... Big integer numerator "a" is bigger than divisor "b".
c.... Remainder "rw" is initially "a" (reversed).
150 continue
cbugcbugc***DEBUG begins.
cbugcbug abug = '"a" is bigger than "b". Residual is initially "a".'
cbugcbug write ( 3, 9100) abug
cbugcbug abug = 'Store reversed "a" into "rw".'
cbugcbug write ( 3, 9100) abug
cbugcbugc***DEBUG ends.
ndigqw = ndiga
ndigrw = ndiga
do n = 1, ndigrw
na = ndiga + 1 - n
idigrw(n) = idiga(na)
enddo
c.... Divisor "bw" is initially "b" (reversed).
ndigbw = ndigb
do n = 1, ndigbw
nb = ndigb + 1 - n
idigbw(n) = idigb(nb)
enddo
cbugcbugc***DEBUG begins.
cbugcbug kbug = 10
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9303) (n, idigbw(n), n = 1, ndigbw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9305) (n, idigqw(n), n = 1, ndigqw)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9307) (n, idigrw(n), n = 1, ndigrw)
cbugcbugc***DEBUG ends.
c.... Shift the digits of divisor "bw" right to line up with remainder "rw",
c.... and fill the new left digits with zeros.
ndiff = ndigrw - ndigbw
if (ndiff .gt. 0) then
do n = ndigbw, 1, -1
idigbw(n+ndiff) = idigbw(n)
enddo
do n = 1, ndiff
idigbw(n) = 0
enddo
endif
ndigbw = ndigrw
nq = ndiff + 1
ndigqw = nq
c.... Return here for each iteration until done.
c.... Put "rw" and "bw" in normal order (most significant digit first).
160 call aptbrev (idigrw, ndigrw, nerrbrev)
call aptbrev (idigbw, ndigbw, nerrbrev)
c.... Find the new remainder "z" = "rw" - "bw" (normal order).
call aptbsub (nbase, idigrw, ndigrw, idigbw, ndigbw, idigy, ndigm,
& idigz, ndigz, isign, nerra)
c.... Put "rw", "bw" and "z" in reversed order (least significant digit first).
call aptbrev (idigrw, ndigrw, nerrbrev)
call aptbrev (idigbw, ndigbw, nerrbrev)
call aptbrev (idigz, ndigz, nerrbrev)
cbugcbugc***DEBUG begins.
cbugcbug nit = nit + 1
cbugcbugc***DEBUG ends.
c.... If the new remainder "z" is positive, move it from "z" to "rw",
c.... and subtract "bw" again.
if (isign .eq. 1) then
idigqw(nq) = idigqw(nq) + 1
ndigrw = ndigz
do n = 1, ndigrw
idigrw(n) = idigz(n)
enddo
go to 160
endif
c.... If the new remainder is negative, ignore it and
c.... shift divisor "bw" left one digit.
if (isign .eq. -1) then
ndigbw = ndigbw - 1
do n = 1, ndigbw
idigbw(n) = idigbw(n+1)
enddo
idigbw(ndigbw+1) = 0
if (ndigbw .gt. 0) then
nq = nq - 1
cbugcbugc***DEBUG begins.
cbugcbug 9932 format (/ ' nq=',i3)
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9932) nq
cbugcbugc***DEBUG ends.
if (nq .le. 0) go to 210
go to 160
endif
endif ! Can not subtract "rw" from "r".
210 continue
cbugcbugc***DEBUG begins.
cbugcbug 9321 format ('nerr=',i3)
cbugcbug abug = 'Jumped to 210. Found an error, or found right answer.'
cbugcbug write ( 3, 9100) abug
cbugcbug write ( 3, 9321) nerr
cbugcbugc***DEBUG ends.
if (nerr .eq. 0) then
c.... Remove trailing zeros from quotient "qw".
if (ndigqw .gt. 1) then
ndigqwx = 1
do n = 1, ndigqw
if (idigqw(n) .gt. 0) ndigqwx = n
enddo
ndigqw = ndigqwx
endif
cbugcbugc***DEBUG begins.
cbugcbug kbug = 20
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9303) (n, idigbw(n), n = 1, ndigbw)
cbugcbugc***DEBUG ends.
c.... Remove trailing zeros from remainder "rw".
if (ndigrw .gt. 1) then
ndigrwx = 1
do n = 1, ndigrw
if (idigrw(n) .gt. 0) ndigrwx = n
enddo
ndigrw = ndigrwx
endif
c.... Move quotient "qw" (reversed) to "q" (normal).
ndigq = ndigqw
do n = 1, ndigq
nq = ndigq + 1 - n
idigq(nq) = idigqw(n)
enddo
c.... Move remainder "rw" (reversed) to "r" (normal).
ndigr = ndigrw
do n = 1, ndigr
nr = ndigr + 1 - n
idigr(nr) = idigrw(n)
enddo
endif
c.... Zero out working arrays.
do n = 1, ndigm
idigqw(n) = 0
idigrw(n) = 0
idigbw(n) = 0
idigy(n) = 0
idigz(n) = 0
enddo
if (nerr .ne. 0) then
ndigq = 0
ndigr = 0
do n = 1, ndigm
idigq(n) = -999999
idigr(n) = -999999
enddo
endif
cbugcbugc***DEBUG begins.
cbugcbug kbug = 21
cbugcbug write ( 3, 9701) kbug, nit
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9303) (n, idigbw(n), n = 1, ndigbw)
cbugcbugc***DEBUG ends.
cbugc***DEBUG begins.
cbug 9915 format (/ 'aptbdiv results: nerr=',i2,', ndigqw =',i7,',',
cbug & ' ndigrw =',i7,'.')
cbug write ( 3, 9900)
cbug write ( 3, 9915) nerr, ndigqw, ndigrw
cbug write ( 3, 9900)
cbug write ( 3, 9305) (n, idigqw(n), n = 1, ndigqw)
cbug write ( 3, 9900)
cbug write ( 3, 9307) (n, idigrw(n), n = 1, ndigrw)
cbugc***DEBUG ends.
return
c.... End of subroutine aptbdiv. (+1 line.)
end
UCRL-WEB-209832