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 */