matrix矩阵求逆 与解方程模板 留做备用 (有bug,待补充)

//
//  main.cpp
//  矩阵求逆
//
//  Created by 唐 锐 on 13-6-20.
//  Copyright (c) 2013年 唐 锐. All rights reserved.
//

#include<iostream>
#include<algorithm>
#include<iomanip>
#include<string>
#include<sstream>
#include<cmath>
#include<vector>
using namespace std;
namespace MATRIX {
    const double eps = 1e-5;
    template<typename T>
    class matrix
    {
        vector<vector<T>>m;
        int check()
        {
            if(m.empty())return 0;
            for(int i=1;i<m.size();i++)
                if(m[i].size()!=m[0].size())
                    return 0;
            return 1;
        };
        matrix wrong()
        {
            matrix mat;
            mat.m.clear();
            vector<T>vec(1,-1);
            mat.m.push_back(vec);
            return mat;
        }
    public:
        T at(int i,int j)
        {
            return m[i][j];
        }
        vector<T>& operator [](const int i) const
        {
            return m[i];
        }
        int size()
        {
            return m.size();
        }
        istream &input(istream& in)
        {
            do {
                m.clear();
                string s;
                while(getline(in,s))
                {
                    if(s=="")break;
                    istringstream temp(s);
                    T t;
                    vector<T> v ;
                    while(temp>>t)
                        v.push_back(t);
                    m.push_back(v);
                }
            } while(!check());
            return in;
        }
        ostream &output(ostream& out) const
        {
            for(int i=0;i<m.size();i++)
            {
                for(int j=0;j<m[i].size();j++)
                {
                    out<<m[i][j]<<' ';
                }
                out<<endl;
            }
            return out;
        }
        friend inline istream& operator>>(istream& in,matrix& m)
        {
            return m.input(in);
        }
        friend inline ostream& operator<<(ostream& out,const matrix& m)
        {
            return m.output(out);
        }
        void eye(int n,int k=-1)
        {
            if(k==-1)k=n;
            m.clear();
            vector<T> vec(k,0);
            vector<vector<T>> mat(n,vec);
            for(int i=0;i<n&&i<k;i++) mat[i][i]=1;
            m=mat;
        }
        
        matrix inv()
        {
            if(m.size()!=m[0].size())return wrong();
            matrix ans;
            vector<vector<T>>cp=m;
            ans.eye((int)m.size());
            for(int i=0;i<cp.size();i++)
            {
                if(fabs(cp[i][i])<eps)
                    for(int j=i+1;j<cp.size();j++)
                    {
                        if(fabs(cp[j][i])>eps)
                        {
                            swap(cp[i],cp[j]);
                            break;
                        }
                        if(j==cp.size()-1) return wrong();
                    }
                for(int j=i+1;j<cp.size();j++)
                {
                    T k=cp[j][i]/cp[i][i];
                    for(int l=0;l<cp[i].size();l++)
                    {
                        cp[j][l]-=k*cp[i][l];
                        ans.m[j][l]-=k*ans.m[i][l];
                    }
                }
            }
            
            for(int i=(int)cp.size()-1;i>=0;i--)
            {
                for(int j=i-1;j>=0;j--)
                {
                    T k=cp[j][i]/cp[i][i];
                    
                    for(int l=0;l<cp[i].size();l++)
                    {
                        cp[j][l]-=k*cp[i][l];
                        ans.m[j][l]-=k*ans.m[i][l];
                    }
                }
            }
            for(int i=0;i<cp.size();i++)
                for(int j=0;j<cp[i].size();j++)
                {
                    ans.m[i][j]/=cp[i][i];
                }
            return ans;
        }
        vector<T>solve(vector<T>v)
        {
            vector<T> wrong;
            vector<int>turn;
            for(int i=0;i<v.size();i++)
                turn.push_back(i);
            if(m[0].size()!=v.size())return wrong;
            vector<T> ans(v.size(),0);
            vector<vector<T>>cp=m;
            for(int i=0;i<cp.size();i++)
            {
                if(fabs(cp[i][i])<eps)
                    for(int j=i+1;j<cp.size();j++)
                    {
                        if(fabs(cp[j][i])>eps)
                        {
                            swap(cp[i],cp[j]);
                            swap(v[i],v[j]);
                            swap(turn[i],turn[j]);
                            break;
                        }
                        if(j==cp.size()-1) return wrong;
                    }
                for(int j=i+1;j<cp.size();j++)
                {
                    T k=cp[j][i]/cp[i][i];
                    for(int l=0;l<cp[i].size();l++)
                    {
                        cp[j][l]-=k*cp[i][l];
                    }
                    v[j]-=k*v[i];
                }
            }
            for(int i=(int)cp.size()-1;i>=0;i--)
            {
                for(int j=i-1;j>=0;j--)
                {
                    T k=cp[j][i]/cp[i][i];
                    
                    for(int l=0;l<cp[i].size();l++)
                    {
                        cp[j][l]-=k*cp[i][l];
                    }
                    v[j]-=k*v[i];
                }
            }
            for(int i=0;i<v.size();i++)
            {
                ans[turn[i]]=v[i]/cp[i][i];
                if(fabs(ans[turn[i]])<eps)
                    ans[turn[i]]=0;
            }
            return ans;
        }
        
    };
    template <typename T>
    ostream &operator<<(ostream& out,const vector<T>& v)
    {
        for(int i=0;i<v.size();i++)
        {
            if(i)out<<' ';
            out<<v[i];
        }
        return out;
    }
    
}
using namespace MATRIX;
int main()
{
    matrix<double> m;
    double t;
    vector<double> v;
    while(cin>>m)
    {
        for(int i=0;i<m.size();i++)
        {
            cin>>t;
            v.push_back(t);
        }
        vector<double>ans=m.solve(v);
        cout<<ans<<endl;
    }
}

/* matrix 1
 
 1 2 0 0
 3 4 0 0
 0 0 4 1
 0 0 3 2
 
 matrix 2
 
 1 0 -1 2 1
 3 2 -3 5 3
 2 2 1 4 -2
 0 4 3 3 1
 1 0 8 -11 4
 
 matrix 3
 
 1 0 -1 2 1 0 2
 1 2 -1 3 1 -1 4
 2 2 1 6 2 1 6
 -1 4 1 4 0 0 0
 4 0 -1 21 9 9 9
 2 4 4 12 5 6 11
 7 -1 -4 22 7 8 18
 
 matrix 4
 
 4 2 -3 -1 2 1 0 0 0 0
 8 6 -5 -3 6 5 0 1 0 0
 4 2 -2 -1 3 2 -1 0 3 1
 0 -2 1 5 -1 3 -1 1 9 4
 -4 2 6 -1 6 7 -3 3 2 3
 8 6 -8 5 7 17 2 6 -3 5
 0 2 -1 3 -4 2 5 3 0 1
 16 10 -11 -9 17 34 2 -1 2 2
 4 6 2 -7 13 9 2 0 12 4
 0 0 -1 8 -3 -24 -8 6 3 -1
 
 vector 1
 
 5 12 3 2 3 46 13 38 19 -21
 
 */