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