!*************************************************************
!*    Interpolate a function F(x) by continuous fractions    *
!* --------------------------------------------------------- *
!* SAMPLE RUN:                                               *
!* (Interpolate function e(x) between x=0 and x=2)           *
!*                                                           *
!* Number of points: 3                                       *
!* X( 0), Y( 0): 0 1                                         *
!* X( 1), Y( 1): 1 2.71828                                   *
!* X( 2), Y( 2): 2 7.38906                                   *
!*                                                           *
!* Coefficients D(K):                                        *
!* D(0) =      1.000000                                      *
!* D(1) =      0.581977                                      *
!* D(2) =     -3.718271                                      *
!*                                                           *
!* X = 1.5                                                   *
!*                                                           *
!* For X =   1.5  Y =   4.351909                             *
!*                                                           *
!* --------------------------------------------------------- *
!* Ref.: "Méthodes de calcul numérique, Tome 2 By Claude     *
!*        Nowakowski, PSI Edition, 1984" [BIBLI 04].         *
!*                                                           *
!*                         F90 Release By J-P Moreau, Paris. *
!*                                (www.jpmoreau.fr)          *
!*************************************************************
Program Confract

real*8, parameter :: SIZE = 50

integer K,L,M,N,N1
real*8  X(SIZE), Y(SIZE), D(SIZE)
real*8  DD, DL, S, XX

  print *,' '
  write(*,10,advance='no'); Read *, N1

  N = N1 - 1
  M = N

!Read data from screen
  do K = 0, N
    Write(*,20,advance='no') K, K; Read *, X(K), Y(K)
  end do

!Calculate coefficients D(K)
  do K = 0, N 
    D(K) = Y(K)
  end do

  do L = 1, M
    do K = L, N
      DD = (X(K) - X(L - 1)) / (D(K) - D(L - 1))
      IF (K.ne.L) THEN
        D(K) = DD
      ELSE
        DL = DD
      end if
    end do
    D(L) = DL
  end do

!print coefficients
  print *,' '
  print *,' Coefficients D(K):'
  do K = 0, N
    Write(*,40)  K, D(K)
  end do

!Interpolate for X=XX
  print *,' '
  write(*,50,advance='no'); Read *, XX

!Evaluate continuous fraction
  S = (XX - X(N - 1)) / D(N)
  do K = N - 1, 1, -1
    S = (XX - X(K - 1)) / (D(K) + S)
  end do

  S = S + D(0)

  print *,' '
  Write(*,60)  XX, S
  print *,' '

  stop

10 format(' Number of points: ')
20 format(' X(',I2,'), Y(',I2,'): ')
40 format(' D(',I2,') = ', F10.6)
50 format(' X= ')
60 format(' For X= ',F10.6,'  Y= ',F10.6)

END

!end of file confract.f90