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