Problem 966: Triangle Circle Intersection

Problem 966 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>

namespace {

constexpr long double kPi = 3.141592653589793238462643383279502884L;
constexpr long double kEps = 1e-12L;

struct Vec2 {
    long double x;
    long double y;

    Vec2 operator+(const Vec2& other) const {
        return {x + other.x, y + other.y};
    }

    Vec2 operator-(const Vec2& other) const {
        return {x - other.x, y - other.y};
    }

    Vec2 operator*(long double s) const {
        return {x * s, y * s};
    }
};

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

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

long double norm2(const Vec2& v) {
    return dot(v, v);
}

long double dist(const Vec2& a, const Vec2& b) {
    return std::sqrt(norm2(a - b));
}

struct Triangle {
    Vec2 a;
    Vec2 b;
    Vec2 c;
    long double area;
    long double radius;
    long double max_side;
};

Triangle make_triangle(int s1, int s2, int s3) {
    std::array<int, 3> sides{s1, s2, s3};
    std::sort(sides.begin(), sides.end());

    const long double a = static_cast<long double>(sides[0]);
    const long double b = static_cast<long double>(sides[1]);
    const long double c = static_cast<long double>(sides[2]);

    const Vec2 A{0.0L, 0.0L};
    const Vec2 B{c, 0.0L};
    const long double x = (b * b + c * c - a * a) / (2.0L * c);
    long double y2 = b * b - x * x;
    if (y2 < 0.0L && y2 > -1e-18L) {
        y2 = 0.0L;
    }
    const Vec2 C{x, std::sqrt(std::max(0.0L, y2))};

    const long double area = std::abs(cross(B - A, C - A)) * 0.5L;
    const long double radius = std::sqrt(area / kPi);
    const long double max_side = std::max({dist(A, B), dist(B, C), dist(C, A)});
    return {A, B, C, area, radius, max_side};
}

bool point_in_triangle(const Triangle& tri, const Vec2& p) {
    const long double c1 = cross(tri.b - tri.a, p - tri.a);
    const long double c2 = cross(tri.c - tri.b, p - tri.b);
    const long double c3 = cross(tri.a - tri.c, p - tri.c);
    return c1 >= -1e-14L && c2 >= -1e-14L && c3 >= -1e-14L;
}

Vec2 nearest_on_segment(const Vec2& p, const Vec2& u, const Vec2& v) {
    const Vec2 uv = v - u;
    const long double d = norm2(uv);
    if (d <= 0.0L) {
        return u;
    }
    long double t = dot(p - u, uv) / d;
    if (t < 0.0L) {
        t = 0.0L;
    } else if (t > 1.0L) {
        t = 1.0L;
    }
    return u + uv * t;
}

Vec2 project_to_triangle(const Triangle& tri, const Vec2& p) {
    if (point_in_triangle(tri, p)) {
        return p;
    }

    const Vec2 candidates[6] = {
        tri.a,
        tri.b,
        tri.c,
        nearest_on_segment(p, tri.a, tri.b),
        nearest_on_segment(p, tri.b, tri.c),
        nearest_on_segment(p, tri.c, tri.a),
    };

    Vec2 best = candidates[0];
    long double best_d = norm2(p - best);
    for (int i = 1; i < 6; ++i) {
        const long double cur = norm2(p - candidates[i]);
        if (cur < best_d) {
            best_d = cur;
            best = candidates[i];
        }
    }
    return best;
}

long double segment_contribution(const Vec2& a, const Vec2& b, long double r2) {
    Vec2 pts[4];
    int n = 0;
    pts[n++] = a;

    const Vec2 d = b - a;
    const long double qa = dot(d, d);
    const long double qb = 2.0L * dot(a, d);
    const long double qc = dot(a, a) - r2;
    const long double disc = qb * qb - 4.0L * qa * qc;

    if (disc > 1e-18L) {
        const long double sdisc = std::sqrt(disc);
        long double t1 = (-qb - sdisc) / (2.0L * qa);
        long double t2 = (-qb + sdisc) / (2.0L * qa);
        if (t1 > t2) {
            std::swap(t1, t2);
        }
        if (t1 > kEps && t1 < 1.0L - kEps) {
            pts[n++] = a + d * t1;
        }
        if (t2 > kEps && t2 < 1.0L - kEps && std::abs(t2 - t1) > 1e-12L) {
            pts[n++] = a + d * t2;
        }
    }

    pts[n++] = b;

    long double acc = 0.0L;
    for (int i = 0; i + 1 < n; ++i) {
        const Vec2 p = pts[i];
        const Vec2 q = pts[i + 1];
        const Vec2 mid{(p.x + q.x) * 0.5L, (p.y + q.y) * 0.5L};
        if (norm2(mid) <= r2 + 1e-13L) {
            acc += 0.5L * cross(p, q);
        } else {
            acc += 0.5L * r2 * std::atan2(cross(p, q), dot(p, q));
        }
    }
    return acc;
}

long double overlap_area(const Triangle& tri, const Vec2& center) {
    const long double r = tri.radius;
    const long double r2 = r * r;

    const Vec2 p0 = tri.a - center;
    const Vec2 p1 = tri.b - center;
    const Vec2 p2 = tri.c - center;

    long double area = 0.0L;
    area += segment_contribution(p0, p1, r2);
    area += segment_contribution(p1, p2, r2);
    area += segment_contribution(p2, p0, r2);
    return std::abs(area);
}

struct EvalPoint {
    Vec2 p;
    long double v;
};

EvalPoint nelder_mead(const Triangle& tri, const Vec2& start, long double step, int iterations) {
    auto eval = [&](const Vec2& raw) -> EvalPoint {
        const Vec2 projected = project_to_triangle(tri, raw);
        return {projected, overlap_area(tri, projected)};
    };

    std::array<EvalPoint, 3> simplex = {
        eval(start),
        eval({start.x + step, start.y}),
        eval({start.x, start.y + step}),
    };

    for (int it = 0; it < iterations; ++it) {
        std::sort(simplex.begin(), simplex.end(),
                  [](const EvalPoint& lhs, const EvalPoint& rhs) { return lhs.v > rhs.v; });

        const Vec2 centroid{
            (simplex[0].p.x + simplex[1].p.x) * 0.5L,
            (simplex[0].p.y + simplex[1].p.y) * 0.5L,
        };

        const Vec2 reflected{
            centroid.x + (centroid.x - simplex[2].p.x),
            centroid.y + (centroid.y - simplex[2].p.y),
        };
        const EvalPoint er = eval(reflected);

        if (er.v > simplex[0].v) {
            const Vec2 expanded{
                centroid.x + 2.0L * (er.p.x - centroid.x),
                centroid.y + 2.0L * (er.p.y - centroid.y),
            };
            const EvalPoint ee = eval(expanded);
            simplex[2] = (ee.v > er.v ? ee : er);
            continue;
        }

        if (er.v > simplex[1].v) {
            simplex[2] = er;
            continue;
        }

        const Vec2 contracted{
            centroid.x + 0.5L * (simplex[2].p.x - centroid.x),
            centroid.y + 0.5L * (simplex[2].p.y - centroid.y),
        };
        const EvalPoint ec = eval(contracted);
        if (ec.v > simplex[2].v) {
            simplex[2] = ec;
            continue;
        }

        simplex[1] = eval({
            simplex[0].p.x + 0.5L * (simplex[1].p.x - simplex[0].p.x),
            simplex[0].p.y + 0.5L * (simplex[1].p.y - simplex[0].p.y),
        });
        simplex[2] = eval({
            simplex[0].p.x + 0.5L * (simplex[2].p.x - simplex[0].p.x),
            simplex[0].p.y + 0.5L * (simplex[2].p.y - simplex[0].p.y),
        });
    }

    std::sort(simplex.begin(), simplex.end(),
              [](const EvalPoint& lhs, const EvalPoint& rhs) { return lhs.v > rhs.v; });
    return simplex[0];
}

long double maximize_overlap(int s1, int s2, int s3) {
    const Triangle tri = make_triangle(s1, s2, s3);

    std::vector<Vec2> seeds;
    seeds.reserve(24);

    const int grid = 5;
    for (int i = 0; i <= grid; ++i) {
        for (int j = 0; j + i <= grid; ++j) {
            const long double u = static_cast<long double>(i) / static_cast<long double>(grid);
            const long double v = static_cast<long double>(j) / static_cast<long double>(grid);
            const long double w = 1.0L - u - v;
            seeds.push_back({
                u * tri.a.x + v * tri.b.x + w * tri.c.x,
                u * tri.a.y + v * tri.b.y + w * tri.c.y,
            });
        }
    }

    const Vec2 centroid{
        (tri.a.x + tri.b.x + tri.c.x) / 3.0L,
        (tri.a.y + tri.b.y + tri.c.y) / 3.0L,
    };
    seeds.push_back(centroid);

    const long double side_a = dist(tri.b, tri.c);
    const long double side_b = dist(tri.c, tri.a);
    const long double side_c = dist(tri.a, tri.b);
    const long double ws = side_a + side_b + side_c;
    seeds.push_back({
        (side_a * tri.a.x + side_b * tri.b.x + side_c * tri.c.x) / ws,
        (side_a * tri.a.y + side_b * tri.b.y + side_c * tri.c.y) / ws,
    });

    std::vector<EvalPoint> evaluated;
    evaluated.reserve(seeds.size());
    for (const Vec2& p : seeds) {
        const Vec2 projected = project_to_triangle(tri, p);
        evaluated.push_back({projected, overlap_area(tri, projected)});
    }

    std::sort(evaluated.begin(), evaluated.end(),
              [](const EvalPoint& lhs, const EvalPoint& rhs) { return lhs.v > rhs.v; });

    EvalPoint best = evaluated[0];
    const int tries = std::min(4, static_cast<int>(evaluated.size()));
    const long double step = tri.max_side * 0.18L;
    for (int i = 0; i < tries; ++i) {
        const EvalPoint cur = nelder_mead(tri, evaluated[static_cast<std::size_t>(i)].p, step, 90);
        if (cur.v > best.v) {
            best = cur;
        }
    }

    const EvalPoint refined = nelder_mead(tri, best.p, tri.max_side * 0.05L, 50);
    if (refined.v > best.v) {
        best = refined;
    }

    if (best.v < 0.0L) {
        return 0.0L;
    }
    if (best.v > tri.area) {
        return tri.area;
    }
    return best.v;
}

bool nearly_equal(long double x, long double y, long double eps) {
    return std::abs(x - y) <= eps;
}

bool run_checkpoints() {
    const long double i345 = maximize_overlap(3, 4, 5);
    if (!nearly_equal(i345, 4.593049L, 2e-6L)) {
        std::cerr << "Checkpoint failed for I(3,4,5): got " << std::setprecision(12) << i345 << '\n';
        return false;
    }

    const long double i346 = maximize_overlap(3, 4, 6);
    if (!nearly_equal(i346, 3.552564L, 2e-6L)) {
        std::cerr << "Checkpoint failed for I(3,4,6): got " << std::setprecision(12) << i346 << '\n';
        return false;
    }

    const long double i_perm = maximize_overlap(5, 3, 4);
    if (!nearly_equal(i_perm, i345, 1e-10L)) {
        std::cerr << "Permutation checkpoint failed for I(5,3,4)" << '\n';
        return false;
    }

    return true;
}

std::vector<std::array<int, 3>> enumerate_triangles() {
    std::vector<std::array<int, 3>> all;
    for (int a = 1; a <= 200; ++a) {
        for (int b = a; b <= 200 - a; ++b) {
            const int max_c = 200 - a - b;
            for (int c = b; c <= max_c; ++c) {
                if (a + b <= c) {
                    continue;
                }
                all.push_back({a, b, c});
            }
        }
    }
    return all;
}

long double solve() {
    const std::vector<std::array<int, 3>> triangles = enumerate_triangles();

    unsigned int thread_count = std::thread::hardware_concurrency();
    if (thread_count == 0) {
        thread_count = 4;
    }
    thread_count = std::min<unsigned int>(thread_count, static_cast<unsigned int>(triangles.size()));
    if (thread_count == 0) {
        return 0.0L;
    }

    std::atomic<std::size_t> next_idx{0};
    std::vector<long double> partial(thread_count, 0.0L);

    auto worker = [&](unsigned int tid) {
        long double sum = 0.0L;
        long double comp = 0.0L;

        while (true) {
            const std::size_t idx = next_idx.fetch_add(1, std::memory_order_relaxed);
            if (idx >= triangles.size()) {
                break;
            }
            const auto& t = triangles[idx];
            const long double value = maximize_overlap(t[0], t[1], t[2]);

            const long double y = value - comp;
            const long double z = sum + y;
            comp = (z - sum) - y;
            sum = z;
        }

        partial[tid] = sum;
    };

    std::vector<std::thread> threads;
    threads.reserve(thread_count);
    for (unsigned int tid = 0; tid < thread_count; ++tid) {
        threads.emplace_back(worker, tid);
    }
    for (std::thread& th : threads) {
        th.join();
    }

    long double total = 0.0L;
    long double comp = 0.0L;
    for (long double part : partial) {
        const long double y = part - comp;
        const long double z = total + y;
        comp = (z - total) - y;
        total = z;
    }
    return total;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    const long double answer = solve();
    std::cout << std::fixed << std::setprecision(2) << std::round(answer * 100.0L) / 100.0L << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


if __name__ == "__main__":
    print(solve())

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;

public class Euler966 {

    static final double K_PI = 3.141592653589793238462643383279502884;
    static final double K_EPS = 1e-12;

    static class Vec2 {
        double x, y;

        Vec2(double x, double y) {
            this.x = x;
            this.y = y;
        }

        Vec2 add(Vec2 o) {
            return new Vec2(x + o.x, y + o.y);
        }

        Vec2 sub(Vec2 o) {
            return new Vec2(x - o.x, y - o.y);
        }

        Vec2 mul(double s) {
            return new Vec2(x * s, y * s);
        }
    }

    static double dot(Vec2 a, Vec2 b) {
        return a.x * b.x + a.y * b.y;
    }

    static double cross(Vec2 a, Vec2 b) {
        return a.x * b.y - a.y * b.x;
    }

    static double norm2(Vec2 v) {
        return dot(v, v);
    }

    static double dist(Vec2 a, Vec2 b) {
        return Math.sqrt(norm2(a.sub(b)));
    }

    static class Triangle {
        Vec2 a, b, c;
        double area, radius, maxSide;

        Triangle(Vec2 a, Vec2 b, Vec2 c, double area, double radius, double maxSide) {
            this.a = a;
            this.b = b;
            this.c = c;
            this.area = area;
            this.radius = radius;
            this.maxSide = maxSide;
        }
    }

    static Triangle makeTriangle(int s1, int s2, int s3) {
        int[] sides = { s1, s2, s3 };
        Arrays.sort(sides);

        double a = sides[0];
        double b = sides[1];
        double c = sides[2];

        Vec2 A = new Vec2(0.0, 0.0);
        Vec2 B = new Vec2(c, 0.0);
        double x = (b * b + c * c - a * a) / (2.0 * c);
        double y2 = b * b - x * x;
        if (y2 < 0.0 && y2 > -1e-18) {
            y2 = 0.0;
        }
        Vec2 C = new Vec2(x, Math.sqrt(Math.max(0.0, y2)));

        double area = Math.abs(cross(B.sub(A), C.sub(A))) * 0.5;
        double radius = Math.sqrt(area / K_PI);
        double maxSide = Math.max(Math.max(dist(A, B), dist(B, C)), dist(C, A));

        return new Triangle(A, B, C, area, radius, maxSide);
    }

    static boolean pointInTriangle(Triangle tri, Vec2 p) {
        double c1 = cross(tri.b.sub(tri.a), p.sub(tri.a));
        double c2 = cross(tri.c.sub(tri.b), p.sub(tri.b));
        double c3 = cross(tri.a.sub(tri.c), p.sub(tri.c));
        return c1 >= -1e-14 && c2 >= -1e-14 && c3 >= -1e-14;
    }

    static Vec2 nearestOnSegment(Vec2 p, Vec2 u, Vec2 v) {
        Vec2 uv = v.sub(u);
        double d = norm2(uv);
        if (d <= 0.0)
            return u;
        double t = dot(p.sub(u), uv) / d;
        if (t < 0.0)
            t = 0.0;
        else if (t > 1.0)
            t = 1.0;
        return u.add(uv.mul(t));
    }

    static Vec2 projectToTriangle(Triangle tri, Vec2 p) {
        if (pointInTriangle(tri, p))
            return p;

        Vec2[] candidates = {
                tri.a, tri.b, tri.c,
                nearestOnSegment(p, tri.a, tri.b),
                nearestOnSegment(p, tri.b, tri.c),
                nearestOnSegment(p, tri.c, tri.a),
        };

        Vec2 best = candidates[0];
        double bestD = norm2(p.sub(best));
        for (int i = 1; i < 6; ++i) {
            double cur = norm2(p.sub(candidates[i]));
            if (cur < bestD) {
                bestD = cur;
                best = candidates[i];
            }
        }
        return best;
    }

    static double segmentContribution(Vec2 a, Vec2 b, double r2) {
        Vec2[] pts = new Vec2[4];
        int n = 0;
        pts[n++] = a;

        Vec2 d = b.sub(a);
        double qa = dot(d, d);
        double qb = 2.0 * dot(a, d);
        double qc = dot(a, a) - r2;
        double disc = qb * qb - 4.0 * qa * qc;

        if (disc > 1e-18) {
            double sdisc = Math.sqrt(disc);
            double t1 = (-qb - sdisc) / (2.0 * qa);
            double t2 = (-qb + sdisc) / (2.0 * qa);
            if (t1 > t2) {
                double tmp = t1;
                t1 = t2;
                t2 = tmp;
            }

            if (t1 > K_EPS && t1 < 1.0 - K_EPS) {
                pts[n++] = a.add(d.mul(t1));
            }
            if (t2 > K_EPS && t2 < 1.0 - K_EPS && Math.abs(t2 - t1) > 1e-12) {
                pts[n++] = a.add(d.mul(t2));
            }
        }

        pts[n++] = b;

        double acc = 0.0;
        for (int i = 0; i + 1 < n; ++i) {
            Vec2 p = pts[i];
            Vec2 q = pts[i + 1];
            Vec2 mid = new Vec2((p.x + q.x) * 0.5, (p.y + q.y) * 0.5);

            if (norm2(mid) <= r2 + 1e-13) {
                acc += 0.5 * cross(p, q);
            } else {
                acc += 0.5 * r2 * Math.atan2(cross(p, q), dot(p, q));
            }
        }
        return acc;
    }

    static double overlapArea(Triangle tri, Vec2 center) {
        double r = tri.radius;
        double r2 = r * r;

        Vec2 p0 = tri.a.sub(center);
        Vec2 p1 = tri.b.sub(center);
        Vec2 p2 = tri.c.sub(center);

        double area = 0.0;
        area += segmentContribution(p0, p1, r2);
        area += segmentContribution(p1, p2, r2);
        area += segmentContribution(p2, p0, r2);
        return Math.abs(area);
    }

    static class EvalPoint implements Comparable<EvalPoint> {
        Vec2 p;
        double v;

        EvalPoint(Vec2 p, double v) {
            this.p = p;
            this.v = v;
        }

        @Override
        public int compareTo(EvalPoint o) {
            return Double.compare(o.v, this.v);
        }
    }

    static EvalPoint nelderMead(Triangle tri, Vec2 start, double step, int iterations) {
        List<EvalPoint> simplex = new ArrayList<>(Arrays.asList(
                new EvalPoint(projectToTriangle(tri, start), overlapArea(tri, projectToTriangle(tri, start))),
                new EvalPoint(projectToTriangle(tri, new Vec2(start.x + step, start.y)),
                        overlapArea(tri, projectToTriangle(tri, new Vec2(start.x + step, start.y)))),
                new EvalPoint(projectToTriangle(tri, new Vec2(start.x, start.y + step)),
                        overlapArea(tri, projectToTriangle(tri, new Vec2(start.x, start.y + step))))));

        for (int it = 0; it < iterations; ++it) {
            Collections.sort(simplex);

            Vec2 centroid = new Vec2(
                    (simplex.get(0).p.x + simplex.get(1).p.x) * 0.5,
                    (simplex.get(0).p.y + simplex.get(1).p.y) * 0.5);

            Vec2 reflected = new Vec2(
                    centroid.x + (centroid.x - simplex.get(2).p.x),
                    centroid.y + (centroid.y - simplex.get(2).p.y));
            Vec2 pr = projectToTriangle(tri, reflected);
            EvalPoint er = new EvalPoint(pr, overlapArea(tri, pr));

            if (er.v > simplex.get(0).v) {
                Vec2 expanded = new Vec2(
                        centroid.x + 2.0 * (er.p.x - centroid.x),
                        centroid.y + 2.0 * (er.p.y - centroid.y));
                Vec2 pe = projectToTriangle(tri, expanded);
                EvalPoint ee = new EvalPoint(pe, overlapArea(tri, pe));
                simplex.set(2, ee.v > er.v ? ee : er);
                continue;
            }

            if (er.v > simplex.get(1).v) {
                simplex.set(2, er);
                continue;
            }

            Vec2 contracted = new Vec2(
                    centroid.x + 0.5 * (simplex.get(2).p.x - centroid.x),
                    centroid.y + 0.5 * (simplex.get(2).p.y - centroid.y));
            Vec2 pc = projectToTriangle(tri, contracted);
            EvalPoint ec = new EvalPoint(pc, overlapArea(tri, pc));
            if (ec.v > simplex.get(2).v) {
                simplex.set(2, ec);
                continue;
            }

            Vec2 p1 = new Vec2(
                    simplex.get(0).p.x + 0.5 * (simplex.get(1).p.x - simplex.get(0).p.x),
                    simplex.get(0).p.y + 0.5 * (simplex.get(1).p.y - simplex.get(0).p.y));
            Vec2 pp1 = projectToTriangle(tri, p1);
            simplex.set(1, new EvalPoint(pp1, overlapArea(tri, pp1)));

            Vec2 p2 = new Vec2(
                    simplex.get(0).p.x + 0.5 * (simplex.get(2).p.x - simplex.get(0).p.x),
                    simplex.get(0).p.y + 0.5 * (simplex.get(2).p.y - simplex.get(0).p.y));
            Vec2 pp2 = projectToTriangle(tri, p2);
            simplex.set(2, new EvalPoint(pp2, overlapArea(tri, pp2)));
        }

        Collections.sort(simplex);
        return simplex.get(0);
    }

    static double maximizeOverlap(int s1, int s2, int s3) {
        Triangle tri = makeTriangle(s1, s2, s3);
        List<Vec2> seeds = new ArrayList<>();

        int grid = 5;
        for (int i = 0; i <= grid; ++i) {
            for (int j = 0; j + i <= grid; ++j) {
                double u = (double) i / grid;
                double v = (double) j / grid;
                double w = 1.0 - u - v;
                seeds.add(new Vec2(
                        u * tri.a.x + v * tri.b.x + w * tri.c.x,
                        u * tri.a.y + v * tri.b.y + w * tri.c.y));
            }
        }

        seeds.add(new Vec2(
                (tri.a.x + tri.b.x + tri.c.x) / 3.0,
                (tri.a.y + tri.b.y + tri.c.y) / 3.0));

        double sideA = dist(tri.b, tri.c);
        double sideB = dist(tri.c, tri.a);
        double sideC = dist(tri.a, tri.b);
        double ws = sideA + sideB + sideC;

        seeds.add(new Vec2(
                (sideA * tri.a.x + sideB * tri.b.x + sideC * tri.c.x) / ws,
                (sideA * tri.a.y + sideB * tri.b.y + sideC * tri.c.y) / ws));

        List<EvalPoint> evaluated = new ArrayList<>();
        for (Vec2 p : seeds) {
            Vec2 proj = projectToTriangle(tri, p);
            evaluated.add(new EvalPoint(proj, overlapArea(tri, proj)));
        }

        Collections.sort(evaluated);

        EvalPoint best = evaluated.get(0);
        int tries = Math.min(4, evaluated.size());
        double step = tri.maxSide * 0.18;

        for (int i = 0; i < tries; ++i) {
            EvalPoint cur = nelderMead(tri, evaluated.get(i).p, step, 90);
            if (cur.v > best.v) {
                best = cur;
            }
        }

        EvalPoint refined = nelderMead(tri, best.p, tri.maxSide * 0.05, 50);
        if (refined.v > best.v) {
            best = refined;
        }

        if (best.v < 0.0)
            return 0.0;
        if (best.v > tri.area)
            return tri.area;
        return best.v;
    }

    static boolean nearlyEqual(double x, double y, double eps) {
        return Math.abs(x - y) <= eps;
    }

    static List<int[]> enumerateTriangles() {
        List<int[]> all = new ArrayList<>();
        for (int a = 1; a <= 200; ++a) {
            for (int b = a; b <= 200 - a; ++b) {
                int maxC = 200 - a - b;
                for (int c = b; c <= maxC; ++c) {
                    if (a + b <= c)
                        continue;
                    all.add(new int[] { a, b, c });
                }
            }
        }
        return all;
    }

    public static String solve() {
        List<int[]> triangles = enumerateTriangles();

        double total = triangles.parallelStream()
                .mapToDouble(t -> maximizeOverlap(t[0], t[1], t[2]))
                .sum();

        total = Math.round(total * 100.0) / 100.0;
        return String.format(java.util.Locale.US, "%.2f", total);
    }

    public static void main(String[] args) {
        double i345 = maximizeOverlap(3, 4, 5);
        if (!nearlyEqual(i345, 4.593049, 2e-6)) {
            System.out.println("Validation failed");
            return;
        }

        double i346 = maximizeOverlap(3, 4, 6);
        if (!nearlyEqual(i346, 3.552564, 2e-6)) {
            System.out.println("Validation failed");
            return;
        }

        double iPerm = maximizeOverlap(5, 3, 4);
        if (!nearlyEqual(iPerm, i345, 1e-10)) {
            System.out.println("Validation failed");
            return;
        }

        System.out.println(solve());
    }
}