1 //
2 // main.cpp
3 // 矩阵求逆
4 //
5 // Created by 唐 锐 on 13-6-20.
6 // Copyright (c) 2013年 唐 锐. All rights reserved.
7 //
8
9 #include<iostream>
10 #include<algorithm>
11 #include<iomanip>
12 #include<string>
13 #include<sstream>
14 #include<cmath>
15 #include<vector>
16 using namespace std;
17 namespace MATRIX {
18 const double eps = 1e-5;
19 template<typename T>
20 class matrix
21 {
22 vector<vector<T>>m;
23 int check()
24 {
25 if(m.empty())return 0;
26 for(int i=1;i<m.size();i++)
27 if(m[i].size()!=m[0].size())
28 return 0;
29 return 1;
30 };
31 matrix wrong()
32 {
33 matrix mat;
34 mat.m.clear();
35 vector<T>vec(1,-1);
36 mat.m.push_back(vec);
37 return mat;
38 }
39 public:
40 T at(int i,int j)
41 {
42 return m[i][j];
43 }
44 vector<T>& operator [](const int i) const
45 {
46 return m[i];
47 }
48 int size()
49 {
50 return m.size();
51 }
52 istream &input(istream& in)
53 {
54 do {
55 m.clear();
56 string s;
57 while(getline(in,s))
58 {
59 if(s=="")break;
60 istringstream temp(s);
61 T t;
62 vector<T> v ;
63 while(temp>>t)
64 v.push_back(t);
65 m.push_back(v);
66 }
67 } while(!check());
68 return in;
69 }
70 ostream &output(ostream& out) const
71 {
72 for(int i=0;i<m.size();i++)
73 {
74 for(int j=0;j<m[i].size();j++)
75 {
76 out<<m[i][j]<<' ';
77 }
78 out<<endl;
79 }
80 return out;
81 }
82 friend inline istream& operator>>(istream& in,matrix& m)
83 {
84 return m.input(in);
85 }
86 friend inline ostream& operator<<(ostream& out,const matrix& m)
87 {
88 return m.output(out);
89 }
90 void eye(int n,int k=-1)
91 {
92 if(k==-1)k=n;
93 m.clear();
94 vector<T> vec(k,0);
95 vector<vector<T>> mat(n,vec);
96 for(int i=0;i<n&&i<k;i++) mat[i][i]=1;
97 m=mat;
98 }
99
100 matrix inv()
101 {
102 if(m.size()!=m[0].size())return wrong();
103 matrix ans;
104 vector<vector<T>>cp=m;
105 ans.eye((int)m.size());
106 for(int i=0;i<cp.size();i++)
107 {
108 if(fabs(cp[i][i])<eps)
109 for(int j=i+1;j<cp.size();j++)
110 {
111 if(fabs(cp[j][i])>eps)
112 {
113 swap(cp[i],cp[j]);
114 break;
115 }
116 if(j==cp.size()-1) return wrong();
117 }
118 for(int j=i+1;j<cp.size();j++)
119 {
120 T k=cp[j][i]/cp[i][i];
121 for(int l=0;l<cp[i].size();l++)
122 {
123 cp[j][l]-=k*cp[i][l];
124 ans.m[j][l]-=k*ans.m[i][l];
125 }
126 }
127 }
128
129 for(int i=(int)cp.size()-1;i>=0;i--)
130 {
131 for(int j=i-1;j>=0;j--)
132 {
133 T k=cp[j][i]/cp[i][i];
134
135 for(int l=0;l<cp[i].size();l++)
136 {
137 cp[j][l]-=k*cp[i][l];
138 ans.m[j][l]-=k*ans.m[i][l];
139 }
140 }
141 }
142 for(int i=0;i<cp.size();i++)
143 for(int j=0;j<cp[i].size();j++)
144 {
145 ans.m[i][j]/=cp[i][i];
146 }
147 return ans;
148 }
149 vector<T>solve(vector<T>v)
150 {
151 vector<T> wrong;
152 vector<int>turn;
153 for(int i=0;i<v.size();i++)
154 turn.push_back(i);
155 if(m[0].size()!=v.size())return wrong;
156 vector<T> ans(v.size(),0);
157 vector<vector<T>>cp=m;
158 for(int i=0;i<cp.size();i++)
159 {
160 if(fabs(cp[i][i])<eps)
161 for(int j=i+1;j<cp.size();j++)
162 {
163 if(fabs(cp[j][i])>eps)
164 {
165 swap(cp[i],cp[j]);
166 swap(v[i],v[j]);
167 swap(turn[i],turn[j]);
168 break;
169 }
170 if(j==cp.size()-1) return wrong;
171 }
172 for(int j=i+1;j<cp.size();j++)
173 {
174 T k=cp[j][i]/cp[i][i];
175 for(int l=0;l<cp[i].size();l++)
176 {
177 cp[j][l]-=k*cp[i][l];
178 }
179 v[j]-=k*v[i];
180 }
181 }
182 for(int i=(int)cp.size()-1;i>=0;i--)
183 {
184 for(int j=i-1;j>=0;j--)
185 {
186 T k=cp[j][i]/cp[i][i];
187
188 for(int l=0;l<cp[i].size();l++)
189 {
190 cp[j][l]-=k*cp[i][l];
191 }
192 v[j]-=k*v[i];
193 }
194 }
195 for(int i=0;i<v.size();i++)
196 {
197 ans[turn[i]]=v[i]/cp[i][i];
198 if(fabs(ans[turn[i]])<eps)
199 ans[turn[i]]=0;
200 }
201 return ans;
202 }
203
204 };
205 template <typename T>
206 ostream &operator<<(ostream& out,const vector<T>& v)
207 {
208 for(int i=0;i<v.size();i++)
209 {
210 if(i)out<<' ';
211 out<<v[i];
212 }
213 return out;
214 }
215
216 }
217 using namespace MATRIX;
218 int main()
219 {
220 matrix<double> m;
221 double t;
222 vector<double> v;
223 while(cin>>m)
224 {
225 for(int i=0;i<m.size();i++)
226 {
227 cin>>t;
228 v.push_back(t);
229 }
230 vector<double>ans=m.solve(v);
231 cout<<ans<<endl;
232 }
233 }
234
235 /* matrix 1
236
237 1 2 0 0
238 3 4 0 0
239 0 0 4 1
240 0 0 3 2
241
242 matrix 2
243
244 1 0 -1 2 1
245 3 2 -3 5 3
246 2 2 1 4 -2
247 0 4 3 3 1
248 1 0 8 -11 4
249
250 matrix 3
251
252 1 0 -1 2 1 0 2
253 1 2 -1 3 1 -1 4
254 2 2 1 6 2 1 6
255 -1 4 1 4 0 0 0
256 4 0 -1 21 9 9 9
257 2 4 4 12 5 6 11
258 7 -1 -4 22 7 8 18
259
260 matrix 4
261
262 4 2 -3 -1 2 1 0 0 0 0
263 8 6 -5 -3 6 5 0 1 0 0
264 4 2 -2 -1 3 2 -1 0 3 1
265 0 -2 1 5 -1 3 -1 1 9 4
266 -4 2 6 -1 6 7 -3 3 2 3
267 8 6 -8 5 7 17 2 6 -3 5
268 0 2 -1 3 -4 2 5 3 0 1
269 16 10 -11 -9 17 34 2 -1 2 2
270 4 6 2 -7 13 9 2 0 12 4
271 0 0 -1 8 -3 -24 -8 6 3 -1
272
273 vector 1
274
275 5 12 3 2 3 46 13 38 19 -21
276
277 */
浙公网安备 33010602011771号