Go: package main

import (
	"fmt"
	"math"
)

func main() {
	var num int
	fmt.Scan(&num)
	for i := 0; i < num; i++ {
		var x1, y1, x2, y2, x3, y3 float64
		fmt.Scanf("%f%f%f%f%f%f", &x1, &y1, &x2, &y2, &x3, &y3)

		a := 2 * (x2 - x1)
		b := 2 * (y2 - y1)
		c := 2 * (x3 - x1)
		d := 2 * (y3 - y1)
		e := x2*x2 - x1*x1 + y2*y2 - y1*y1
		f := x3*x3 - x1*x1 + y3*y3 - y1*y1

		// 直線1: ax + by - e = 0
		// 直線2: cx + dy - f = 0
		px := (d*e - b*f) / (a*d - b*c)
		py := (c*e - a*f) / (b*c - a*d)
		r := math.Sqrt((px-x1)*(px-x1) + (py-y1)*(py-y1))
		fmt.Printf("%.3f %.3f %.3f\n", px, py, r)
	}
}

C++: #include <vector>
#include <iostream>
#include <vector>
#include <queue>
#include <algorithm>
#include <map>
#include <climits>
#include <cmath>
#include <cfloat>
#include <cstdio>
#include <ctime>
#include <functional>

using namespace std;

#define EPS 1e-6

class xy{
	
public:
	double x;
	double y;

	xy(){

	}
	
	xy(double xx, double yy){
		x = xx;
		y = yy;
	}

	xy(const xy &v){
		x = v.x;
		y = v.y;
	}

	xy& operator=(const xy &v){
		x = v.x;
		y = v.y;
		return *this;
	}

	xy operator+(const xy &v) const{
		return xy(this->x+v.x, this->y+v.y);
		
	}
	
	xy operator-(const xy &v) const{
		return xy(this->x-v.x, this->y-v.y);
	}
	
	void operator+=(const xy &v){
		x+=v.x;
		y+=v.y;
	}
	void operator-=(const xy &v){
		x-=v.x;
		y-=v.y;
	}
/* something wrong this operator for sorting
	bool operator<(const xy &v){
		if(x!=v.x) return x < v.x;
		return y < v.y;
	}
	bool operator>(const xy &v){
		if(x!=v.x) return x > v.x;
		return y > v.y;
	}
*/

};
//for sorting
bool comp_xy(const xy &a, const xy &b){
	if(a.x != b.x) return a.x < b.x;
	return a.y < b.y;
}

//u,v : vector O = (0,0)
double cross(const xy &u, const xy &v){
	return u.x*v.y - u.y*v.x;
}

//u,v : vector O = (0,0)
double dot(const xy &u, const xy &v){
	return u.x*v.x + u.y*v.y;
}

//distance between two points
double dist_p_p(const xy &a, const xy &b){
	return sqrt( fabs(dot(a-b, a-b)) );
}

//distance between a point and a line segment
double dist_p_ls(const xy &p, const xy &s1, const xy &s2){
	xy vl = s2 - s1;
	xy vp = p - s1;

	return fabs( cross(vl, vp) / sqrt( dot(vl, vl) ) );
}


class Smallest_Enclosing_Disk{
	const static int INF = (1<<30);
	priority_queue< pair<int,int> > pq;
	vector<xy> v;

	typedef struct{
		double r;
		xy O;
	}disk;
	
	pair<vector<int>, disk> make_disk_from_2_points(vector<int> R){
		disk ret;
		ret.O.x = (v[R[0]].x + v[R[1]].x) * 0.5;
		ret.O.y = (v[R[0]].y + v[R[1]].y) * 0.5;
		ret.r = dist_p_p(v[R[0]], v[R[1]]) * 0.5;
		if(ret.r == 0) R.erase(R.begin()+1);
		return make_pair(R, ret);
	}

	pair<vector<int>, disk> make_disk_from_3_points(vector<int> R){
		disk ret;
		
		if(dot(v[R[1]] - v[R[0]], v[R[2]] - v[R[0]]) <= 0){
			pair<vector<int>, disk> d = make_disk_from_2_points( vector<int>{R[1], R[2]});
			R = d.first;
			ret = d.second;
			
		}else if(dot(v[R[0]] - v[R[1]], v[R[2]] - v[R[1]]) <= 0){
			pair<vector<int>, disk> d = make_disk_from_2_points( vector<int>{R[0], R[2]});
			R = d.first;
			ret = d.second;
			
		}else if(dot(v[R[0]] - v[R[2]], v[R[1]] - v[R[2]]) <= 0){
			pair<vector<int>, disk> d = make_disk_from_2_points( vector<int>{R[0], R[1]});
			R = d.first;
			ret = d.second;
			
		}else{
			double A = dot(v[R[0]], v[R[0]]);
			double B = dot(v[R[1]], v[R[1]]);
			double C = dot(v[R[2]], v[R[2]]);

			//r(x,y) = { A(B+C-A)r_a + B(C+A-B)r_b + C(A+B-C)r_c } / 16S^2
			//16S^2 = A(B+C-A) + B(C+A-B) + C(A+B-C)

			double S_ = A*(B+C-A) + B*(C+A-B) + C*(A+B-C);
			ret.O.x = (A*(B+C-A) * v[R[0]].x + B*(C+A-B) * v[R[1]].x + C*(A+B-C) * v[R[2]].x) / S_;
			ret.O.y = (A*(B+C-A) * v[R[0]].y + B*(C+A-B) * v[R[1]].y + C*(A+B-C) * v[R[2]].y) / S_;
			ret.r = dist_p_p( v[R[0]], ret.O ) * 0.5;
		}

		return make_pair( R, ret);
	}
	
	//R := index of points which are on disk border
	//S := smallest enclosing disk
	//P := new added point
	pair<vector<int>, disk> D(vector<int> R, disk S, int P){
		//P is inside of S
		if( (dist_p_p(v[P], S.O) <= S.r) && (R.size() > 0) ){
			return make_pair(R, S);
		}
		disk ret;
		if(R.size() == 0){
			R.push_back(P);
			ret.O = v[P];
			ret.r = 0.0;
		}else if(R.size() == 1){
			R.push_back(P);
			
			pair<vector<int>, disk> d = make_disk_from_2_points(R);
			R = d.first;
			ret = d.second;
			
		}else if(R.size() == 2){
			R.push_back(P);
			pair<vector<int>, disk> d = make_disk_from_3_points(R);
			R = d.first;
			ret = d.second;
		}else if(R.size() == 3){
			vector<int> tmp_R;
			disk tmp_S;
			for(int i=0; i<3; i++){
				pair<vector<int>, disk> d = make_disk_from_3_points(vector<int>{R[i%3],R[(i+1)%3], P});

				//disk can't enclose every point
				if(dist_p_p(v[R[(i+2)%3]], d.second.O) > d.second.r){
					d.second.r = INF;
				}
				if(i==0 || d.second.r < tmp_S.r){
					tmp_R = d.first;
					tmp_S = d.second;
				}
			}
			R = tmp_R;
			ret = tmp_S;
		}
		return make_pair(R,ret);
	}

public:
	//r, xy
	Smallest_Enclosing_Disk(vector<xy> &v){
		this->v = v;
		//random key, index
		for(int i=0; i<v.size(); i++){
			pq.push( make_pair(rand(), i) );
		}
	}

	double calc(){
		vector<int> R;
		disk ret;
		while(!pq.empty()){
			pair<vector<int>, disk> d = D(R, ret, pq.top().second);
			pq.pop();
			R = d.first;
			ret = d.second;
		}
		return ret.r;
	}
};



//test
int main(){
	srand(time(NULL));
	int n;
	cin >> n;
	for(int i=0; i<n; i++){
		vector<xy> v(3);
		for(int j=0; j<3; j++){
			double x,y;
			cin >> x >> y;
			v[j] = xy(x,y);
		}
		xy O;
		double r;
		// rA = v[0]
		// rB = v[1]
		// rC = v[2]
		xy a = v[2] - v[1];
		xy b = v[2] - v[0];
		xy c = v[1] - v[0];
		double A = dot(a, a);
		double B = dot(b, b);
		double C = dot(c, c);

		//r(x,y) = { A(B+C-A)r_a + B(C+A-B)r_b + C(A+B-C)r_c } / 16S^2
		//16S^2 = A(B+C-A) + B(C+A-B) + C(A+B-C)

		double S_ = A*(B+C-A) + B*(C+A-B) + C*(A+B-C);
		O.x = (A*(B+C-A) * v[0].x + B*(C+A-B) * v[1].x + C*(A+B-C) * v[2].x) / S_;
		O.y = (A*(B+C-A) * v[0].y + B*(C+A-B) * v[1].y + C*(A+B-C) * v[2].y) / S_;
		r = dist_p_p( v[0], O );

		printf("%.3f %.3f %.3f\n", O.x, O.y, r);
	}
}

C++: #include <bits/stdc++.h>

using namespace std;

#define EPS (1e-5)
#define equals(a, b) (fabs(a - b) < EPS)

struct Point {
    double x, y;

    Point(){}
    Point(double x, double y) : x(x), y(y) {}

    Point operator + (const Point &p) const {
        return Point(x + p.x, y + p.y);
    }

    Point operator - (const Point &p) const {
        return Point(x - p.x, y - p.y);
    }
    
    Point operator * (const double &k) const {
        return Point(x * k, y * k);
    }
};

typedef Point Vector;

double dot(const Point &a, const Point &b)
{
    return a.x * b.x + a.y * b.y;
}

double cross(const Point &a, const Point &b)
{
    return a.x * b.y - b.x * a.y;
}

double norm(const Point &p)
{
    return dot(p, p);
} 

double abs(const Point &p)
{
    return sqrt(norm(p));
}

struct Line {
    Point s, t;
    Line(){}
    Line(Point s, Point t) : s(s), t(t) {}
};

Line getPerpendicularBisector(const Point &a, const Point &b)
{
    double cx = (a.x + b.x) / 2.0;
    double cy = (a.y + b.y) / 2.0;
    Point p = Point(cx + (a.y - b.y), cy + (b.x - a.x));
    return Line(Point(cx, cy),p);
}

Point crosspointLL(const Line &a, const Line &b)
{
    Vector va = a.t - a.s, vb = b.t - b.s;
    double d = cross(vb, va);
    if (abs(d) < EPS) return b.s;
    return a.s+va*cross(vb, b.t - a.s) * (1.0 / d);
}

int main()
{
    Point p[3];
    int N;
    cin >> N;
    while (N--) {
        for (int i = 0; i < 3; i++) {
            cin >> p[i].x >> p[i].y;
        }
        Line l[2];
        for (int i = 0; i < 2; i++) {
            l[i] = getPerpendicularBisector(p[i], p[i+1]);
        }
        Point c = crosspointLL(l[0], l[1]);
        double r = abs(c - p[0]);
        printf("%.3f %.3f %.3f\n", c.x, c.y, r);
    }
    return 0;
}