//-------------------------------------------------------------
// Linear Equation Solver for Teplitz Matrix
//   O(n^2)
//
// Faster solver for matrices satisfying
//   M[i][j] == M[i+1][j+1]
// for all i, j>0
// Be careful that this implementation does not support the
// case M[0][0]=M[1][1]=...=0. (zero-division)
//
// TODO: support M[0][0]=M[1][1]=...=0
// QUICKHACK:
//    find the least degenerated k and 
//    calc initial vectors by gauss_jordan...
//
// Verified by
//   - (Checked by random data to match with linearEq.cpp)
//-------------------------------------------------------------

vector<double> solve_teplitz( int n, vector< vector<double> > M, const vector<double>& V )
{
	vector<double> f, b, x;

	// invariants
	//   M|k f|k = (1 0 ... 0)
	//   M|k b|k = (0 ... 0 1)
	//   M|k x|k = V|k
	f.push_back(    1 / M[0][0] );
	b.push_back(    1 / M[0][0] );
	x.push_back( V[0] / M[0][0] );

	for(int k=2; k<=n; ++k)
	{
		// M|k (f|k-1 0) = (1 0 ... 0 ea)
		double ea = 0;
		for(int i=0; i<k-1; ++i)
			ea += M[k-1][i] * f[i];

		// M|k (0 b|k-1) = (eb 0 ... 0 1)
		double eb = 0;
		for(int i=1; i<k; ++i)
			eb += M[0][i] * b[i-1];

		// u(1 ea) + v(eb 1) = (1 0)
		//==> u = 1 / (1-ea.eb),  v = -ea /(1-ea.eb)
		double u = 1/(1-ea*eb), v = -ea/(1-ea*eb);
		vector<double> ff;
		for(int i=0; i<k; ++i)
			ff.push_back( u*(i<k-1 ? f[i] : 0) + v*(i>0 ? b[i-1] : 0) );

		// similarly...
		u = -eb/(1-ea*eb), v = 1/(1-ea*eb);
		vector<double> bb;
		for(int i=0; i<k; ++i)
			bb.push_back( u*(i<k-1 ? f[i] : 0) + v*(i>0 ? b[i-1] : 0) );

		// update
		f.swap(ff);
		b.swap(bb);

		// M|k (x|k 0) (V|k ec)
		double ec = 0;
		for(int i=0; i<k-1; ++i)
			ec += M[k-1][i] * x[i];
		x.push_back(0);
		for(int i=0; i<k; ++i)
			x[i] -= b[i] * (ec-V[k-1]);
	}

	return x;
}
