#include "iostream.h"
#include "math.h"

void toeplz(double r[], double x[], double y[], int n)
{
	int k,j,m,m1,m2;
	double sxn,sd,sgn1,shn,sgd,pt1,pt2,qt1,qt2,pp,qq;
    const int nmax = 100;
    double g[100], h[100];
    if (r[n] == 0.0) 
	{
		cout<<"levinson method fails: singular principal minor"<<endl;
		return;
	}
	x[1] = y[1] / r[n];
    if (n == 1)
	{
		return;
	}
    g[1] = r[n - 1] / r[n];
    h[1] = r[n + 1] / r[n];
    for (m = 1; m<=n; m++)
	{
        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 == 0.0)
		{
			cout<<"levinson method fails: singular principal minor"<<endl;
			return;
		}
        x[m1] = sxn / sd;
        for (j = 1; j<=m; j++)
		{
            x[j] = x[j] - x[m1] * g[m - j + 1];
        }
        if (m1 == n)
		{
			return;
		}
        sgn1 = -r[n - m1];
        shn = -r[n + m1];
        sgd = -r[n];
        for (j = 1; j<=m; j++)
		{
            sgn1 = sgn1 + r[n + j - m1] * g[j];
            shn = shn + r[n + m1 - j] * h[j];
            sgd = sgd + r[n + j - m1] * h[m - j + 1];
        }
        if ((sd == 0.0) || (sgd == 0.0))
		{
			cout<<"levinson method fails: singular principal minor"<<endl;
			return;
		}
        g[m1] = sgn1 / sgd;
        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 = k - 1;
        }
    }
    cout<<"never get here"<<endl;
}

void main()
{
    //program d1r7
    //driver program for routine toeplz
    int i,j,n2,n = 5;
	double sum1;
    n2 = 2 * n;
    double x[6], y[6], r[11];
    //輸入已知的方程組的系數矩陣的r
    for (i = 1; i<=2 * n - 1; i++)
	{
		r[i] = 1 / double(i);
    }
    //輸入已知的方程組的右端向量
    for (i = 1; i<=n; i++)
	{
        y[i] = 0.1 * i;
    }
	cout.setf(ios::fixed|ios::right);
	cout.precision(5);
    cout<<"已知的方程組的右端向量"<<endl;
	cout.width(12);    cout<<y[1]<<endl;
	cout.width(12);    cout<<y[2]<<endl;
	cout.width(12);    cout<<y[3]<<endl;
	cout.width(12);    cout<<y[4]<<endl;
	cout.width(12);    cout<<y[5]<<endl;
    toeplz(r, x, y, n);
    //輸出方程組的解x
    cout<<endl;
    cout<<"計算出的方程組的的解"<<endl;
	cout.width(15);    cout<<x[1]<<endl;
	cout.width(15);    cout<<x[2]<<endl;
	cout.width(15);    cout<<x[3]<<endl;
	cout.width(15);    cout<<x[4]<<endl;
	cout.width(15);    cout<<x[5]<<endl;
    //將計算出的解x乘以系數矩陣，以驗證計算結果正確
    cout<<endl;
    cout<<"將計算出的解乘以系數矩陣，以驗證計算結果正確"<<endl;
    cout<<endl;
    cout<<"解乘以系數矩陣"<<"        "<<"方程組的右端向量"<<endl;
    for (i = 1; i<=n; i++)
	{
        sum1 = 0;
        for (j = 1; j<=n; j++)
		{
           sum1 = sum1 + r[n + i - j] * x[j];
        }
        cout<<"    "<<sum1<<"              "<<y[i]<<endl;
    }
}