/********************************************************************
* This program solves a Toeplitz linear system:                     *
*               N                                                   *
*              S     R      x  = y   (i=1,...,N)                    *
*               j=1   N+i-j  j    i                                 *
* ----------------------------------------------------------------- *
* SAMPLE RUN:                                                       *
*   N=10                                                            *
*   R = (-10,-12,3,4,5,6,7,8,9,5,11,12,13,14,15,16,17,18,19)        *
*   Y = (1,1,1,1,1,1,1,1,30,1)                                      *
*                                                                   *
* The solution vector is:                                           *
*  1   4.145258                                                     *
*  2   -7.202165                                                    *
*  3   0.136163                                                     *
*  4   0.473885                                                     *
*  5   0.811607                                                     *
*  6   1.149329                                                     *
*  7   1.487051                                                     *
*  8   1.824773                                                     *
*  9   -3.637505                                                    *
* 10   2.500217                                                     *
*                                                                   *
* ----------------------------------------------------------------- *
* Reference: "Numerical Recipes by W.H. Press, B.P. Flannery, S.A.  *
*             Teukolsky, W.T. Vetterling, Cambridge University      *
*             Press, 1987"                                          *
*                                                                   *
*                             C++ Version By J-P Moreau, Paris      *
*                                    (www.jpmoreau.fr)              *
********************************************************************/
#include <stdio.h>


void TOEPLZ(double *R, double *X, double *Y, int N)  {
/*--------------------------------------------------------------------
 * This program solves a Toeplitz linear system:                     *
 *               N                                                   *
 *              S     R      x  = y   (i=1,...,N)                    *
 *               j=1   N+i-j  j    i                                 *
 *                                                                   *
 *  Input consists of vectors R of size 2N-1 and Y of size N.        *
 *  Output is solution vector X of size N.                           *
 *  N is the size of the linear system.                              *
 --------------------------------------------------------------------*/
const int NMAX=100; double ZERO=0.0;
// NMAX is the maximum expected value of N.
double G[NMAX], H[NMAX];
double PP, PT1, PT2, QQ, QT1, QT2, SD, SGD, SGN, SHN, SXN;
int J, K, M, M1, M2;

if (R[N] == ZERO) goto e99;
X[1]=Y[1]/R[N];             // initialize for the recursion.
if(N == 1) return;
G[1]=R[N-1]/R[N];
H[1]=R[N+1]/R[N];

for (M=1; M<=N; M++) {     // main loop over the recursion. 
  M1=M+1;
  SXN=-Y[M1];
  SD=-R[N];
  for (J=1; J<=M; J++) {
	SXN=SXN+R[N+M1-J]*X[J];
	SD=SD + R[N+M1-J]*G[M-J+1];
  }
  if (SD == ZERO) goto e99;
  X[M1]=SXN/SD;             // whence x.
  for (J=1; J<=M; J++)
    X[J] -= X[M1]*G[M-J+1];
  if (M1 == N)  return;     // normal exit; 
  SGN=-R[N-M1];             // compute numerator and denominator for G and H.
  SHN=-R[N+M1];
  SGD=-R[N];
  for (J=1; J<=M; J++) {
    SGN += R[N+J-M1]*G[J];
    SHN += R[N-J+M1]*H[J];
	SGD += R[N+J-M1]*H[M-J+1];
  }

  if(SD == ZERO || SGD == ZERO) goto e99;
  G[M1]=SGN/SGD;            // whence G and H
  H[M1]=SHN/SD;
  K=M;
  M2=(M+1)/2;
  PP=G[M1];
  QQ=H[M1];
  for (J=1; J<=M2; J++) {
    PT1=G[J];
	PT2=G[K];
    QT1=H[J];
	QT2=H[K];
    G[J]=PT1-PP*QT2;
    G[K]=PT2-PP*QT1;
	H[J]=QT1-QQ*PT2;
    H[K]=QT2-QQ*PT1;
	K--;
  }
} // M loop  	 

e99:printf(" Levinson method fails: singular principal minor.\n");

} // TOEPLZ


// main program
void main() {

  int i, N;

  double R[20] = { 0.0,-10.0, -12.0,  3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 5.0,
	  11.0,  12.0, 13.0,14.0,15.0,16.0,17.0,18.0,19.0};
  double Y[11] = { 0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,30.0, 1.0};

  double X[11];

  N = 10;
  
  TOEPLZ(R,X,Y,N);
  
  printf("\n The solution vector is:\n");
  for (i=1; i<=N; i++)
    printf(" %2d   %f\n", i, X[i]);

}

// end of file toeplitz.cpp