subroutine aptbrtn (nroot, nbase, idiga, ndiga, idmax, iemax,
& idigs, idigt, idigu, idigv, idigw, idigx,
& idigy, idigz, idigw1, idigw2, idigw3,
& idigbw, ndigm, idigb, ndigb, isign, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBRTN
c
c call aptbrtn (nroot, nbase, idiga, ndiga, idmax, iemax,
c & idigs, idigt, idigu, idigv, idigw, idigx,
c & idigy, idigz, idigw1, idigw2, idigw3,
c & idigbw, ndigm, idigb, ndigb, isign, nerr)
c
c Version: aptbrtn Updated 2006 May 12 13:40.
c aptbrtn Originated 2005 June 1 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the approximate nroot'th root of the big integer "a",
c stored as the array of non-negative base nbase digits idiga,
c of length ndiga, and store the result as the array of
c non-negative base nbase digits idigb, of length ndigb.
c The nroot'th power "s" of "b", and the difference "t" between
c "s" and "a" are also returned.
c
c Integer arrays idigs thru idigz and idigbs are needed for
c intermediate results. Each must have its own memory space,
c at least of length ndigm.
c
c The digits of idiga and idigbw are the base nbase digits
c representing the non-negative decimal integers ideca and idecb,
c in order from most significant to least significant, using the
c 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^nroot.
c See aptbrev to reverse the order of the digits.
c
c Note: the nroot'th root is the (1 / nroot)'th power.
c
c Flag nerr indicates any input error.
c
c See aptbtod, aptdtob, aptbadd, aptbsub, aptbmul, aptbdiv,
c aptbrev, aptbpow, aptbfac.
c
c Method: The method used here is to calculate the root directly if the
c integer value of "a" will fit in one machine word, or to
c iterate on an initial guess b to find a better guess b':
c b' = ((n - 1) * b + (a / b^(n-1))) / n, or
c b' - b = (a - b^n) / (a * b^(n-1)).
c stopping when b' = b. Due to numerical truncation and integer
c round-off, the result in only accurate to within 1 or 2 digits.
c The initial guess is obtained from the floating point value of
c b, if it fits in a machine word, or from the log of b.
c
c Input: nbase, idiga, ndiga, idigs, idigt, idigu, idigv, idigw, idigx,
c idigy, idigz, idigbw, ndigm.
c
c Output: idigb, ndigb, isign, nerr.
c
c Calls: aptbadd, aptbmul, aptbrev, aptdtob, aptbdiv, aptbpow
c
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 be positive.
c
c idigb Output See idiga.
c Big integer "b" is the nroot'th root of "a", within
c a digit or two, depending on integer truncation.
c
c idigbw Input See idiga. Temporary working space for "b".
c Must have its own working space, ndigm or more.
c
c idigs Output See idiga. Array idigs is returned with the value
c of the nroot'th power of "b", which if accurate,
c would be equal to "a". The number of digits in
c array idigs is returned as idigu(1).
c
c idigt Output See idiga. Array idigt is returned with the
c absolute value of the difference between the values
c of "a" and the nroot'th power of "b". The actual
c sign of the difference is returned as isign.
c The number of digits in array idigt is returned
c as idigu(2).
c
c idigu Output The number of digits in array idigs is returned as
c idigu(1). The number of digits in array idigt
c is returned as idigu(2).
c
c idigv-z Input See idiga. Arrays idigs-z must not be idiga or
c idigbw, but must each have its own memory space,
c which must be at least ndiga.
c
c isign Output Indicates the sign of the difference between the values
c of "a" and the nroot'th power of "b".
c
c nbase Input The number base for which the digits of integer arrays
c idiga, idigw and idigbw are the coefficients of
c the powers of nbase, in order from most to least
c significant.
c
c ndiga Input The length of the integer array idiga.
c
c ndigb Output The length of the integer array idiga. Memory space
c must be at least ndiga.
c
c ndigm Input The maximum number of words allowed in integer arrays
c idigs thru idigz. Memory space for idigbw must
c be at least ndigm, and ndigm must be at least as
c big as nroot * ndiga.
c
c nerr Output Indicates an input or a result error, if not zero.
c 1 if nroot is not positive.
c 2 if nbase is less than 2.
c 3 if ndiga is not positive.
c 4 if any digits idiga are negative.
c 5 if ndigm is less than ndiga.
c 6 if the root can not be found, too close to 1.
c 7 if iteration fails to converge, but loops.
c 8 if subroutine aptbpow fails.
c 9 if subroutine aptbdiv fails.
c 10 if subroutine aptdtob fails.
c 11 if subroutine aptbmul fails.
c 14 if subroutine aptbadd or aptbsub fails.
c
c nroot Input Specifies that the nroot'th or the (1 / nroot)'th
c power of big number "a" is to be found.
c Must be positive.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
integer idiga(1) ! Input "a". Find a^(1/nroot).
integer idigb(1) ! Output "b" = a^(1/nroot).
integer idigbw(1) ! Temporary for "b" = a^(1/nroot).
integer idigs(1) ! Residual, x / z. Output: b^nroot.
integer idigt(1) ! t = a - s = a - b^nroot.
integer idigu(1) ! u = b^(nroot-1). Output Ns, Nt.
integer idigv(1) ! a / u = a / b^(nroot-1). May be zero.
integer idigw(1) ! Temporary, then b * (nroot - 1).
integer idigx(1) ! Temporary, then nroot - 1.
integer idigy(1) ! New b.
integer idigz(1) ! Temporary: nroot. Output: y - b.
c.... Local variables.
integer irootx(10) ! Last ten guesses of "b(1)".
cbugc***DEBUG begins.
cbug character*8 asub ! Subroutine name.
cbug character*16 avarx ! DEBUG label.
cbug character*80 action ! Action label.
cbug avarx = 'NO ARRAY'
cbug action = 'UNDEFINED'
cbug 9900 format (/)
cbug 9901 format (/ 'aptbrtn finding array idigbw, root',i7,
cbug & ' of array idiga' /
cbug & ' both in number base ',i7,'.' /
cbug & ' Length of idiga is ',i7,' with max array length ',i7,'.')
cbug 9902 format ('nit=',i5,2x,a,2x,i6,' = ',25i1)
cbug write ( 3, 9901) nroot, nbase, ndiga, ndigm
cbug write ( 3, 9900)
cbug nit = 0
cbug avarx = 'a input'
cbug do nin = 1, ndiga, 25
cbug nlast = min (nin + 24, ndiga)
cbug write ( 3, 9902) nit, avarx, nin,
cbug & (idiga(n), n = nin, nlast)
cbug enddo
cbugc***DEBUG ends.
c.... Initialize.
c.... Change to least significant digit first.
call aptbrev (idiga, ndiga, nerrb)
c.... Test for input errors.
nerr = 0
if (nroot .le. 0) then
nerr = 1
go to 210
endif
if (nbase .lt. 2) then
nerr = 2
go to 210
endif
if ((ndiga .le. 0) .or.
& ((ndiga .eq. 1) .and. (idiga(1) .eq. 0))) then
nerr = 3 ! Input "a" is zero.
go to 210
endif
ndigax = 0
do n = 1, ndiga
if (idiga(n) .gt. 0) ndigax = n
enddo
if (ndigax .le. 0) then
nerr = 3 ! Input "a" is zero.
go to 210
endif
do n = 1, ndiga
if (idiga(n) .lt. 0) then
nerr = 4 ! Input "a" has negative digits.
go to 210
endif
enddo
if (ndigm .lt. ndiga) then
nerr = 5 ! Not enough memory for "a".
go to 210
endif
c.... Initialize.
ndigbw = 1
do n = 1, ndigm
idigbw(n) = 0
idigs(n) = 0
idigt(n) = 0
idigu(n) = 0
idigv(n) = 0
idigw(n) = 0
idigx(n) = 0
idigy(n) = 0
idigz(n) = 0
enddo
c.... See if nroot = 1, so "b" = "a", "s" = "a" and "r" = "a" - "s" = 0.
if (nroot .eq. 1) then
ndigbw = ndiga
do n = 1, ndigbw
idigbw(n) = idiga(n)
idigs(n) = idiga(n)
enddo
ndigs = ndiga
ndigt = 1
idigt(1) = 0
idigu(1) = ndigs ! The number of digits in idigs.
idigu(2) = ndigt ! The number of digits in idigt.
isign = 1
go to 210
endif
c.... Find the base 10 logarithm of the root.
c.... Change to most significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbtod (nbase, idiga, ndiga, idmax, iemax,
& ideca, fdeca, floga, nerra)
c.... Change to least significant digit first.
call aptbrev (idiga, ndiga, nerrb)
if (floga .eq. -99.0) then
nerr = 3
go to 210
endif
fnroot = nroot
flogb = floga / fnroot
c.... Find the approximate value of the root.
call aptltob (flogb, nbase, ndigm, idigbw, ndigbw, nerra)
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb) ! DO NOT NEED.
c.... Find the integer value of the root, if possible.
c.... Change to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb) ! DO NOT NEED.
call aptbtod (nbase, idigbw, ndigbw, idmax, iemax,
& idecb, fdecb, flogb, nerra)
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
c.... If the integer value is accurate, use it as the root "b".
c.... Find "s", the nroot'th power of the root "b", and the difference "t"
c.... between "a" and "s".
if (flogb .le. 16) then
cbugcbugc***DEBUG begins.
cbugcbug 9301 format ('APTBRTN DEBUG. Used integer value directly.')
cbugcbug write ( 3, 9301)
cbugcbugc***DEBUG ends.
c.... Change to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbpow (nroot, nbase, idigbw, ndigbw, idigw, idigw1,
& ndigm, idigs, ndigs, nerra)
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbrev (idigs, ndigs, nerrb) ! DO NOT NEED.
c.... Change to most significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbrev (idigs, ndigs, nerrb) ! DO NOT NEED.
call aptbsub (nbase, idiga, ndiga, idigs, ndigs,
& idigz, ndigm, idigt, ndigt, isign, nerra)
c.... Change to least significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbrev (idigs, ndigs, nerrb)
call aptbrev (idigt, ndigt, nerrb)
idigu(1) = ndigs
idigu(2) = ndigt
cbugcbugc***DEBUG begins.
cbugcbug 9333 format ('aptbrtn DEBUG. ndigs =',i5,' ndigt =',i5)
cbugcbug write ( 3, 9333) ndigs, ndigt
cbugcbugc***DEBUG ends.
go to 210
endif
c.... Initialize for finding the root by iteration.
do n = 1, 10
irootx(n) = -999999
enddo
c.... Find the nroot'th root of "a" by iteration.
nit = 0
120 nit = nit + 1
cbugcbugc***DEBUG begins.
cbugcbug 9931 format ('APTBRTN NBUG =',i3)
cbugcbug 9932 format (' # abstuvwxy =',9i5)
cbugcbug 9933 format (/ a80)
cbugcbug 9971 format ('nit=',i5,2x,a,2x,i6,' = ',25i1)
cbugcbug nbug = 1
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Start with a, initial guess b.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'a'
cbugcbug call aptbrev (idiga, ndiga, nerrb)
cbugcbug do nin = 1, ndiga, 25
cbugcbug nlast = min (nin + 24, ndiga)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idiga(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idiga, ndiga, nerrb)
cbugcbug
cbugcbug avarx = 'b'
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbugc***DEBUG ends.
c.... Find "v" = b^(nroot - 1). Use work array "w".
if (nroot .eq. 2) then
ndigv = ndigbw
do n = 1, ndigv
idigv(n) = idigbw(n)
enddo
cbugcbugc***DEBUG begins.
cbugcbug nbug = 21
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found v = b^(N-1), using work array v.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'b arg^(1/N)'
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbug
cbugcbug avarx = 'v b^(N-1)'
cbugcbug call aptbrev (idigv, ndigv, nerrb)
cbugcbug do nin = 1, ndigv, 25
cbugcbug nlast = min (nin + 24, ndigv)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigv(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idigv, ndigv, nerrb)
cbugcbugc***DEBUG ends.
else
npow = nroot - 1
c.... Change to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbpow (npow, nbase, idigbw, ndigbw, idigw, idigw1,
& ndigm, idigv, ndigv, nerra)
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptbpow1'
cbugcbug 9905 format ('Called ',a,'. nerra =',i5)
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug 9972 format ('aptbrtn nerra =',i3,' ndigv =',i20)
cbugcbug write ( 3, 9972) nerra, ndigv
cbugcbug nbug = 3
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found v = b^(N-1), using work array w.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'b arg^(1/N)'
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'v b^(N-1)'
cbugcbug do nin = 1, ndigv, 25
cbugcbug nlast = min (nin + 24, ndigv)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigv(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbrev (idigv, ndigv, nerrb)
if (nerra .ne. 0) then
nerr = 8
go to 210
endif
endif
c.... Find "s" = b^n = b * b^(n-1), which should be equal to "a".
c.... Change to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbmul (nbase, idigbw, ndigbw, idigv, ndigv,
& idigz, ndigm, idigs, ndigs, nerra)
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptbpow2'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 12
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found s = b^N, stored length in u(1).'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'b a^(1/N)'
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 's b^N ~ a'
cbugcbug do nin = 1, ndigs, 25
cbugcbug nlast = min (nin + 24, ndigs)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigs(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbrev (idigs, ndigs, nerrb)
c.... Find "t" = |"a" - "s"| = |a - b^nroot|.
c.... Change to most significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbrev (idigs, ndigs, nerrb)
call aptbsub (nbase, idiga, ndiga, idigs, ndigs,
& idigz, ndigm, idigt, ndigt, isign, nerra)
ndigu = 2
idigu(1) = ndigs ! Number of digits in "s" = b^nroot.
idigu(2) = ndigt ! Number of digits in "t" = "a" - "s".
cbugcbugc***DEBUG begins.
cbugcbug 9970 format ('nit=',i5,2x,a,2x,i6,' = ',25i5)
cbugcbug asub = 'aptbsub2'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 13
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found t = |a - s|. Stored sizes in u(2).'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'a arg'
cbugcbug do nin = 1, ndiga, 25
cbugcbug nlast = min (nin + 24, ndiga)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idiga(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 's b^N'
cbugcbug do nin = 1, ndigs, 25
cbugcbug nlast = min (nin + 24, ndigs)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigs(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 't |a-b^N|'
cbugcbug do nin = 1, ndigt, 25
cbugcbug nlast = min (nin + 24, ndigt)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigt(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'u Nt, Ns '
cbugcbug do nin = 1, ndigu, 25
cbugcbug nlast = min (nin + 24, ndigu)
cbugcbug write ( 3, 9970) nit, avarx, nin,
cbugcbug & (idigu(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbrev (idigs, ndigs, nerrb)
call aptbrev (idigt, ndigt, nerrb)
c.... Test for convergence. Stop if "t" is zero.
if ((ndigt .eq. 0) .or.
& (ndigt .eq. 1) .and. (idigt(1) .eq. 0)) then
go to 210
endif
c.... Find the big number array nroot, store in "x".
c.... Change to most significant digit first.
call aptbrev (idigx, ndigx, nerrb)
call aptdtob (nroot, nbase, ndigm, idigx, ndigx, nerra)
cbugcbugc***DEBUG begins.
cbugcbugc***NOTE: THIS HAD BEEN RETURNING "X" IN REVERSE ORDER.
cbugcbug 9903 format ('aptbrtn DEBUG. nbug=',i5)
cbugcbug 9904 format (' x n =',i7,', idigx =',i7,'.')
cbugcbug nbig = 101
cbugcbug write ( 3, 9903) nbug, nroot, nbase
cbugcbug write ( 3, 9900)
cbugcbug write ( 3, 9904) (n, idigx(n), n = 1, ndigx)
cbugcbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptdtob1'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 62
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found x = N.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'x N'
cbugcbug do nin = 1, ndigx, 25
cbugcbug nlast = min (nin + 24, ndigx)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigx(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
c***NOTE: THIS HAD BEEN PUTTING "X" IN NORMAL ORDER.
call aptbrev (idigx, ndigx, nerrb)
if (nerra .ne. 0) then
nerr = 10
go to 210
endif
c.... Find "w" = "x" * "v" = N * b^(N-1).
c.... Change to most significant digit first.
c***NOTE: THIS HAD BEEN PUTTING "X" IN REVERSE ORDER.
call aptbrev (idigx, ndigx, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbmul (nbase, idigx, ndigx, idigv, ndigv,
& idigz, ndigm, idigw, ndigw, nerra)
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptbmul'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 7
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found w = x * v = N * b^(N-1).'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'x N'
cbugcbug do nin = 1, ndigx, 25
cbugcbug nlast = min (nin + 24, ndigx)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigx(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'b arg^(1/N)'
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'w N*b^(N-1)'
cbugcbug do nin = 1, ndigw, 25
cbugcbug nlast = min (nin + 24, ndigw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigw(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idigx, ndigx, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbrev (idigw, ndigw, nerrb)
if (nerra .ne. 0) then
nerr = 11
go to 210
endif
c.... Find "v" = t / w = |a - b^N| / (N * b^(N-1)). Use work arrays y, z.
c.... Change to most significant digit first.
call aptbrev (idigt, ndigt, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbdiv (nbase, idigt, ndigt, idigw, ndigw,
& idigy, idigz, idigw1, idigw2, idigw3, ndigm,
& idigv, ndigv, idigx, ndigx, nerra)
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptbdiv1'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 5
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found v = t/w, x = residue, with y and z working.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'a arg'
cbugcbug do nin = 1, ndiga, 25
cbugcbug nlast = min (nin + 24, ndiga)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idiga(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 's b^N'
cbugcbug do nin = 1, ndigs, 25
cbugcbug nlast = min (nin + 24, ndigs)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigs(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 't |a-b^N|'
cbugcbug do nin = 1, ndigt, 25
cbugcbug nlast = min (nin + 24, ndigt)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigt(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'w N * b^(N-1)'
cbugcbug do nin = 1, ndigw, 25
cbugcbug nlast = min (nin + 24, ndigw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'v t / w'
cbugcbug do nin = 1, ndigv, 25
cbugcbug nlast = min (nin + 24, ndigv)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigv(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'x remainder'
cbugcbug do nin = 1, ndigx, 25
cbugcbug nlast = min (nin + 24, ndigx)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigx(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idigt, ndigt, nerrb)
call aptbrev (idigw, ndigw, nerrb)
call aptbrev (idigy, ndigy, nerrb)
call aptbrev (idigz, ndigz, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbrev (idigx, ndigx, nerrb)
if (nerra .ne. 0) then
nerr = 9
cbugcbugc***DEBUG begins.
cbugcbug 9974 format ('aptbdiv nerra =',i3,' ndigu =',i20)
cbugcbug 9973 format ('APTBDIV BAD. nerr =',i3,' u is zero.')
cbugcbug write ( 3, 9974) nerra, ndigu
cbugcbug write ( 3, 9973) nerr
cbugcbugc***DEBUG ends.
go to 210
endif
c.... See if correction to "bw" is zero. If so, done.
if ((ndigv .eq. 0) .or.
& ((ndigv .eq. 1) .and. (idigv(1) .eq. 0))) then
go to 210
endif
c.... Find "y" = "bw" + "v". This is a better guess of "bw" = a^(1/n).
c.... Work array "z".
c.... Change to most significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigv, ndigv, nerrb)
if (isign .eq. 1) then
call aptbadd (nbase, idigbw, ndigbw, idigv, ndigv,
& idigz, ndigm, idigy, ndigy, nerra)
else
call aptbsub (nbase, idigbw, ndigbw, idigv, ndigv,
& idigz, ndigm, idigy, ndigy, isign, nerra)
endif
cbugcbugc***DEBUG begins.
cbugcbug asub = 'aptbdiv2'
cbugcbug write ( 3, 9905) asub, nerra
cbugcbug nbug = 9
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found y = b + v = new b.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'b A^(1/n)'
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'v corr b'
cbugcbug do nin = 1, ndigv, 25
cbugcbug nlast = min (nin + 24, ndigv)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigv(n), n = nin, nlast)
cbugcbug enddo
cbugcbug
cbugcbug avarx = 'y new b'
cbugcbug do nin = 1, ndigy, 25
cbugcbug nlast = min (nin + 24, ndigy)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigy(n), n = nin, nlast)
cbugcbug enddo
cbugcbugc***DEBUG ends.
c.... Change to least significant digit first.
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigv, ndigv, nerrb)
call aptbrev (idigy, ndigy, nerrb)
if (nerra .ne. 0) then
nerr = 14
go to 210
endif
c.... Replace "bw" with "y", iterate.
ndigbw = ndigy
do n = 1, ndigy
idigbw(n) = idigy(n)
enddo
cbugcbugc***DEBUG begins.
cbugcbug nbug = 11
cbugcbug write ( 3, 9931) nbug
cbugcbug action = 'Found next guess for b, using y.'
cbugcbug write ( 3, 9933) action
cbugcbug write ( 3, 9932) ndiga, ndigbw, ndigs, ndigt,
cbugcbug & ndigu, ndigv, ndigw, ndigx, ndigy
cbugcbug
cbugcbug avarx = 'y new'
cbugcbug call aptbrev (idigy, ndigy, nerrb)
cbugcbug do nin = 1, ndigy, 25
cbugcbug nlast = min (nin + 24, ndigy)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigy(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idigy, ndigy, nerrb)
cbugcbug
cbugcbug avarx = 'b new'
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbug do nin = 1, ndigbw, 25
cbugcbug nlast = min (nin + 24, ndigbw)
cbugcbug write ( 3, 9971) nit, avarx, nin,
cbugcbug & (idigbw(n), n = nin, nlast)
cbugcbug enddo
cbugcbug call aptbrev (idigbw, ndigbw, nerrb)
cbugcbugc***DEBUG ends.
c.... See if iteration is in a loop.
if (nit .gt. 1000) then
do n = 1, 10
if (idigbw(1) .eq. irootx(n)) then
nerr = 7
go to 210
endif
enddo
isave = 1 + mod (isave, 10)
irootx(isave) = idigbw(1)
endif
go to 120
210 continue
c.... Change "a", "bw", "s" and "t" to most significant digit first.
call aptbrev (idiga, ndiga, nerrb)
call aptbrev (idigbw, ndigbw, nerrb)
call aptbrev (idigs, ndigs, nerrb)
call aptbrev (idigt, ndigt, nerrb)
c.... Store working array "bw" into "b".
ndigb = ndigbw
do n = 1, ndigb
idigb(n) = idigbw(n)
enddo
if (nerr .ne. 0) then
ndigs = 0
ndigt = 0
ndigu = 0
ndigb = 0
do n = 1, ndigm
idigs(n) = -999999
idigt(n) = -999999
idigu(n) = -999999
idigb(n) = -999999
enddo
endif
cbugc***DEBUG begins.
cbug 9911 format (/ 'aptbrtn results: nerr=',i2,', ndigbw =',i7,'.')
cbug write ( 3, 9911) nerr, ndigbw
cbug
cbug write ( 3, 9900)
cbug avarx = 'a arg "a"'
cbug do nin = 1, ndiga, 25
cbug nlast = min (nin + 24, ndiga)
cbug write ( 3, 9902) nit, avarx, nin,
cbug & (idiga(n), n = nin, nlast)
cbug enddo
cbug
cbug write ( 3, 9900)
cbug avarx = 'b a^(1/N)'
cbug do nin = 1, ndigbw, 25
cbug nlast = min (nin + 24, ndigbw)
cbug write ( 3, 9902) nit, avarx, nin,
cbug & (idigbw(n), n = nin, nlast)
cbug enddo
cbug
cbug write ( 3, 9900)
cbug avarx = 's b^N ~ a'
cbug do nin = 1, ndigs, 25
cbug nlast = min (nin + 24, ndigs)
cbug write ( 3, 9902) nit, avarx, nin,
cbug & (idigs(n), n = nin, nlast)
cbug enddo
cbug write ( 3, 9900)
cbug
cbug avarx = 'a - b^N'
cbug do nin = 1, ndigt, 25
cbug nlast = min (nin + 24, ndigt)
cbug write ( 3, 9902) nit, avarx, nin,
cbug & (idigt(n), n = nin, nlast)
cbug enddo
cbugc***DEBUG ends.
return
c.... End of subroutine aptbrtn. (+1 line.)
end
UCRL-WEB-209832