Artifact bf7f1d164fa346363d8bc257b4591b1250f97cc5
//-------------------------------------------------------------
// 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;
}