Go: package main

import (
	"bufio"
	"fmt"
	"io"
	"math"
	"os"
	"strconv"
	"strings"

)


func main() {
	solve(os.Stdin, os.Stdout)
}

type ReaderEx struct {
	baseScanner *bufio.Scanner
	words       []string
	nextWordIdx func() int
	initWordIdx func(int)
}

func NewScannerEx(stdin io.Reader) *ReaderEx {
	var this ReaderEx
	var idx int
	var wordsCnt int
	this.baseScanner = bufio.NewScanner(stdin)
	this.baseScanner.Buffer(make([]byte, 0), 1E+11)
	this.nextWordIdx = func() int {
		cur := idx
		idx = (idx + 1) % wordsCnt
		return cur
	}
	this.initWordIdx = func(wc int) {
		idx = 0
		wordsCnt = wc
	}
	return &this
}

func (me *ReaderEx) Scan(doSplit bool) (isScanned bool) {
	isScanned = me.baseScanner.Scan()
	if doSplit {
		me.words = strings.Fields(me.baseScanner.Text())
	} else {
		me.words = []string{me.baseScanner.Text()}
	}
	me.initWordIdx(len(me.words))
	return
}

func (me *ReaderEx) nextBytes() []byte {
	return []byte(me.words[me.nextWordIdx()])
}
func (me *ReaderEx) nextString() string {
	return me.words[me.nextWordIdx()]
}

type StringArray []string

func (me *ReaderEx) nextStringArray() StringArray {
	return me.words[me.nextWordIdx():]
}

func (me StringArray) toStrings() []string {
	return me
}

func (me StringArray) toInts() []int {
	strValues := me.toStrings()
	values := make([]int, len(strValues))
	for i, v := range strValues {
		n, err := strconv.Atoi(v)
		if err != nil {
			panic(err)
		}
		values[i] = n
	}
	return values
}
func (me *ReaderEx) nextInt() int {
	n, err := strconv.Atoi(me.words[me.nextWordIdx()])

	defer func() {
		if err != nil {
			panic(fmt.Sprintf("getInt(%d)", n))
		}
	}()
	return n
}

func pow(x, y, mod int) (ans int) {

	ans = 1
	for bit := 1 << uint(math.Floor(math.Log2(float64(y)))); ; bit >>= 1 {
		if y&bit != 0 {
			ans *= x
			if mod > 0 {
				ans %= mod
			}
		}
		if bit <= 1 {
			return
		}
		ans *= ans
		if mod > 0 {
			ans %= mod
		}
	}
}

type factType struct {
	prime int
	count int
}

func fact(n int) (facts []*factType) {

	for i := 2; i*i <= n; i++ {
		if n%i == 0 {
			fa := factType{prime: i}
			for n%fa.prime == 0 {
				n /= fa.prime
				fa.count++
			}
			facts = append(facts, &fa)
		}
	}
	if n > 1 {
		facts = append(facts, &factType{prime: n, count: 1})
	}
	return facts
}
func dec2bin(n int) (ret []bool) {
	for bit := 1 << uint(math.Floor(math.Log2(float64(n)))); ; bit >>= 1 {
		ret = append(ret, n&bit != 0)

		if bit <= 1 {
			return
		}
	}

}

func abs(v int) int {
	if v < 0 {
		v = -v
	}
	return v
}

func b2i(b bool) (i int) {
	if b {
		i = 1
	} else {
		i = 0
	}
	return
}

func solve(stdin io.Reader, stdout io.Writer) {

	r := NewScannerEx(stdin)
	r.Scan(true)
	n := r.nextInt()

	type Vec struct {
		x, y float64
	}
	vec := func(a, b Vec) Vec { return Vec{b.x - a.x, b.y - a.y} }
	unv := func(a Vec) Vec { return Vec{-a.y, a.x} }
	cdot := func(a, b Vec) float64 { return a.x*b.x + a.y*b.y }
	times := func(t float64, a Vec) Vec { return Vec{t * a.x, t * a.y} }
	add := func(a, b Vec) Vec { return Vec{a.x + b.x, a.y + b.y} }
	for i := 0; i < n; i++ {

		r.Scan(true)
		nextFloat64 := func() (ret float64) { fmt.Sscanf(r.nextString(), "%f", &ret); return }
		a1, a2, a3 := Vec{nextFloat64(), nextFloat64()}, Vec{nextFloat64(), nextFloat64()}, Vec{nextFloat64(), nextFloat64()}
		unv23 := unv(vec(a2, a3))
		v12 := vec(a1, a2)
		v31 := vec(a3, a1)
		s := cdot(v12, v31) / 2 / cdot(v12, unv23)
		c := add(times(.5, add(a2, a3)), times(s, unv23))
		v1c := vec(a1, c)
		fmt.Fprintf(stdout, "%.3f %.3f %.3f\n", c.x, c.y, math.Sqrt(v1c.x*v1c.x+v1c.y*v1c.y))
	}
}

C++: #include <iostream>
#include <complex>
#include <cstdio>
#include <cmath>
using namespace std;

double area(const complex<double> &p0, const complex<double> &p1, const complex<double> &p2) {
    complex<double> _p0 (p0.real() - p1.real(), p0.imag() - p1.imag());
    complex<double> _p1 (p1.real() - p2.real(), p1.imag() - p2.imag());

    double a = (_p0.real() * _p1.imag() - _p0.imag() * _p1.real()) * 0.5;
    return (a > 0) ? a : -a;
}

double dist2(const complex<double> &p0, const complex<double> &p1) {
    return (p0.real() - p1.real()) * (p0.real() - p1.real())
        + (p0.imag() - p1.imag()) * (p0.imag() - p1.imag());
}

int main () {
    int n;
    cin >> n;
    while (n--) {
        complex<double> p[3];
        for (int i = 0; i < 3; i++) {
            double x, y;
            cin >> x >> y;
            p[i] = complex<double> (x, y);
        }

        double d2[3], d[3];
        for (int i = 0; i < 3; i++) {
            d2[i] = dist2(p[(i+1)%3], p[(i+2)%3]);
            d[i] = sqrt(d2[i]);
        }

        double s = area(p[0], p[1], p[2]);
        double ss = s * s * 16;

        double e[3];
        e[0] = d2[0] * (-d2[0] + d2[1] + d2[2]) / ss;
        e[1] = d2[1] * (d2[0] - d2[1] + d2[2]) / ss;
        e[2] = d2[2] * (d2[0] + d2[1] - d2[2]) / ss;

        double x = 0, y = 0;
        for (int i = 0; i < 3; i++) {
            x += e[i] * p[i].real();
            y += e[i] * p[i].imag();
        }
        double r = d[0] * d[1] * d[2] /
            sqrt((d[0]+d[1]+d[2]) * (-d[0]+d[1]+d[2]) * (d[0]-d[1]+d[2]) * (d[0]+d[1]-d[2]));

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

C++: #include<iostream>
#include<cmath>
using namespace std;

int main()
{
	int n;
	double x1,y1,x2,y2,x3,y3;
	cin >> n;

	for(int i=0;i<n;i++)
	{
		cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3;

		double a1 = 2*(x2 - x1);
		double b1 = 2*(y2 - y1);
		double c1 = x1*x1 - x2*x2 + y1*y1 - y2*y2;

		double a2 = 2*(x3 -x1);
		double b2 = 2*(y3 - y1);
		double c2 = x1*x1 - x3*x3 + y1*y1 - y3*y3;

		double x = (b1*c2 - b2*c1) / (a1*b2 - a2*b1);
		double y = (c1*a2 - c2*a1) / (a1*b2 - a2*b1);

		double r_2 = (x1 - x)*(x1 - x) + (y1 - y)*(y1 - y);
		double r = sqrt(r_2);

		cout << x << " " << y << " " << r << endl;
	}

	return 0;
}

C++: #include<iostream>
#include<vector>
#include<algorithm>
#include<iomanip>
using namespace std;

const double EPS = 1e-8;

vector<double> gauss_jordan(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]);

        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> x(n);
    for(int i = 0 ; i < n ; ++i){
        x[i] = B[i][n];
    }

    return x;
}

int main(){
    int n;
    cin >> n;

    while(n--){
        vector< vector<double> > A(2, vector<double>(2));
        vector<double> b(2);
        vector<double> x(3), y(3);
        for(int i = 0 ; i < 3 ; ++i){
            cin >> x[i] >> y[i];
        }
        A[0][0] = 2.0 * (x[0] - x[1]);
        A[0][1] = 2.0 * (y[0] - y[1]);
        b[0] = x[0] * x[0] + y[0] * y[0] - x[1] * x[1] - y[1] * y[1];
        A[1][0] = 2.0 * (x[1] - x[2]);
        A[1][1] = 2.0 * (y[1] - y[2]);
        b[1] = x[1] * x[1] + y[1] * y[1] - x[2] * x[2] - y[2] * y[2];
        vector<double> p = gauss_jordan(A, b);
        double r = sqrt((x[0] - p[0]) * (x[0] - p[0]) + (y[0] - p[1]) * (y[0] - p[1]));
        cout << fixed << setprecision(3) << p[0] << ' ' << p[1] << ' ' << r << endl;
    }

    return 0;
}