/************************************************************
*    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].         *
*                                                           *
*                         C++ Release By J-P Moreau, Paris. *
*                                 (www.jpmoreau.fr)         *
************************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  50

int K,L,M,N,N1;
double X[SIZE], Y[SIZE], D[SIZE];
double DD, DL, S, XX;


void main() {

  printf("\n Number of points: "); scanf("%d", &N1);

  N = N1 - 1;
  M = N;

// Read data from screen
  for (K=0; K<=N; K++) {
    printf(" X(%d), Y(%d): ",K, K); scanf("%lf %lf", &X[K], &Y[K]);
  }

// Calculate coefficients D(K)
  for (K=0; K<=N; K++)  D[K] = Y[K];

  for (L=1; L<=M; L++) {
	for (K=L; K<=N; K++) {
      DD = (X[K] - X[L - 1]) / (D[K] - D[L - 1]);
      if (K != L)
        D[K] = DD;
      else
        DL = DD;
    }
    D[L] = DL;
  }

// print coefficients
  printf("\n Coefficients D(K):\n");
  for (K=0; K<=N; K++)
    printf(" D(%d) = %10.6f\n", K, D[K]);

// Interpolate for X=XX
  printf("\n X = "); scanf("%lf", &XX);

// Evaluate continuous fraction
  S = (XX - X[N - 1]) / D[N];
  for (K=N-1; K>0; K--)
    S = (XX - X[K - 1]) / (D[K] + S);

  S += D[0];

  printf(" For X= %10.6f  Y= %10.6f\n\n", XX, S);

}

// end of file confract.cpp