Artifact bf7f1d164fa346363d8bc257b4591b1250f97cc5:
0000: 2f 2f 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d //--------------
0010: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d ----------------
0020: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d ----------------
0030: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 0a ---------------.
0040: 2f 2f 20 4c 69 6e 65 61 72 20 45 71 75 61 74 69 // Linear Equati
0050: 6f 6e 20 53 6f 6c 76 65 72 20 66 6f 72 20 54 65 on Solver for Te
0060: 70 6c 69 74 7a 20 4d 61 74 72 69 78 0a 2f 2f 20 plitz Matrix.//
0070: 20 20 4f 28 6e 5e 32 29 0a 2f 2f 0a 2f 2f 20 46 O(n^2).//.// F
0080: 61 73 74 65 72 20 73 6f 6c 76 65 72 20 66 6f 72 aster solver for
0090: 20 6d 61 74 72 69 63 65 73 20 73 61 74 69 73 66 matrices satisf
00a0: 79 69 6e 67 0a 2f 2f 20 20 20 4d 5b 69 5d 5b 6a ying.// M[i][j
00b0: 5d 20 3d 3d 20 4d 5b 69 2b 31 5d 5b 6a 2b 31 5d ] == M[i+1][j+1]
00c0: 0a 2f 2f 20 66 6f 72 20 61 6c 6c 20 69 2c 20 6a .// for all i, j
00d0: 3e 30 0a 2f 2f 20 42 65 20 63 61 72 65 66 75 6c >0.// Be careful
00e0: 20 74 68 61 74 20 74 68 69 73 20 69 6d 70 6c 65 that this imple
00f0: 6d 65 6e 74 61 74 69 6f 6e 20 64 6f 65 73 20 6e mentation does n
0100: 6f 74 20 73 75 70 70 6f 72 74 20 74 68 65 0a 2f ot support the./
0110: 2f 20 63 61 73 65 20 4d 5b 30 5d 5b 30 5d 3d 4d / case M[0][0]=M
0120: 5b 31 5d 5b 31 5d 3d 2e 2e 2e 3d 30 2e 20 28 7a [1][1]=...=0. (z
0130: 65 72 6f 2d 64 69 76 69 73 69 6f 6e 29 0a 2f 2f ero-division).//
0140: 0a 2f 2f 20 54 4f 44 4f 3a 20 73 75 70 70 6f 72 .// TODO: suppor
0150: 74 20 4d 5b 30 5d 5b 30 5d 3d 4d 5b 31 5d 5b 31 t M[0][0]=M[1][1
0160: 5d 3d 2e 2e 2e 3d 30 0a 2f 2f 20 51 55 49 43 4b ]=...=0.// QUICK
0170: 48 41 43 4b 3a 0a 2f 2f 20 20 20 20 66 69 6e 64 HACK:.// find
0180: 20 74 68 65 20 6c 65 61 73 74 20 64 65 67 65 6e the least degen
0190: 65 72 61 74 65 64 20 6b 20 61 6e 64 20 0a 2f 2f erated k and .//
01a0: 20 20 20 20 63 61 6c 63 20 69 6e 69 74 69 61 6c calc initial
01b0: 20 76 65 63 74 6f 72 73 20 62 79 20 67 61 75 73 vectors by gaus
01c0: 73 5f 6a 6f 72 64 61 6e 2e 2e 2e 0a 2f 2f 0a 2f s_jordan....//./
01d0: 2f 20 56 65 72 69 66 69 65 64 20 62 79 0a 2f 2f / Verified by.//
01e0: 20 20 20 2d 20 28 43 68 65 63 6b 65 64 20 62 79 - (Checked by
01f0: 20 72 61 6e 64 6f 6d 20 64 61 74 61 20 74 6f 20 random data to
0200: 6d 61 74 63 68 20 77 69 74 68 20 6c 69 6e 65 61 match with linea
0210: 72 45 71 2e 63 70 70 29 0a 2f 2f 2d 2d 2d 2d 2d rEq.cpp).//-----
0220: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d ----------------
0230: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d ----------------
0240: 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d 2d ----------------
0250: 2d 2d 2d 2d 2d 2d 2d 2d 0a 0a 76 65 63 74 6f 72 --------..vector
0260: 3c 64 6f 75 62 6c 65 3e 20 73 6f 6c 76 65 5f 74 <double> solve_t
0270: 65 70 6c 69 74 7a 28 20 69 6e 74 20 6e 2c 20 76 eplitz( int n, v
0280: 65 63 74 6f 72 3c 20 76 65 63 74 6f 72 3c 64 6f ector< vector<do
0290: 75 62 6c 65 3e 20 3e 20 4d 2c 20 63 6f 6e 73 74 uble> > M, const
02a0: 20 76 65 63 74 6f 72 3c 64 6f 75 62 6c 65 3e 26 vector<double>&
02b0: 20 56 20 29 0a 7b 0a 09 76 65 63 74 6f 72 3c 64 V ).{..vector<d
02c0: 6f 75 62 6c 65 3e 20 66 2c 20 62 2c 20 78 3b 0a ouble> f, b, x;.
02d0: 0a 09 2f 2f 20 69 6e 76 61 72 69 61 6e 74 73 0a ..// invariants.
02e0: 09 2f 2f 20 20 20 4d 7c 6b 20 66 7c 6b 20 3d 20 .// M|k f|k =
02f0: 28 31 20 30 20 2e 2e 2e 20 30 29 0a 09 2f 2f 20 (1 0 ... 0)..//
0300: 20 20 4d 7c 6b 20 62 7c 6b 20 3d 20 28 30 20 2e M|k b|k = (0 .
0310: 2e 2e 20 30 20 31 29 0a 09 2f 2f 20 20 20 4d 7c .. 0 1)..// M|
0320: 6b 20 78 7c 6b 20 3d 20 56 7c 6b 0a 09 66 2e 70 k x|k = V|k..f.p
0330: 75 73 68 5f 62 61 63 6b 28 20 20 20 20 31 20 2f ush_back( 1 /
0340: 20 4d 5b 30 5d 5b 30 5d 20 29 3b 0a 09 62 2e 70 M[0][0] );..b.p
0350: 75 73 68 5f 62 61 63 6b 28 20 20 20 20 31 20 2f ush_back( 1 /
0360: 20 4d 5b 30 5d 5b 30 5d 20 29 3b 0a 09 78 2e 70 M[0][0] );..x.p
0370: 75 73 68 5f 62 61 63 6b 28 20 56 5b 30 5d 20 2f ush_back( V[0] /
0380: 20 4d 5b 30 5d 5b 30 5d 20 29 3b 0a 0a 09 66 6f M[0][0] );...fo
0390: 72 28 69 6e 74 20 6b 3d 32 3b 20 6b 3c 3d 6e 3b r(int k=2; k<=n;
03a0: 20 2b 2b 6b 29 0a 09 7b 0a 09 09 2f 2f 20 4d 7c ++k)..{...// M|
03b0: 6b 20 28 66 7c 6b 2d 31 20 30 29 20 3d 20 28 31 k (f|k-1 0) = (1
03c0: 20 30 20 2e 2e 2e 20 30 20 65 61 29 0a 09 09 64 0 ... 0 ea)...d
03d0: 6f 75 62 6c 65 20 65 61 20 3d 20 30 3b 0a 09 09 ouble ea = 0;...
03e0: 66 6f 72 28 69 6e 74 20 69 3d 30 3b 20 69 3c 6b for(int i=0; i<k
03f0: 2d 31 3b 20 2b 2b 69 29 0a 09 09 09 65 61 20 2b -1; ++i)....ea +
0400: 3d 20 4d 5b 6b 2d 31 5d 5b 69 5d 20 2a 20 66 5b = M[k-1][i] * f[
0410: 69 5d 3b 0a 0a 09 09 2f 2f 20 4d 7c 6b 20 28 30 i];....// M|k (0
0420: 20 62 7c 6b 2d 31 29 20 3d 20 28 65 62 20 30 20 b|k-1) = (eb 0
0430: 2e 2e 2e 20 30 20 31 29 0a 09 09 64 6f 75 62 6c ... 0 1)...doubl
0440: 65 20 65 62 20 3d 20 30 3b 0a 09 09 66 6f 72 28 e eb = 0;...for(
0450: 69 6e 74 20 69 3d 31 3b 20 69 3c 6b 3b 20 2b 2b int i=1; i<k; ++
0460: 69 29 0a 09 09 09 65 62 20 2b 3d 20 4d 5b 30 5d i)....eb += M[0]
0470: 5b 69 5d 20 2a 20 62 5b 69 2d 31 5d 3b 0a 0a 09 [i] * b[i-1];...
0480: 09 2f 2f 20 75 28 31 20 65 61 29 20 2b 20 76 28 .// u(1 ea) + v(
0490: 65 62 20 31 29 20 3d 20 28 31 20 30 29 0a 09 09 eb 1) = (1 0)...
04a0: 2f 2f 3d 3d 3e 20 75 20 3d 20 31 20 2f 20 28 31 //==> u = 1 / (1
04b0: 2d 65 61 2e 65 62 29 2c 20 20 76 20 3d 20 2d 65 -ea.eb), v = -e
04c0: 61 20 2f 28 31 2d 65 61 2e 65 62 29 0a 09 09 64 a /(1-ea.eb)...d
04d0: 6f 75 62 6c 65 20 75 20 3d 20 31 2f 28 31 2d 65 ouble u = 1/(1-e
04e0: 61 2a 65 62 29 2c 20 76 20 3d 20 2d 65 61 2f 28 a*eb), v = -ea/(
04f0: 31 2d 65 61 2a 65 62 29 3b 0a 09 09 76 65 63 74 1-ea*eb);...vect
0500: 6f 72 3c 64 6f 75 62 6c 65 3e 20 66 66 3b 0a 09 or<double> ff;..
0510: 09 66 6f 72 28 69 6e 74 20 69 3d 30 3b 20 69 3c .for(int i=0; i<
0520: 6b 3b 20 2b 2b 69 29 0a 09 09 09 66 66 2e 70 75 k; ++i)....ff.pu
0530: 73 68 5f 62 61 63 6b 28 20 75 2a 28 69 3c 6b 2d sh_back( u*(i<k-
0540: 31 20 3f 20 66 5b 69 5d 20 3a 20 30 29 20 2b 20 1 ? f[i] : 0) +
0550: 76 2a 28 69 3e 30 20 3f 20 62 5b 69 2d 31 5d 20 v*(i>0 ? b[i-1]
0560: 3a 20 30 29 20 29 3b 0a 0a 09 09 2f 2f 20 73 69 : 0) );....// si
0570: 6d 69 6c 61 72 6c 79 2e 2e 2e 0a 09 09 75 20 3d milarly......u =
0580: 20 2d 65 62 2f 28 31 2d 65 61 2a 65 62 29 2c 20 -eb/(1-ea*eb),
0590: 76 20 3d 20 31 2f 28 31 2d 65 61 2a 65 62 29 3b v = 1/(1-ea*eb);
05a0: 0a 09 09 76 65 63 74 6f 72 3c 64 6f 75 62 6c 65 ...vector<double
05b0: 3e 20 62 62 3b 0a 09 09 66 6f 72 28 69 6e 74 20 > bb;...for(int
05c0: 69 3d 30 3b 20 69 3c 6b 3b 20 2b 2b 69 29 0a 09 i=0; i<k; ++i)..
05d0: 09 09 62 62 2e 70 75 73 68 5f 62 61 63 6b 28 20 ..bb.push_back(
05e0: 75 2a 28 69 3c 6b 2d 31 20 3f 20 66 5b 69 5d 20 u*(i<k-1 ? f[i]
05f0: 3a 20 30 29 20 2b 20 76 2a 28 69 3e 30 20 3f 20 : 0) + v*(i>0 ?
0600: 62 5b 69 2d 31 5d 20 3a 20 30 29 20 29 3b 0a 0a b[i-1] : 0) );..
0610: 09 09 2f 2f 20 75 70 64 61 74 65 0a 09 09 66 2e ..// update...f.
0620: 73 77 61 70 28 66 66 29 3b 0a 09 09 62 2e 73 77 swap(ff);...b.sw
0630: 61 70 28 62 62 29 3b 0a 0a 09 09 2f 2f 20 4d 7c ap(bb);....// M|
0640: 6b 20 28 78 7c 6b 20 30 29 20 28 56 7c 6b 20 65 k (x|k 0) (V|k e
0650: 63 29 0a 09 09 64 6f 75 62 6c 65 20 65 63 20 3d c)...double ec =
0660: 20 30 3b 0a 09 09 66 6f 72 28 69 6e 74 20 69 3d 0;...for(int i=
0670: 30 3b 20 69 3c 6b 2d 31 3b 20 2b 2b 69 29 0a 09 0; i<k-1; ++i)..
0680: 09 09 65 63 20 2b 3d 20 4d 5b 6b 2d 31 5d 5b 69 ..ec += M[k-1][i
0690: 5d 20 2a 20 78 5b 69 5d 3b 0a 09 09 78 2e 70 75 ] * x[i];...x.pu
06a0: 73 68 5f 62 61 63 6b 28 30 29 3b 0a 09 09 66 6f sh_back(0);...fo
06b0: 72 28 69 6e 74 20 69 3d 30 3b 20 69 3c 6b 3b 20 r(int i=0; i<k;
06c0: 2b 2b 69 29 0a 09 09 09 78 5b 69 5d 20 2d 3d 20 ++i)....x[i] -=
06d0: 62 5b 69 5d 20 2a 20 28 65 63 2d 56 5b 6b 2d 31 b[i] * (ec-V[k-1
06e0: 5d 29 3b 0a 09 7d 0a 0a 09 72 65 74 75 72 6e 20 ]);..}...return
06f0: 78 3b 0a 7d 0a x;.}.