23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: //-------------------------------------------------------------
23dfcca431 2011-02-23        kinaba: // Linear Equation Solver
23dfcca431 2011-02-23        kinaba: //   O(n^3)
23dfcca431 2011-02-23        kinaba: //
23dfcca431 2011-02-23        kinaba: // Verified by
23dfcca431 2011-02-23        kinaba: //   - SRM 398 Div1 LV3
23dfcca431 2011-02-23        kinaba: //   - AOJ 0004
23dfcca431 2011-02-23        kinaba: //-------------------------------------------------------------
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: vector<double> solve_linear_eq( int n, vector< vector<double> > M, const vector<double>& V )
23dfcca431 2011-02-23        kinaba: {
23dfcca431 2011-02-23        kinaba: 	vector<double> A(V);
23dfcca431 2011-02-23        kinaba: 	for(int i=0; i<n; ++i)
23dfcca431 2011-02-23        kinaba: 	{
23dfcca431 2011-02-23        kinaba: 		// pivot
23dfcca431 2011-02-23        kinaba: 		if( M[i][i] == 0 )
23dfcca431 2011-02-23        kinaba: 			for(int j=i+1; j<n; ++j)
23dfcca431 2011-02-23        kinaba: 				if( M[j][i] != 0 )
23dfcca431 2011-02-23        kinaba: 					{swap(M[i], M[j]); swap(A[i], A[j]); break;}
23dfcca431 2011-02-23        kinaba: 		if( M[i][i] == 0 )
23dfcca431 2011-02-23        kinaba: 			throw "no anser";
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: 		// M[i][i] <-- 1
23dfcca431 2011-02-23        kinaba: 		double p = M[i][i];
23dfcca431 2011-02-23        kinaba: 		for(int j=i; j<n; ++j)
23dfcca431 2011-02-23        kinaba: 			M[i][j] /= p;
23dfcca431 2011-02-23        kinaba: 		A[i] /= p;
23dfcca431 2011-02-23        kinaba: 
23dfcca431 2011-02-23        kinaba: 		// M[*][i] <-- 0
23dfcca431 2011-02-23        kinaba: 		for(int j=0; j<n; ++j) if(j!=i)
23dfcca431 2011-02-23        kinaba: 		{
23dfcca431 2011-02-23        kinaba: 			double r = M[j][i];
23dfcca431 2011-02-23        kinaba: 			for(int k=i; k<n; ++k)
23dfcca431 2011-02-23        kinaba: 				M[j][k] -= M[i][k] * r;
23dfcca431 2011-02-23        kinaba: 			A[j] -= A[i] * r;
23dfcca431 2011-02-23        kinaba: 		}
23dfcca431 2011-02-23        kinaba: 	}
23dfcca431 2011-02-23        kinaba: 	return A;
23dfcca431 2011-02-23        kinaba: }