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

• 問題: No.3675 偏光板
• ユーザー: TKTYI
• 提出日時: 2026-09-03 02:56:32
• 言語: C++23 (gcc 15.3.0 + boost 1.92.0)
• 結果: AC
• 実行時間: 93 ms / 2,000 ms
• コード長: 11,450 bytes
• コンパイル時間: 2,257 ms
• コンパイル使用メモリ: 257,600 KB
• 実行使用メモリ: 9,796 KB
• 最終ジャッジ日時: 2026-09-04 23:11:47
• 合計ジャッジ時間: 5,796 ms
• ジャッジサーバーID （参考情報）: judge1_0 / judge6_0

ファイルパターン	結果
sample	AC * 6
other	AC * 54

ソースコード

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

using namespace std;

struct Point2D
{
  long double x, y;
};
struct Circle
{
  long double x, y, r;
};
struct Point3D
{
  long double x, y, z;
  int id;

  Point3D operator-(const Point3D &p) const
  {
    return {x - p.x, y - p.y, z - p.z, id};
  }
};
struct Face
{
  int i, j, k;
  bool visible = true;
  vector<int> pts;
};

vector<Point2D> clip_polygon_halfplane(const vector<Point2D> &poly, long double a, long double b, long 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];
    long double d1 = a * p1.x + b * p1.y + c;
    long double d2 = a * p2.x + b * p2.y + c;

    if (d1 >= -1e-9L)
    {
      res.push_back(p1);
      if (d2 < -1e-9L)
      {
        long double t = d1 / (d1 - d2);
        t = max(0.0L, min(1.0L, t));
        res.push_back({p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y)});
      }
    }
    else if (d2 >= -1e-9L)
    {
      long double t = d1 / (d1 - d2);
      t = max(0.0L, min(1.0L, 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, long double X, long double Y)
{
  poly = clip_polygon_halfplane(poly, 1.0L, 0.0L, 0.0L);
  poly = clip_polygon_halfplane(poly, -1.0L, 0.0L, X);
  poly = clip_polygon_halfplane(poly, 0.0L, 1.0L, 0.0L);
  poly = clip_polygon_halfplane(poly, 0.0L, -1.0L, Y);
  return poly;
}

long double circle_polygon_area(long double cx, long double cy, long double R, vector<Point2D> poly)
{
  int n = poly.size();
  if (n < 3)
    return 0.0L;
  long double area = 0.0L;
  long double R2 = R * R;
  for (int i = 0; i < n; i++)
  {
    Point2D p1 = poly[i];
    Point2D p2 = poly[(i + 1) % n];
    long double x1 = p1.x - cx, y1 = p1.y - cy;
    long double x2 = p2.x - cx, y2 = p2.y - cy;
    long double dx = x2 - x1, dy = y2 - y1;
    long double a = dx * dx + dy * dy;
    long double b = 2.0L * (x1 * dx + y1 * dy);
    long double c = x1 * x1 + y1 * y1 - R2;
    vector<long double> ts = {0.0L, 1.0L};

    if (std::abs(a) > 1e-9L)
    {
      long double disc = b * b - 4.0L * a * c;
      if (disc > 0.0L)
      {
        long double sq = std::sqrt(disc);
        long double t1 = (-b - sq) / (2.0L * a);
        long double t2 = (-b + sq) / (2.0L * a);
        if (t1 > 0.0L && t1 < 1.0L)
          ts.push_back(t1);
        if (t2 > 0.0L && t2 < 1.0L)
          ts.push_back(t2);
      }
    }
    sort(ts.begin(), ts.end());
    for (size_t j = 0; j < ts.size() - 1; j++)
    {
      long double ta = ts[j], tb = ts[j + 1];
      if (tb - ta < 1e-12L)
        continue;
      Point2D sub_p1 = {x1 + ta * dx, y1 + ta * dy};
      Point2D sub_p2 = {x1 + tb * dx, y1 + tb * dy};
      long double mx = (sub_p1.x + sub_p2.x) / 2.0L;
      long double my = (sub_p1.y + sub_p2.y) / 2.0L;

      if (mx * mx + my * my <= R2 + 1e-9L)
      {
        area += 0.5L * (sub_p1.x * sub_p2.y - sub_p1.y * sub_p2.x);
      }
      else
      {
        long double ang1 = std::atan2(sub_p1.y, sub_p1.x);
        long double ang2 = std::atan2(sub_p2.y, sub_p2.x);
        long double dtheta = std::atan2(std::sin(ang2 - ang1), std::cos(ang2 - ang1));
        area += 0.5L * R2 * dtheta;
      }
    }
  }
  return std::abs(area);
}

Point3D cross(const Point3D &a, const Point3D &b)
{
  return {
      a.y * b.z - a.z * b.y,
      a.z * b.x - a.x * b.z,
      a.x * b.y - a.y * b.x,
      0};
}

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

long double orient3d(const Point3D &a, const Point3D &b, const Point3D &c, const Point3D &d)
{
  return dot(cross(b - a, c - a), d - a);
}

const long double EPS = 1e-9L;

vector<Face> get_convex_hull_3d(vector<Point3D> pts)
{
  int n = pts.size();
  if (n < 4)
    return {};

  int p0 = 0, p1 = 1, p2 = -1, p3 = -1;
  while (p1 < n && std::abs(pts[p1].x - pts[p0].x) < EPS &&
         std::abs(pts[p1].y - pts[p0].y) < EPS &&
         std::abs(pts[p1].z - pts[p0].z) < EPS)
  {
    p1++;
  }
  if (p1 == n)
    return {};

  for (int i = p1 + 1; i < n; i++)
  {
    Point3D crs = cross(pts[p1] - pts[p0], pts[i] - pts[p0]);
    if (std::abs(crs.x) > EPS || std::abs(crs.y) > EPS || std::abs(crs.z) > EPS)
    {
      p2 = i;
      break;
    }
  }
  if (p2 == -1)
    return {};

  for (int i = p2 + 1; i < n; i++)
  {
    if (std::abs(orient3d(pts[p0], pts[p1], pts[p2], pts[i])) > EPS)
    {
      p3 = i;
      break;
    }
  }
  if (p3 == -1)
    return {};

  vector<int> remaining;
  for (int i = 0; i < n; i++)
  {
    if (i != p0 && i != p1 && i != p2 && i != p3)
    {
      remaining.push_back(i);
    }
  }
  mt19937 rng(1337);
  shuffle(remaining.begin(), remaining.end(), rng);

  vector<Face> faces;
  auto add_face = [&](int a, int b, int c)
  {
    faces.push_back({a, b, c, true, {}});
  };

  if (orient3d(pts[p0], pts[p1], pts[p2], pts[p3]) > 0.0L)
  {
    add_face(p0, p2, p1);
    add_face(p0, p1, p3);
    add_face(p1, p2, p3);
    add_face(p2, p0, p3);
  }
  else
  {
    add_face(p0, p1, p2);
    add_face(p0, p3, p1);
    add_face(p1, p3, p2);
    add_face(p2, p3, p0);
  }

  for (int p_idx : remaining)
  {
    for (auto &f : faces)
    {
      if (orient3d(pts[f.i], pts[f.j], pts[f.k], pts[p_idx]) > EPS)
      {
        f.pts.push_back(p_idx);
        break;
      }
    }
  }

  for (int p_idx : remaining)
  {
    vector<int> vis_faces;
    for (int i = 0; i < (int)faces.size(); i++)
    {
      if (faces[i].visible && orient3d(pts[faces[i].i], pts[faces[i].j], pts[faces[i].k], pts[p_idx]) > EPS)
      {
        vis_faces.push_back(i);
      }
    }

    if (vis_faces.empty())
      continue;

    map<pair<int, int>, int> edge_cnt;
    for (int f_idx : vis_faces)
    {
      const auto &f = faces[f_idx];
      edge_cnt[{f.i, f.j}]++;
      edge_cnt[{f.j, f.k}]++;
      edge_cnt[{f.k, f.i}]++;
    }

    vector<pair<int, int>> horizon;
    for (int f_idx : vis_faces)
    {
      const auto &f = faces[f_idx];
      pair<int, int> edges[3] = {{f.i, f.j}, {f.j, f.k}, {f.k, f.i}};
      for (auto &e : edges)
      {
        if (edge_cnt.find({e.second, e.first}) == edge_cnt.end())
        {
          horizon.push_back(e);
        }
      }
    }

    vector<int> orphaned;
    for (int f_idx : vis_faces)
    {
      faces[f_idx].visible = false;
      for (int pt_i : faces[f_idx].pts)
      {
        if (pt_i != p_idx)
        {
          orphaned.push_back(pt_i);
        }
      }
    }
    sort(orphaned.begin(), orphaned.end());
    orphaned.erase(unique(orphaned.begin(), orphaned.end()), orphaned.end());

    int new_start_idx = faces.size();
    for (auto &e : horizon)
    {
      add_face(e.first, e.second, p_idx);
    }

    for (int pt_i : orphaned)
    {
      for (int i = new_start_idx; i < (int)faces.size(); i++)
      {
        if (orient3d(pts[faces[i].i], pts[faces[i].j], pts[faces[i].k], pts[pt_i]) > EPS)
        {
          faces[i].pts.push_back(pt_i);
          break;
        }
      }
    }
  }

  vector<Face> result;
  for (const auto &f : faces)
  {
    if (f.visible)
    {
      result.push_back(f);
    }
  }
  return result;
}

vector<Face> get_lower_hull(const vector<Point3D> &pts)
{
  auto all_faces = get_convex_hull_3d(pts);
  vector<Face> lower_faces;
  for (const auto &f : all_faces)
  {
    Point3D norm = cross(pts[f.j] - pts[f.i], pts[f.k] - pts[f.i]);
    if (norm.z < -EPS)
    {
      lower_faces.push_back(f);
    }
  }
  return lower_faces;
}

long double f(long double X, long double Y, vector<Circle> C)
{
  vector<Circle> unique_C;
  for (auto c : C)
  {
    bool dup = false;
    for (auto &uc : unique_C)
    {
      if (std::abs(c.x - uc.x) < 1e-7L && std::abs(c.y - uc.y) < 1e-7L)
      {
        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.0L;

  long double min_x = 0.0L, max_x = X, min_y = 0.0L, 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);
  }
  long double CX = (min_x + max_x) / 2.0L;
  long double CY = (min_y + max_y) / 2.0L;

  long double max_dim = max({max_x - min_x, max_y - min_y, X, Y, 100.0L});
  long double M = max_dim * 10.0L;

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

  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, 0});
  }

  vector<Face> lower = get_lower_hull(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];

    long double A = 2.0L * (p_j.x - p_i.x);
    long double B_coef = 2.0L * (p_j.y - p_i.y);
    long double E = pts[face.j].z - pts[face.i].z;

    long double C_coef = 2.0L * (p_k.x - p_i.x);
    long double D_coef = 2.0L * (p_k.y - p_i.y);
    long double F = pts[face.k].z - pts[face.i].z;

    long double det = A * D_coef - B_coef * C_coef;
    if (std::abs(det) < 1e-9L)
      continue;

    long double cx = (E * D_coef - B_coef * F) / det;
    long 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);
  }

  long double area = 0.0L;
  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 (std::hypot(c.x - u.x, c.y - u.y) < 1e-6L)
        {
          dup = true;
          break;
        }
      }
      if (!dup)
        unique_centers.push_back(c);
    }
    if (unique_centers.size() < 3)
      continue;

    long double gx = 0.0L, gy = 0.0L;
    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 std::atan2(a.y - gy, a.x - gx) < std::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()
{
  ios::sync_with_stdio(false);
  cin.tie(nullptr);

  long double X, Y;
  int N;
  if (!(cin >> X >> Y >> N))
    return 0;
  vector<Circle> v, h;
  mt19937 rnd(1337);
  uniform_real_distribution<long double> dist(-1e-8L, 1e-8L);

  for (int i = 0; i < N; i++)
  {
    long double x, y, r;
    char d;
    cin >> x >> y >> r >> d;

    x += dist(rnd);
    y += dist(rnd);
    r += dist(rnd);

    (d == 'V' ? v : h).push_back({x, y, r});
  }

  long double ans = X * Y;
  ans -= f(X, Y, v) * 0.5L;
  ans -= f(X, Y, h) * 0.5L;

  cout << fixed << setprecision(20) << (double)ans << "\n";
  return 0;
}