subroutine aptcofr (x, ncmax, tol, kc, inum, iden, xc, nc, nerr) c. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTCOFR c c call aptcofr (x, ncmax, tol, kc, inum, iden, xc, nc, nerr) c c Version: aptcofr Updated 2001 July 16 14:00. c aptcofr 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 continued c fraction equivalent of x, with a relative accuracy of tol or c better, and the numerator inum(n), denominator iden(n), and c partial sum xc(n) of the n'th convergent, for n = 1, nc. c c The continued fraction form is: c x = kc(1) + 1/(kc(2) + 1/(kc(3) + 1/(kc(4) + ...))) c c and is obtained by starting with the integer <= x, then the c integer <= the reciprocal of the remainder, etc. c All coefficients will have the same sign as x, and none will c be zero except possibly the first. c c The continued fraction convergent xc(n) is given by: c xc(n) = kc(1) + 1/(kc(2) + 1/(kc(3) + ... + 1 / kc(n))) c = inum(n) / iden(n). c c and is obtained by forward recursion: c inum(0) = 1, inum(1) = kc(1), c inum(n) = kc(n) * inum(n-1) + inum(n-2), n = 2, nc, c iden(0) = 0, iden(1) = 1, c iden(n) = kc(n) * iden(n-1) + iden(n-2), n = 2, nc. c xc(n) = inum(n) / iden(n), n = 1, nc. c c Notes: c c See subroutine aptcofv for the inverse problem (given the c coefficients, find the value, convergents, and partial 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, x is the reciprocal of the c x evaluated from the remaining coefficients. c c Any coefficient and all subsequent coefficients may be replaced c by a single coefficient equal to the x (not necessarily an c integer) evaluated for that and all subsequent coefficients. c c Input: x, ncmax, tol. c c Output: kc, inum, iden, nc, nerr. c c Glossary: c c iden(n) Output The denominator of the n'th convergent of the c continued fraction. x(n) = inum(n) / iden(n). c Size nc. c c inum(n) Output The numerator of the n'th convergent of the c continued fraction. x(n) = inum(n) / iden(n). c Size nc. c c kc(n) Output The n'th integer coefficient of the continued fraction c equal to x, n = 1, nc <= ncmax. c c nc Output The actual number of coefficients kc, numerators inum, c denominators iden and partial sums 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 continued fraction terminates. c Arrays kc, inum, iden and xc must be sized at least c as big as 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 denominator of convergent is 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 x Input The value of the continued fraction c (kc1, kc2, kc3, ...). c c xc(n) Output The n'th partial sum of the continued fraction c (kc1, kc2, kc3, ...). Size nc. c x(n) = inum(n) / iden(n) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. implicit none c.... Arguments. integer iden(1) ! Denominators of convergents. integer inum(1) ! Numerators of convergents. integer kc(1) ! Coeffficients of continued fraction. 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 continued fraction. real xc(1) ! Partial sum of continued fraction. c.... Local variables. integer idenss ! Denominator of a convergent. integer inumss ! Numerator of a convergent. integer n ! Index in arrays kc, inum, iden, xc. real pn ! Numerator of convergent, inum. real qn ! Denominator of convergent, iden. real rem ! Decimal part of xnew. real xint ! The integer part of xnew. real xnew ! Inverse of reminder rem. cbugc***DEBUG begins. cbug 9901 format (/ 'aptcofr finding coeffs of continued fraction.' / 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 do 110 n = 1, ncmax kc(n) = -1000000000 inum(n) = -1000000000 iden(n) = -1000000000 xc(n) = -1.e99 110 continue c.... Test for input errors. if (ncmax .le. 0) then nerr = 1 go to 210 endif xnew = x xint = int (xnew) if (abs (xnew - xint) .ge. 1.0) then nerr = 2 go to 210 endif c.... Find all significant coefficients kc(1) to kc(nc) in the continued c.... fraction expression for x, and the numerators and denominators of the c.... convergents, and the partial sums. nc = 1 xint = int (xnew) kc(1) = xint inumss = 1.0 inum(1) = kc(1) idenss = 0.0 iden(1) = 1.0 xc(1) = inum(1) rem = mod (xnew, 1.0) if (abs (rem) .le. tol) rem = 0.0 if (rem .eq. 0.0) go to 210 xnew = 1.0 / rem do 120 nc = 2, ncmax xint = int (xnew) if (abs (xnew - xint) .ge. 1.0) then nerr = 2 go to 210 endif kc(nc) = xint inum(nc) = kc(nc) * inum(nc-1) + inumss iden(nc) = kc(nc) * iden(nc-1) + idenss if (iden(nc) .eq. 0) then nerr = 3 go to 210 endif pn = inum(nc) qn = iden(nc) xc(nc) = pn / qn inumss = inum(nc-1) idenss = iden(nc-1) rem = mod (xnew, 1.0) if (abs (1.0 - rem) .le. tol) then kc(nc) = kc(nc) + 1 rem = rem - 1.0 endif if (abs (rem) .le. tol) rem = 0.0 if (abs (rem) .le. 0.0) go to 210 if (abs (xc(nc) - xc(nc-1)) .le. tol * abs (xc(nc))) go to 210 if (nc .eq. ncmax) go to 210 xnew = 1.0 / rem 120 continue c=============================================================================== 210 continue cbugc***DEBUG begins. cbug 9902 format (/ 'aptcofr results: nerr = ',i3,' nc = ',i3 ) cbug 9903 format (' n, kc(n) =',i3,i16 / cbug & ' xc,inum,iden=',1pe22.14,2i16) cbug write ( 3, 9902) nerr, nc cbug if (nc .gt. 0) then cbug write ( 3, 9903) (n, kc(n), xc(n), inum(n), iden(n), n = 1, nc) cbug endif cbugc***DEBUG ends. return c.... End of subroutine aptcofr. (+1 line.) end UCRL-WEB-209832