Paste code

#include <iostream>
#include <cmath>
#include <complex>
#include <limits>
#include <cstdio>
#include <iomanip>
#include <cstdlib>

using namespace std;
typedef long double ld;
typedef std::complex <ld> r2vec;

ld R(r2vec A, r2vec B, r2vec C) {
    ld a(abs(B-C)), b(abs(A-C)), c(abs(A-B));
    ld S = sqrt((a+b+c)*(-a+b+c)*(a-b+c)*(a+b-c))/4;
    return a*b*c/(4*S);
}

ld scalar(r2vec a, r2vec b) {
    return real(a)*real(b) + imag(a)*imag(b);
}

r2vec U(r2vec A, r2vec B, r2vec C) {
    ld a(abs(B-C)), b(abs(A-C)), c(abs(A-B));
    ld S = sqrt((a+b+c)*(-a+b+c)*(a-b+c)*(a+b-c))/4;
    ld aa = a*a*scalar(A-B, A-C);
    ld ab = b*b*scalar(B-A, B-C);
    ld ac = c*c*scalar(C-A, C-B);
    return (aa*A+ab*B+ac*C)/(8*S*S);
}

ld S(ld R, int n) {
    return (n*R*R*sin(2*M_PI/n))/2;
}

bool eps_eq(r2vec a, r2vec b) {
    return abs(a-b) < 1e-4;
}

ld find_sq(r2vec A, r2vec B, r2vec C, int n) {
    ld r = R(A, B, C);
    r2vec u = U(A, B, C);
    ld p = acos(real((A-u)/r));
    int count = 0;
    A -= u; B -= u; C -= u;
    for (int i = 0; i < n; ++i) {
        r2vec v = r*r2vec(cos(p+2*M_PI*i/n), sin(p+2*M_PI*i/n));
        count += eps_eq(v, A) + eps_eq(v, B) + eps_eq(v, C);
    }
    if (count > 2) return S(r, n);
    else return numeric_limits <ld>::max();
}

double find_answer(r2vec A, r2vec B, r2vec C) {
    ld answer = numeric_limits <ld>::max();
    for (int i = 3; i < 100; ++i) {
        answer = min(find_sq(A, B, C, i), answer);
    }
    return answer;
}

ld gen_rand(ld low, ld high) {
    ld range = high-low;
    return  rand() * range / RAND_MAX + low;
}

void test_it() {
    while (true) {
        ld low = -1e3;
        ld high = 1e3;
        r2vec u = r2vec(gen_rand(low, high),gen_rand(low, high));
        ld r = gen_rand(0, high);
        int n = rand() % 98 + 3;
        ld a(rand()%n), b(rand()%n), c(rand()%n);
        if (a == b || a == c || b == c) {
            continue;
        }
        r2vec im_base(0, 2*M_PI/n);
        r2vec A(u+r*exp(a*im_base)),
              B(u+r*exp(b*im_base)),
              C(u+r*exp(c*im_base));
        ld s = S(r, n);
        ld sapp = find_answer(A, B, C);
        if (abs(s - sapp) > 1e-6) {
            cout << s << ' ' << sapp << '\n';
            cout << fixed << setprecision(6) << u << ' ' << r << ' ' << n << ' ' << A << ' ' << B << ' ' << C << '\n';
        }
    }
}

int main() {
    ld x, y;
    cin >> x >> y;
    r2vec A(x, y);
    cin >> x >> y;
    r2vec B(x, y);
    cin >> x >> y;
    r2vec C(x, y);
    cout << fixed << setprecision(8) << find_answer(A, B, C);
//    test_it();
}