Gauss-Jordanの消去法。蟻本のパクリ。
/* 実行方法: ./a.out < input.txt input.txtの中身: 3 1 -2 3 4 -5 6 7 -8 10 6 12 21 先頭行は元数(n)。 続くn行は左辺の係数。 最終行は右辺の定数。 上の例だと x - 2y + 3z = 6 4x - 5y + 6z = 12 7x - 8y + 10z = 21 を意味する。 これを解くと (x, y, z) = (1, 2, 3) となるので、 1, 2, 3, と出力する。 */ #include <iostream> #include <vector> #include <cmath> using namespace std; const double EPS = 1e-9; vector<double> solve(const vector<vector<double> >& A, const vector<double>& b) { int n = A.size(); vector<vector<double> > B(n, vector<double>(n+1)); for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { B[i][j] = A[i][j]; } } for(int i=0; i<n; i++) { B[i][n] = b[i]; } for(int i=0; i<n; i++) { int pivot = i; for(int j=i; j<n; j++) { if(abs(B[j][i]) > abs(B[pivot][i])) { pivot = j; } } swap(B[i], B[pivot]); // 解がないか、一意でない場合は長さ0の配列を返すことにする。 if(abs(B[i][i]) < EPS) { return vector<double>(); } for(int j=i+1; j<=n; j++) { B[i][j] /= B[i][i]; } for(int j=0; j<n; j++) { if(i!=j) { for(int k=i+1; k<=n; k++) { B[j][k] -= B[j][i] * B[i][k]; } } } } vector<double> res(n); for(int i=0; i<n; i++) { res[i] = B[i][n]; } return res; } int main(void) { int n; cin >> n; vector<vector<double> > A(n, vector<double>(n+1)); vector<double> b(n); for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { cin >> A[i][j]; } } for(int i=0; i<n; i++) { cin >> b[i]; } vector<double> res = solve(A, b); int sz = res.size(); for(int i=0; i<sz; i++) { cout << res[i] << ", "; } cout << endl; return 0; }