subroutine aptcofv (c, nc, tol, xnum, xden, xc, nerr)
c.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCOFV
c
c call aptcofv (c, nc, tol, xnum, xden, xc, nerr)
c
c Version: aptcofv Updated 2001 July 17 14:20.
c aptcofv Originated 2000 October 20 14:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Given the first nc coefficients c of a continued fraction,
c to find the partial sums xc(n), and the numerators xnum(n) and
c denominators xden(n) of the n'th convergent to xc(n),
c for n = 1, nc.
c
c The n'th partial sum xc(n) is given by:
c xc(n) = c(1) + 1/(c(2) + 1/(c(3) + ... + 1/c(n)))
c = xnum(n) / xden(n).
c
c and is obtained by forward recursion:
c xnum(0) = 1.0, xnum(1) = c(1),
c xnum(n) = c(n) * xnum(n-1) + xnum(n-2),
c xden(0) = 0.0, xden(1) = 1.0,
c xden(n) = c(n) * xden(n-1) + xden(n-2),
c xc(n) = xnum(n) / xden(n), n = 1, nc.
c
c Notes:
c
c The coefficients may have any values. Trailing zeros cause
c the partial sums to oscillate. A final coefficient of 1 may be
c added to the preceding coefficient.
c
c See subroutine aptcofr for the inverse operation (given the
c value, find the coefficients, the convergents and the partial
c sums).
c For a rational number (an integer or a ratio of two integers),
c the coefficients terminate. For an irrational number (a root
c of a polynomial equation), the coefficients repeat in groups.
c For a transcendental number (neither rational or irrational,
c such as pi or e), the coefficients do not repeat in groups, but
c may still form a predictable pattern.
c
c If the first coefficient is zero, xc is the reciprocal of the
c xc evaluated from the remaining coefficients.
c
c Any coefficient and all subsequent coefficients may be replaced
c by a single coefficient equal to the xc evaluated for that and
c all subsequent coefficients.
c
c Input: c, nc, tol.
c
c Output: xc, xnum, xden, nerr.
c
c Glossary:
c
c c(n) Input The n'th coefficient of the continued fraction
c equivalent of xc. Size nc.
c Normally only the first coefficient may be zero, and
c the rest must be positive integers, but you may
c experiment with any desired values, at the risk of
c generating infinities.
c
c nc Input The number of coefficients c.
c
c nerr Output Error indicator.
c 1 if nc = 0.
c 2 if denominator becomes 0.
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 xc(n) Output The n'th partial sum of the continued fraction
c (c(1), c(2), c(3), ... c(n)). Size nc.
c xc(n) = xnum(n) / xden(n)
c
c xden(n) Output The denominator of the n'th convergent of the
c continued fraction. xc(n) = xnum(n) / xden(n).
c Size nc.
c
c xnum(n) Output The numerator of the n'th convergent of the
c continued fraction. xc(n) = xnum(n) / xden(n).
c Size nc.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
implicit none
c.... Arguments.
integer nc ! Number of calculated coefficients.
integer nerr ! Error flag, 1 if infinity occurs.
real c(1) ! Coeffficients of continued fraction.
real tol ! Floating point precision.
real xc(1) ! Partial sum of continued fraction.
real xnum(1) ! Numerator of convergent.
real xnumnew ! Numerator of convergent.
real xden(1) ! Denominator of convergent.
c.... Local variables.
integer n ! Index in array c.
real err ! Rel error in numerator or denominator.
real xnumss ! 2nd preceding numerator.
real xdenss ! 2nd preceding denominator.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcofv finding value of continued fraction.',
cbug & ' nc = ',i3 /
cbug & (' c(',i3,') =',1pe22.14 ))
cbug write ( 3, 9901) nc, (n, c(n), n = 1, nc)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
if (nc .le. 0) then
nerr = 1
go to 210
endif
do 110 n = 1, nc
xc(n) = -1.e99
xnum(n) = -1.e99
xden(n) = -1.e99
110 continue
c.... Find the partial sums of the continued fraction, and the corresponding
c.... numerators and denominators of the convergents..
xnumss = 1.0
xnum(1) = c(1)
xdenss = 0.0
xden(1) = 1.0
xc(1) = c(1)
do 280 n = 2, nc
xnum(n) = c(n) * xnum(n-1) + xnumss
err = tol * (abs (c(n) * xnum(n-1)) + abs (xnumss))
if (abs (xnum(n)) .le. err) xnum(n) = 0.0
xden(n) = c(n) * xden(n-1) + xdenss
err = tol * (abs (c(n) * xden(n-1)) + abs (xdenss))
if (abs (xden(n)) .le. err) xden(n) = 0.0
if (xden(n) .eq. 0.0) then
nerr = 2
go to 210
endif
xc(n) = xnum(n) / xden(n)
xnumss = xnum(n-1)
xdenss = xden(n-1)
280 continue
c===============================================================================
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptcofv results: nerr = ',i3,' nc = ',i3 /
cbug & (' xc,xnum,xden=',1p3e22.14 ))
cbug write ( 3, 9902) nerr, nc, (xc(n), xnum(n), xden(n),
cbug & n = 1, nc)
cbugc***DEBUG ends.
return
c.... End of subroutine aptcofv. (+1 line.)
end
UCRL-WEB-209832