euler/966.cpp

/* Generated by GPT-5 on 2025-11-05 after a few attempts and telling previous run was wrong. Tested on Ubuntu 25.04. Runtime on Lenovo ThinkPad P14s: 0m12.971s */

// sum_I_fixed.cpp
// Compile: g++ -O3 -std=c++17 sum_I_fixed.cpp -o sum_I_fixed
// Run: ./sum_I_fixed

#include <bits/stdc++.h>
using namespace std;
struct P { double x,y; P():x(0),y(0){} P(double _x,double _y):x(_x),y(_y){} };
static inline P operator+(const P&a,const P&b){return P(a.x+b.x,a.y+b.y);}
static inline P operator-(const P&a,const P&b){return P(a.x-b.x,a.y-b.y);}
static inline P operator*(const P&a,double t){return P(a.x*t,a.y*t);}
static inline double dot(const P&a,const P&b){return a.x*b.x + a.y*b.y;}
static inline double cross(const P&a,const P&b){return a.x*b.y - a.y*b.x;}
static inline double norm2(const P&a){return dot(a,a);}
static inline double norm(const P&a){return sqrt(max(0.0, norm2(a)));}

const double PI = acos(-1.0);

// Solve quadratic for segment-circle intersections when circle is at origin.
// Return t values clipped to [0,1], sorted, unique.
vector<double> segment_circle_intersections(const P &A, const P &B, double r) {
    P d = B - A;
    double a = dot(d,d);
    vector<double> ts;
    if (a < 1e-18) return ts; // degenerate segment
    double b = 2 * dot(A, d);
    double c = dot(A,A) - r*r;
    double disc = b*b - 4*a*c;
    if (disc < -1e-12) return ts;
    if (disc < 0) disc = 0;
    double sdisc = sqrt(disc);
    double t1 = (-b - sdisc) / (2*a);
    double t2 = (-b + sdisc) / (2*a);
    if (t1 >= -1e-12 && t1 <= 1+1e-12) ts.push_back(max(0.0,min(1.0,t1)));
    if (t2 >= -1e-12 && t2 <= 1+1e-12) {
        double t2c = max(0.0,min(1.0,t2));
        if (ts.empty() || fabs(t2c - ts.back()) > 1e-12) ts.push_back(t2c);
    }
    sort(ts.begin(), ts.end());
    return ts;
}

// signed angle from a to b
static inline double angle_between(const P &a, const P &b) {
    return atan2(cross(a,b), dot(a,b));
}

// contribution for oriented edge A->B when circle centered at origin
double area_contrib_edge(const P &A, const P &B, double r) {
    double ra = norm(A), rb = norm(B);
    if (ra <= r + 1e-12 && rb <= r + 1e-12) {
        // fully inside circle: triangle area
        return 0.5 * cross(A, B);
    }
    // split by intersection points
    vector<P> pts;
    pts.push_back(A);
    auto ts = segment_circle_intersections(A, B, r);
    for (double t : ts) {
        P ip = A + (B - A) * t;
        if (norm(ip - pts.back()) > 1e-12) pts.push_back(ip);
    }
    if (norm(pts.back() - B) > 1e-12) pts.push_back(B);

    double sum = 0.0;
    for (size_t i = 0; i + 1 < pts.size(); ++i) {
        P u = pts[i], v = pts[i+1];
        P mid = P(0.5*(u.x+v.x), 0.5*(u.y+v.y));
        if (norm(mid) <= r + 1e-12) {
            sum += 0.5 * cross(u, v);
        } else {
            double ang = angle_between(u, v);
            sum += 0.5 * r * r * ang;
        }
    }
    return sum;
}

// polygon (triangle) ∩ circle area; circle center 'ctr', radius r
double polygon_circle_intersection_area(const vector<P> &poly, const P &ctr, double r) {
    double total = 0.0;
    int n = (int)poly.size();
    for (int i=0;i<n;++i) {
        P a = poly[i] - ctr;
        P b = poly[(i+1)%n] - ctr;
        total += area_contrib_edge(a,b,r);
    }
    return fabs(total);
}

void triangle_coords_from_sides(int a,int b,int c, P &A, P &B, P &C) {
    A = P(0.0,0.0);
    B = P((double)c,0.0);
    double x = ((double)b*b + (double)c*c - (double)a*a) / (2.0 * c);
    double y2 = (double)b*b - x*x;
    if (y2 < 0 && y2 > -1e-12) y2 = 0;
    double y = sqrt(max(0.0, y2));
    C = P(x,y);
}
double triangle_area(const P &A,const P &B,const P &C) { return fabs(0.5 * cross(B - A, C - A)); }
P centroid(const P&A,const P&B,const P&C){ return P((A.x+B.x+C.x)/3.0, (A.y+B.y+C.y)/3.0); }
P incenter(const P&A,const P&B,const P&C, double a,double b,double c) {
    double per = a+b+c;
    return P( (a*A.x + b*B.x + c*C.x)/per, (a*A.y + b*B.y + c*C.y)/per );
}
bool circumcenter(const P&A,const P&B,const P&C, P &O) {
    double d = 2*(A.x*(B.y - C.y) + B.x*(C.y - A.y) + C.x*(A.y - B.y));
    if (fabs(d) < 1e-15) return false;
    double ux = ( (A.x*A.x + A.y*A.y)*(B.y - C.y) + (B.x*B.x + B.y*B.y)*(C.y - A.y) + (C.x*C.x + C.y*C.y)*(A.y - B.y) ) / d;
    double uy = ( (A.x*A.x + A.y*A.y)*(C.x - B.x) + (B.x*B.x + B.y*B.y)*(A.x - C.x) + (C.x*C.x + C.y*C.y)*(B.x - A.x) ) / d;
    O = P(ux, uy);
    return true;
}

// stronger multi-start optimizer
double maximize_intersection_for_triangle(const P &A, const P &B, const P &C, double r) {
    vector<P> poly = {A,B,C};
    double side_a = norm(B - C);
    double side_b = norm(A - C);
    double side_c = norm(A - B);

    vector<P> starts;
    // geometric starts
    starts.push_back( centroid(A,B,C) );
    starts.push_back( incenter(A,B,C, side_a, side_b, side_c) );
    P cc; if (circumcenter(A,B,C,cc)) starts.push_back(cc);
    starts.push_back(A); starts.push_back(B); starts.push_back(C);
    starts.push_back(P(0.5*(A.x+B.x), 0.5*(A.y+B.y)));
    starts.push_back(P(0.5*(B.x+C.x), 0.5*(B.y+C.y)));
    starts.push_back(P(0.5*(C.x+A.x), 0.5*(C.y+A.y)));

    // bounding box expanded
    double xmin = min({A.x,B.x,C.x}) - 2.0*r;
    double xmax = max({A.x,B.x,C.x}) + 2.0*r;
    double ymin = min({A.y,B.y,C.y}) - 2.0*r;
    double ymax = max({A.y,B.y,C.y}) + 2.0*r;

    // deterministic grid starts (5x5)
    int G = 5;
    for (int i=0;i<G;++i) for (int j=0;j<G;++j) {
        double fx = (i + 0.5) / G;
        double fy = (j + 0.5) / G;
        starts.emplace_back( xmin + fx*(xmax-xmin), ymin + fy*(ymax-ymin) );
    }

    // radial starts from centroid (distances up to 2r, directions 12)
    P cen = centroid(A,B,C);
    for (int k=0;k<12;++k) {
        double ang = 2*PI*k/12.0;
        for (int t = 0; t <= 3; ++t) {
            double d = (t/3.0) * 2.0 * r; // 0..2r
            starts.emplace_back( cen.x + d*cos(ang), cen.y + d*sin(ang) );
        }
    }

    // many deterministic pseudo-random starts (LCG)
    uint64_t seed = 123456789ULL;
    for (int k=0;k<200;++k) {
        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
        double u = ((seed >>  0) & 0xffffffffULL) / double(0xffffffffULL);
        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
        double v = ((seed >> 16) & 0xffffffffULL) / double(0xffffffffULL);
        starts.emplace_back( xmin + u*(xmax-xmin), ymin + v*(ymax-ymin) );
    }

    auto eval = [&](const P &ctr)->double {
        return polygon_circle_intersection_area(poly, ctr, r);
    };

    double bestVal = 0.0;
    P bestCtr = cen;

    // evaluate starts
    for (const P &s0 : starts) {
        P s = s0;
        s.x = min(max(s.x, xmin), xmax);
        s.y = min(max(s.y, ymin), ymax);
        double v = eval(s);
        if (v > bestVal) { bestVal = v; bestCtr = s; }
    }

    // local pattern search from best found; shrink step
    double step = max({xmax - xmin, ymax - ymin, r});
    if (step < 1e-6) step = 1.0;
    const vector<P> dirs = { P(0,0), P(1,0), P(-1,0), P(0,1), P(0,-1),
                             P(1,1), P(-1,1), P(1,-1), P(-1,-1) };
    while (step > 1e-7) {
        bool improved = false;
        for (const P &d : dirs) {
            P cand = P(bestCtr.x + d.x * step, bestCtr.y + d.y * step);
            cand.x = min(max(cand.x, xmin), xmax);
            cand.y = min(max(cand.y, ymin), ymax);
            double v = eval(cand);
            if (v > bestVal + 1e-12) {
                bestVal = v; bestCtr = cand; improved = true;
            }
        }
        if (!improved) step *= 0.5;
    }

    return bestVal;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    const int LIMIT = 200;
    double total = 0.0;
    long long count = 0;
    long long printed = 0;

    for (int a=1;a<=LIMIT;++a){
        for (int b=a;b<=LIMIT;++b){
            int maxc_by_triangle = a + b - 1;
            int maxc_by_sum = LIMIT - a - b;
            int maxc = min(maxc_by_triangle, maxc_by_sum);
            if (maxc < b) continue;
            for (int c=b;c<=maxc;++c){
                P A,B,C;
                triangle_coords_from_sides(a,b,c,A,B,C);
                double tri_area = triangle_area(A,B,C);
                if (tri_area <= 0) continue;
                double r = sqrt(tri_area / PI);
                double I = maximize_intersection_for_triangle(A,B,C, r);
                total += I;
                ++count;
                if ((++printed) % 5000 == 0) {
                    cerr << "processed " << printed << " triangles (so far) ...\n";
                }
            }
        }
    }
    cout.setf(std::ios::fixed);
    cout<<setprecision(2)<<total<<"\n";
    return 0;
}