/************************************************************
*    Integration of a discrete function by the weighting    *
*    coefficients Method                                    *
* --------------------------------------------------------- *
* SAMPLE RUN:                                               *
*                                                           *
* Number of points: 7                                       *
*                                                           *
* Input the 7 points, example: 1.5  0.666666                *
*                                                           *
* 0  ?  1 1                                                 *
* 1  ?  1.5 0.666666                                        *
* 2  ?  2 0.5                                               *
* 3  ?  2.5 0.4                                             *
* 4  ?  3 0.333333                                          *
* 5  ?  3.5 0.285714                                        *
* 6  ?  4 0.25                                              *
*                                                           *
* Coefficients A(I):                                        *
* A(0) =  0.146429                                          *
* A(1) =  0.771429                                          * 
* A(2) =  0.096429                                          *
* A(3) =  0.971429                                          * 
* A(4) =  0.096429                                          *
* A(5) =  0.771429                                          *
* A(6) =  0.146429                                          *
*                                                           *
* Integral of F(X) from  1.00 to  4.00:                     *
*   I =     1.386657                                        *
*                                                           *
* --------------------------------------------------------- *
* Reference: "Méthodes de calcul numérique - Tome 1 By      *
*             Claude Nowakowski, PS1 1981" [BIBLI 07].      *
*                                                           *
*                       C++ Release By J-P Moreau, Paris.   *
*                              (www.jpmoreau.fr)            *
*************************************************************
See explanations in file dinteg.txt
-----------------------------------  */
#include <stdio.h>

#define  SIZE  25

typedef  double MAT[SIZE][SIZE];
typedef  double VEC[SIZE];


VEC   X,Y,A,B;
MAT   T;
int   I,J,K,N;
double  AA,BB,S;


void main() {

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

  printf("\n Input the %d points, example: 1.5  0.666666\n\n", N);

  N = N - 1;

  for (K=0; K<=N; K++) {
    printf("  %d  ?  ", K); scanf("%lf %lf", &X[K], &Y[K]);
  }

// Calculate the coefficients of the linear system
  for (J=0; J<=N; J++) {
    T[0][J] = 1.0;
    for (I=1; I<=N; I++)
      T[I][J] = T[I-1][J] * X[J];
  }

// right side coefficients
  AA = 1.0; BB = 1.0;
  for (I=0; I<=N; I++) {
    AA = AA * X[0]; BB = BB * X[N];
    B[I] = (BB - AA) / (1.0*(I + 1));
  }

// Calculate coefficients A(I) and triangularize
  for (K=0; K<N; K++) {
	for (I=K+1; I<=N; I++) {
      B[I] -=  T[I][K] / T[K][K] * B[K];
      for (J=K+1; J<=N; J++)
        T[I][J] -= T[I][K] / T[K][K] * T[K][J];
    }
  }

// Solve triangular linear system T * A = B
  A[N] = B[N] / T[N][N];
  for (I=N-1; I>=0; I--) {
    S = 0.0;
    for (K=I+1; K<=N; K++)
      S += T[I][K] * A[K];
    A[I] = (B[I] - S) / T[I][I];
  }

  printf("\n Coefficients A(I):\n");
  for (I=0; I<=N; I++)
    printf(" A(%d) = %f\n", I, A[I]);

// Calculate Integral of discrete F(X)
  S = 0.0;
  for (I=0; I<=N; I++)  S += A[I] * Y[I];

// Print final result
  printf("\n Integral of F(X) from %5.2f to %5.2f\n", X[0], X[N]);
  printf("   I = %12.6f\n\n", S);

}

//end of file dinteg.cpp