/********************************************************************
* This program solves a Vandermonde linear system:                  *
*               N      k-1                                          *
*              S     xi    wi = qk  (k=1,...,N)                     *
*               i=1                                                 *
* ----------------------------------------------------------------- *
* SAMPLE RUN:                                                       *
*   N=10                                                            *
*   X(I)=I (I=1,10)                                                 *
*   Q(I)=1, except Q(9)=30                                          *
*                                                                   *
* The solution vector is:                                           *
*    1    1.004315                                                  *
*    2    -0.038120                                                 *
*    3    0.149603                                                  *
*    4    -0.342631                                                 *
*    5    0.503472                                                  *
*    6    -0.493403                                                 *
*    7    0.322222                                                  *
*    8    -0.135218                                                 *
*    9    0.033085                                                  *
*   10    -0.003596                                                 *
*                                                                   *
* ----------------------------------------------------------------- *
* Reference: "Numerical Recipes by W.H. Press, B.P. Flannery, S.A.  *
*             Teukolsky, W.T. Vetterling, Cambridge University      *
*             Press, 1987"                                          *
*                                                                   *
*                         Visual C++ Version By J-P Moreau, Paris   *
*                                   (www.jpmoreau.fr)               *
********************************************************************/
#include <stdio.h>
#include <math.h>

int i,N;
double X[11], Q[11], W[11];   // index zero not used here


void Vander(double *X, double *W, double *Q, int N)  {
/*--------------------------------------------------------------------
!                                        N     k-1
! Solves the Vandermonde linear system S     xi    wi = qk (k=1,...,N)
!                                        i=1
! Input consists of the vectors X and Q, each of length N, the vector
! W is output.
!--------------------------------------------------------------------*/
const int NMAX=100; double ZERO=0.0, ONE=1.0;
// NMAX is the maximum expected value of N.
double B, S, T, XX;
int i,j,k,k1;
double C[11];

if (N==1)
  W[1]=Q[1];
else {
  for (i=1;i<=N; i++) // initialize array C
    C[i]=ZERO;
  C[N]=-X[1];         // coefficients of the master polynomial
                      // are found by recursion.
  for (i=2; i<=N; i++) {
    XX=-X[i];
    for (j=N+1-i; j<N; j++)
      C[j] += XX*C[j+1];
   C[N]=C[N]+XX;
  }
  for (i=1; i<=N; i++) {   // each subfactor in turn
    XX=X[i];
    T=ONE;
    B=ONE;
    S=Q[N];
    k=N;
    for (j=2; j<=N; j++) { // is synthetically divided,
      k1=k-1;
      B=C[k]+XX*B;
      S=S+Q[k1]*B;         // matrix-multiplied by the right-hand side,
      T=XX*T+B;
      k=k1;
    }
    W[i]=S/T;              // and suppliedwith a denominator.
  }
}

}


void main() {

  N=10;

  for (i=1; i<=N; i++)  X[i]=1.0*i;
  for (i=1; i<=N; i++)  Q[i]=1.0;
  Q[9]=30.0;

  Vander(X,W,Q,N);

  printf("\n The solution vector is:\n");
  for (i=1; i<=N; i++)
    printf(" %4d    %f\n", i, W[i]);

}

// end of file tvander.cpp