#1182146 (C++23) No.3675 偏光板

• 問題: No.3675 偏光板
• ユーザー: TKTYI
• 提出日時: 2026-09-03 01:15:07
• 言語: C++23 (gcc 15.3.0 + boost 1.92.0)
• 結果: TLE
• 実行時間: -
• コード長: 7,715 bytes
• コンパイル時間: 1,715 ms
• コンパイル使用メモリ: 226,824 KB
• 実行使用メモリ: 9,888 KB
• 最終ジャッジ日時: 2026-09-04 23:10:47
• 合計ジャッジ時間: 7,763 ms
• ジャッジサーバーID （参考情報）: judge2_0 / judge1_1

ファイルパターン	結果
sample	AC * 6
other	AC * 13 TLE * 1 -- * 40

ソースコード

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <random>

using namespace std;

struct Point2D { double x, y; };
struct Point3D { double x, y, z; };
struct Circle { double x, y, r; };
struct Face { int i, j, k; };

vector<Point2D> clip_polygon_halfplane(const vector<Point2D>& poly, double a, double b, double c) {
  vector<Point2D> res;
  int n = poly.size();
  if (n == 0) return res;

  for (int i = 0; i < n; i++) {
    Point2D p1 = poly[i];
    Point2D p2 = poly[(i + 1) % n];
    double d1 = a * p1.x + b * p1.y + c;
    double d2 = a * p2.x + b * p2.y + c;

    if (d1 >= -1e-9) {
      res.push_back(p1);
      if (d2 < -1e-9) {
        double t = d1 / (d1 - d2);
        t = max(0.0, min(1.0, t));
        res.push_back({p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y)});
      }
    } else if (d2 >= -1e-9) {
      double t = d1 / (d1 - d2);
      t = max(0.0, min(1.0, t));
      res.push_back({p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y)});
    }
  }
  return res;
}

vector<Point2D> clip_to_rect(vector<Point2D> poly, double X, double Y) {
  poly = clip_polygon_halfplane(poly, 1.0, 0.0, 0.0);
  poly = clip_polygon_halfplane(poly, -1.0, 0.0, X);
  poly = clip_polygon_halfplane(poly, 0.0, 1.0, 0.0);
  poly = clip_polygon_halfplane(poly, 0.0, -1.0, Y);
  return poly;
}

double circle_polygon_area(double cx, double cy, double R, vector<Point2D> poly) {
  int n = poly.size();
  if (n < 3) return 0.0;
  double area = 0.0;
  double R2 = R * R;
  for (int i = 0; i < n; i++) {
    Point2D p1 = poly[i];
    Point2D p2 = poly[(i + 1) % n];
    double x1 = p1.x - cx, y1 = p1.y - cy;
    double x2 = p2.x - cx, y2 = p2.y - cy;
    double dx = x2 - x1, dy = y2 - y1;
    double a = dx * dx + dy * dy;
    double b = 2.0 * (x1 * dx + y1 * dy);
    double c = x1 * x1 + y1 * y1 - R2;
    vector<double> ts = {0.0, 1.0};
    if (abs(a) > 1e-12) {
      double disc = b * b - 4 * a * c;
      if (disc > 0) {
        double sq = sqrt(disc);
        double t1 = (-b - sq) / (2 * a);
        double t2 = (-b + sq) / (2 * a);
        if (t1 > 0 && t1 < 1) ts.push_back(t1);
        if (t2 > 0 && t2 < 1) ts.push_back(t2);
      }
    }
    sort(ts.begin(), ts.end());
    for (size_t j = 0; j < ts.size() - 1; j++) {
      double ta = ts[j], tb = ts[j + 1];
      if (tb - ta < 1e-12) continue;
      Point2D sub_p1 = {x1 + ta * dx, y1 + ta * dy};
      Point2D sub_p2 = {x1 + tb * dx, y1 + tb * dy};
      double mx = (sub_p1.x + sub_p2.x) / 2.0;
      double my = (sub_p1.y + sub_p2.y) / 2.0;
      if (mx * mx + my * my <= R2 + 1e-9) {
        area += 0.5 * (sub_p1.x * sub_p2.y - sub_p1.y * sub_p2.x);
      } else {
        double ang1 = atan2(sub_p1.y, sub_p1.x);
        double ang2 = atan2(sub_p2.y, sub_p2.x);
        double dtheta = atan2(sin(ang2 - ang1), cos(ang2 - ang1));
        area += 0.5 * R2 * dtheta;
      }
    }
  }
  return abs(area);
}

// O(N^4) Lower Hull 計算 (精度確保のため内部は long double)
vector<Face> get_lower_hull_naive(const vector<Point3D>& pts) {
  vector<Face> lower_faces;
  int n = pts.size();
  for (int i = 0; i < n; i++) {
    for (int j = i + 1; j < n; j++) {
      for (int k = j + 1; k < n; k++) {
        long double v1x = pts[j].x - pts[i].x, v1y = pts[j].y - pts[i].y, v1z = pts[j].z - pts[i].z;
        long double v2x = pts[k].x - pts[i].x, v2y = pts[k].y - pts[i].y, v2z = pts[k].z - pts[i].z;
        long double nx = v1y * v2z - v1z * v2y;
        long double ny = v1z * v2x - v1x * v2z;
        long double nz = v1x * v2y - v1y * v2x;
        
        if (abs((double)nz) < 1e-9) continue;
        
        int a = i, b = j, c = k;
        if (nz < 0) {
          swap(b, c);
          nx = -nx; ny = -ny; nz = -nz;
        }
        
        long double norm = sqrt(nx*nx + ny*ny + nz*nz);
        bool is_face = true;
        for (int m = 0; m < n; m++) {
          if (m == a || m == b || m == c) continue;
          long double vx = pts[m].x - pts[a].x;
          long double vy = pts[m].y - pts[a].y;
          long double vz = pts[m].z - pts[a].z;
          long double dot = nx * vx + ny * vy + nz * vz;
          
          if (dot < -1e-7 * norm) {
            is_face = false;
            break;
          }
        }
        if (is_face) lower_faces.push_back({a, b, c});
      }
    }
  }
  return lower_faces;
}

double f(int X, int Y, vector<Circle> C) {
  vector<Circle> unique_C;
  for (auto c : C) {
    bool dup = false;
    for (auto& uc : unique_C) {
      if (abs(c.x - uc.x) < 1e-7 && abs(c.y - uc.y) < 1e-7) {
        uc.r = max(uc.r, c.r);
        dup = true;
        break;
      }
    }
    if (!dup) unique_C.push_back(c);
  }
  C = unique_C;

  int N = C.size();
  if (N == 0) return 0.0;

  double min_x = 0, max_x = X, min_y = 0, max_y = Y;
  for (auto c : C) {
    min_x = min(min_x, c.x); max_x = max(max_x, c.x);
    min_y = min(min_y, c.y); max_y = max(max_y, c.y);
  }
  double CX = (min_x + max_x) / 2.0;
  double CY = (min_y + max_y) / 2.0;
  double M = max({max_x - min_x, max_y - min_y, (double)X, (double)Y, 100.0}) * 10.0;

  vector<Circle> C_ext = C;
  C_ext.push_back({CX - M, CY - M, 0.0});
  C_ext.push_back({CX + M, CY - M, 0.0});
  C_ext.push_back({CX + M, CY + M, 0.0});
  C_ext.push_back({CX - M, CY + M, 0.0});

  vector<Point3D> pts;
  for (auto c : C_ext) {
    pts.push_back({c.x, c.y, c.x * c.x + c.y * c.y - c.r * c.r});
  }

  vector<Face> lower = get_lower_hull_naive(pts);
  vector<vector<Point2D>> circle_to_centers(C_ext.size());

  for (auto& face : lower) {
    Circle p_i = C_ext[face.i], p_j = C_ext[face.j], p_k = C_ext[face.k];
    
    double A = 2 * (p_j.x - p_i.x);
    double B_coef = 2 * (p_j.y - p_i.y);
    double E = pts[face.j].z - pts[face.i].z;
    
    double C_coef = 2 * (p_k.x - p_i.x);
    double D_coef = 2 * (p_k.y - p_i.y);
    double F = pts[face.k].z - pts[face.i].z;
    
    double det = A * D_coef - B_coef * C_coef;
    if (abs(det) < 1e-12) continue;
    
    double cx = (E * D_coef - B_coef * F) / det;
    double cy = (A * F - E * C_coef) / det;
    Point2D center = {cx, cy};

    circle_to_centers[face.i].push_back(center);
    circle_to_centers[face.j].push_back(center);
    circle_to_centers[face.k].push_back(center);
  }

  double area = 0.0;
  for (int i = 0; i < N; i++) {
    auto centers = circle_to_centers[i];
    if (centers.size() < 3) continue;

    vector<Point2D> unique_centers;
    for (auto c : centers) {
      bool dup = false;
      for (auto u : unique_centers) {
        if (hypot(c.x - u.x, c.y - u.y) < 1e-6) { dup = true; break; }
      }
      if (!dup) unique_centers.push_back(c);
    }
    if (unique_centers.size() < 3) continue;

    double gx = 0, gy = 0;
    for (auto p : unique_centers) gx += p.x, gy += p.y;
    gx /= unique_centers.size(); gy /= unique_centers.size();
    sort(unique_centers.begin(), unique_centers.end(), [&](const Point2D& a, const Point2D& b) {
      return atan2(a.y - gy, a.x - gx) < atan2(b.y - gy, b.x - gx);
    });

    vector<Point2D> clipped_poly = clip_to_rect(unique_centers, X, Y);

    if (clipped_poly.size() >= 3) {
      area += circle_polygon_area(C[i].x, C[i].y, C[i].r, clipped_poly);
    }
  }
  return area;
}

int main() {
  int X, Y, N;
  if (!(cin >> X >> Y >> N)) return 0;
  vector<Circle> v, h;
  mt19937 rnd(1337);
  uniform_real_distribution<double> dist(-1e-9, 1e-9);
  for (int i = 0; i < N; i++) {
    double x, y, r;
    char d;
    cin >> x >> y >> r >> d;
  
    x += dist(rnd);
    y += dist(rnd);
    r += abs(dist(rnd));
    
    (d == 'V' ? v : h).push_back({x, y, r});
  }
  double ans = (double) X * Y;
  ans -= f(X, Y, v) * 0.5;
  ans -= f(X, Y, h) * 0.5;
  printf("%.20lf\n", ans);
  return 0;
}