4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: //-------------------------------------------------------------
4fd800b3a8 2011-02-23        kinaba: // The circle passing three points
4fd800b3a8 2011-02-23        kinaba: //
4fd800b3a8 2011-02-23        kinaba: // Verified by
4fd800b3a8 2011-02-23        kinaba: //   - AOJ 0010
4fd800b3a8 2011-02-23        kinaba: //-------------------------------------------------------------
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: vector<double> solve_linear_eq( int n, vector< vector<double> > M, const vector<double>& V )
4fd800b3a8 2011-02-23        kinaba: {
4fd800b3a8 2011-02-23        kinaba: 	vector<double> A(V);
4fd800b3a8 2011-02-23        kinaba: 	for(int i=0; i<n; ++i)
4fd800b3a8 2011-02-23        kinaba: 	{
4fd800b3a8 2011-02-23        kinaba: 		// pivot
4fd800b3a8 2011-02-23        kinaba: 		if( M[i][i] == 0 )
4fd800b3a8 2011-02-23        kinaba: 			for(int j=i+1; j<n; ++j)
4fd800b3a8 2011-02-23        kinaba: 				if( M[j][i] != 0 )
4fd800b3a8 2011-02-23        kinaba: 					{swap(M[i], M[j]); swap(A[i], A[j]); break;}
4fd800b3a8 2011-02-23        kinaba: 		if( M[i][i] == 0 )
4fd800b3a8 2011-02-23        kinaba: 			throw "no anser";
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 		// M[i][i] <-- 1
4fd800b3a8 2011-02-23        kinaba: 		double p = M[i][i];
4fd800b3a8 2011-02-23        kinaba: 		for(int j=i; j<n; ++j)
4fd800b3a8 2011-02-23        kinaba: 			M[i][j] /= p;
4fd800b3a8 2011-02-23        kinaba: 		A[i] /= p;
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 		// M[*][i] <-- 0
4fd800b3a8 2011-02-23        kinaba: 		for(int j=0; j<n; ++j) if(j!=i)
4fd800b3a8 2011-02-23        kinaba: 		{
4fd800b3a8 2011-02-23        kinaba: 			double r = M[j][i];
4fd800b3a8 2011-02-23        kinaba: 			for(int k=i; k<n; ++k)
4fd800b3a8 2011-02-23        kinaba: 				M[j][k] -= M[i][k] * r;
4fd800b3a8 2011-02-23        kinaba: 			A[j] -= A[i] * r;
4fd800b3a8 2011-02-23        kinaba: 		}
4fd800b3a8 2011-02-23        kinaba: 	}
4fd800b3a8 2011-02-23        kinaba: 	return A;
4fd800b3a8 2011-02-23        kinaba: }
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: void circle_3pt( CMP p1, CMP p2, CMP p3, CMP* c, double* r )
4fd800b3a8 2011-02-23        kinaba: {
4fd800b3a8 2011-02-23        kinaba: 	p2 -= p1;
4fd800b3a8 2011-02-23        kinaba: 	p3 -= p1;
4fd800b3a8 2011-02-23        kinaba: 	// c == p2/2 + A p2 i == p3/2 + B p3 i
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 	vector< vector<double> > M(2, vector<double>(2));
4fd800b3a8 2011-02-23        kinaba: 	vector<double> V(2);
4fd800b3a8 2011-02-23        kinaba: 	M[0][0] = -p2.imag(); M[0][1] = +p3.imag(); V[0] = (p3.real()-p2.real())/2.0;
4fd800b3a8 2011-02-23        kinaba: 	M[1][0] = +p2.real(); M[1][1] = -p3.real(); V[1] = (p3.imag()-p2.imag())/2.0;
4fd800b3a8 2011-02-23        kinaba: 	V = solve_linear_eq(2, M, V);
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 	double A=V[0], B=V[1];
4fd800b3a8 2011-02-23        kinaba: 	*c = p1 + p2/2.0 + A * p2 * CMP(0,1);
4fd800b3a8 2011-02-23        kinaba: 	*r = abs(*c - p1);
4fd800b3a8 2011-02-23        kinaba: }