subroutine aptiris (x, ncmax, tol, kc, xc, nc, nerr)
c.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTIRIS
c
c call aptiris (x, ncmax, tol, kc, xc, nc, nerr)
c
c Version: aptiris Updated 2002 June 4 14:40.
c aptiris Originated 2000 October 16 12:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the first nc <= ncmax coefficients kc of the
c equivalent of x expressed as an integer plus the sum of
c reciprocals of integers:
c
c x = kc(1) + 1/kc(2) + 1/kc(3) + 1/kc(4) + ...
c
c with a relative accuracy of tol or better,
c and the partial sums xc(n), for n = 1, nc:
c
c xc(n) = kc(1) + 1/kc(2) + 1/kc(3) + ... + 1 / kc(n).
c
c The coefficients kc are obtained by starting with the integer
c nearest x, finding the remainder, then finding the integer
c nearest the reciprocal of the remainder, etc.
c
c Input: x, ncmax, tol.
c
c Output: kc, xc, nc, nerr.
c
c Glossary:
c
c kc(n) Output The n'th integer coefficient of the sum of fractions
c equal to x, n = 1, nc <= ncmax.
c
c nc Output The actual number of coefficients kc and partial sums
c xc calculated.
c
c ncmax Input The maximum number of coefficients kc to be calculated.
c Fewer (nc) will be calculated if additional
c coefficients would not make the accuracy any better
c than tol, or if the sum of fractions terminates.
c Arrays kc and xc must be sized at least as big as
c ncmax.
c
c nerr Output Error indicator. 0 if no errors.
c 1 if ncmax is not positive.
c 2 if x exceeds the largest integer.
c 3 if any kc exceeds the largest integer.
c
c tol Input The estimated relative error in floating point numbers
c on the machine in use. Approximately 1.e-15 for
c 64-bit floating point numbers.
c
c x Input The value of the sum of fractions
c kc(1) + 1 / kc(2) + 1 / kc(3) + ... + 1 / kc(nc).
c
c xc(n) Output The n'th partial sum of the sum of fractions
c kc(1) + 1 / kc(2) + 1 / kc(3) + ... + 1 / kc(n).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
implicit none
c.... Arguments.
integer kc(1) ! Coeffficients of sum of fractions.
integer nc ! Number of calculated coefficients.
integer ncmax ! Max number of calculated coefficients.
integer nerr ! Error flag.
real tol ! Floating point precision.
real x ! Value of sum of fractions.
real xc(1) ! Partial sum of sum of fractions.
c.... Local variables.
integer n ! Index in arrays kc, xc.
integer ngood ! Index of current coefficient.
real rem ! Difference between x and xc(n).
real rint ! The nearest integer to rinv.
real rinv ! Inverse of reminder rem.
real xint ! The nearest integer to x.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptiris finding coeffs of sum of fractions.' /
cbug & ' x =',1pe22.14 /
cbug & ' ncmax =',i3 /
cbug & ' tol =',1pe22.14 )
cbug write ( 3, 9901) x, ncmax, tol
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
ngood = 0
do 110 n = 1, ncmax
kc(n) = -1000000000000000
xc(n) = -1.e99
110 continue
c.... Test for input errors.
if (ncmax .le. 0) then
nerr = 1
go to 210
endif
xint = nint (x)
if (abs (x - xint) .ge. 1.0) then
nerr = 2
go to 210
endif
c.... Find the coefficient kc(1), and the first partial sum xc(1).
nc = 1
kc(1) = xint
xc(1) = xint
rem = x - xc(1)
if (abs (rem) .le. tol * abs (x)) then
ngood = 1
go to 210
endif
rinv = 1.0 / rem
rint = nint (rinv)
cbugcbugc***DEBUG begins.
cbugcbug 9700 format (' rinv=',1pe22.14,' xint=',1pe22.14 /
cbugcbug & ' rem= ',1pe22.14)
cbugcbug write ( 3, 9903) nc, kc(nc), xc(nc)
cbugcbug write ( 3, 9700) rinv, xint, rem
cbugcbugc***DEBUG ends.
if (abs (rinv - rint) .ge. 1.0) then
if (kc(1) .eq. 0) then
nerr = 3
else
ngood = nc
endif
go to 210
endif
c.... Find all significant coefficients kc(2) to kc(nc) in the sum or
c.... fractions expression for x, and the partial sums xc(n) to xc(nc).
do 120 nc = 2, ncmax
kc(nc) = rint
xc(nc) = xc(nc-1) + 1.0 / rint
rem = x - xc(nc)
cbugcbugc***DEBUG begins.
cbugcbug 9701 format (' rinv=',1pe22.14,' rint=',1pe22.14 /
cbugcbug & ' rem= ',1pe22.14)
cbugcbug write ( 3, 9903) nc, kc(nc), xc(nc)
cbugcbug write ( 3, 9701) rinv, rint, rem
cbugcbugc***DEBUG ends.
if (abs (rem) .le. tol * abs (x)) then
ngood = nc
go to 210
endif
rinv = 1.0 / rem
rint = nint (rinv)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9903) nc, kc(nc), xc(nc)
cbugcbug write ( 3, 9701) rinv, rint, rem
cbugcbugc***DEBUG ends.
if (abs (rinv - rint) .ge. 1.0) then
ngood = nc
go to 210
endif
120 continue
ngood = ncmax
c===============================================================================
210 nc = ngood
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptiris results: nerr = ',i3,' nc = ',i3 )
cbug 9903 format (' n, kc(n) =',i3,i22,' xc=',1pe22.14)
cbug write ( 3, 9902) nerr, nc
cbug if (nc .gt. 0) then
cbug write ( 3, 9903) (n, kc(n), xc(n), n = 1, nc)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptiris. (+1 line.)
end
UCRL-WEB-209832