AOJ 1394: Fair Chocolate-Cutting

この記事は帰ってきた AOJ-ICPC Advent Calendar 2022 1 日目の記事です。

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1394

問題概要

n 頂点の凸多角形が与えられる。この多角形の辺上に端点があり、多角形の面積を 2 等分する線分の長さとしてありうる最大値と最小値を求めよ。

解法

多角形の頂点を反時計回りに P _ 0, P _ 1, \dots, P _ {n-1}、面積を S とします。また、 P _ N = P _ 0 であるとします。

線分の端点を A, B (ただし A は辺 P _ i P _ {i+1} 上、 B は辺 P _ j P _ {j+1} 上の点) とします。多角形 P _ {i+1} P _ {i+2} \dots P _ j の面積を S _ 1、多角形 P _ {j+1} P _ {j+2} \dots P _ i の面積を S _ 2 とおくと、明らかに S _ 1, S _ 2 \leq \frac{S}{2} が必要です。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行かどうかで場合分けをします。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行でない場合

直線 P _ i P _ {i+1} と直線 P _ j P _ {j+1} の交点を Q とします。このとき、 Q が半直線 P _ i P _ {i+1} 上にある場合のみ考えるとしてよいです。(そうでない場合、 i と j を逆にして考える)

QP _ {i+1} = a _ 1, QA = a, QP _ i = a _ 2, QP _ j = b _ 1, QB = b, QP _ {j+1} = b _ 2, \angle P _ {i+1} Q P _ j = \theta とおくと、 a, b の取りうる値の範囲は a _ 1 \leq a \leq a _ 2, b _ 1 \leq b \leq b _ 2 となります。また、多角形 AP _ {i+1} P _ {i+2} \dots P _ j B の面積は S _ 1 + \frac{1}{2} (ab - a _ 1 b _ 1) \sin \theta となります。これが \frac{1}{2}S と等しくなればよいので、式変形して

ab = a _ 1 b _ 1 + \frac{S - 2S _ 1}{\sin \theta}

となります。右辺を p とおくと、 a = \frac{p}{b}, b = \frac{p}{a} です。よって、 a _ \text{min} = \max(a _ 1, \frac{p}{b _ 2}), a _ \text{max} = \min(a _ 2, \frac{p}{b _ 1}) とおくと、 a の取りうる値の範囲は a _ \text{min} \leq a \leq a _ \text{max} となります。

線分 AB の長さは余弦定理より \sqrt{a ^ 2 + b ^ 2 - 2ab \cos \theta} です。よって a ^ 2 + b ^ 2 - 2ab \cos \theta の最大値と最小値を求めればよいですが、これは (a + b) ^ 2 - 2p (1 + \cos\theta) に等しいので、 a + b = a + \frac{p}{a} の最大値と最小値を考えればよいです。

a > 0 において、 a + \frac{p}{a} は下に凸な関数で a = \sqrt{p} で極小値を取ります。従って、最大値を与える a の候補は a _ \text{min}, a _ \text{max}、最小値を与える a の候補は a _ \text{min}, a _ \text{max}, \sqrt{p} となります。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行な場合

P _ {i+1} A = a, P _ {i+1} P _ i = a _ 0, P _ j B = b, P _ j P _ {j+1} = b _ 0 とおくと、と、 a, b の取りうる値の範囲は 0 \leq a \leq a _ 0, 0 \leq b \leq b _ 0 となります。また、直線 P _ i P _ {i +1}, P _ j P _ {j+1} の距離を d とおくと、多角形 AP _ {i+1} P _ {i+2} \dots P _ j B の面積は S _ 1 + \frac{1}{2} (a+b)d となります。これが \frac{1}{2}S と等しくなればよいので、式変形して

a+b = \frac{S - 2S _ 1}{d}

となります。右辺を s とおくと、 a = s - b, b = s - a です。よって、 a _ \text{min} = \max(0, s - b _ 0), a _ \text{max} = \min(a _ 0, s) とおくと、 a の取りうる値の範囲は a _ \text{min} \leq a \leq a _ \text{max} となります。

幾何的な考察から、最大値を与える a の候補は a _ \text{min}, a _ \text{max}、最小値の候補は a = a _ \text{min}, a _ \text{max} と AB \perp P _ i P _ {i+1} となるような a です。 AB \perp P _ i P _ {i+1} となるような a が存在するかは、内積の計算で判定することができます。

実装

judge.u-aizu.ac.jp

#include <bits/stdc++.h>
using namespace std;
const long double eps = 0.000001;
const long double INF = 100000;
struct point{
  long double x, y;
  point(){
  }
  point(long double x, long double y): x(x), y(y){
  }
  point operator +(point P){
    return point(x + P.x, y + P.y);
  }
  point operator -(point P){
    return point(x - P.x, y - P.y);
  }
  point operator *(long double k){
    return point(x * k, y * k);
  }
  point operator /(long double k){
    return point(x / k, y / k);
  }
};
long double abs(point P){
  return sqrt(P.x * P.x + P.y * P.y);
}
long double dist(point P, point Q){
  return abs(Q - P);
}
point unit(point P){
  return P / abs(P);
}
long double dot(point P, point Q){
  return P.x * Q.x + P.y * Q.y;
}
long double cross(point P, point Q){
  return P.x * Q.y - P.y * Q.x;
}
long double angle_cosine(point P, point Q){
  return dot(P, Q) / (abs(P) * abs(Q));
}
long double angle_sine(point P, point Q){
  return abs(cross(P, Q)) / (abs(P) * abs(Q));
}
long double angle_one_minus_cosine(point P, point Q){
  return 1 - angle_cosine(P, Q);
}
long double area(point A, point B, point C){
  return abs(cross(B - A, C - A)) / 2;
}
struct line{
  point A, B;
  line(){
  }
  line(point A, point B): A(A), B(B){
  }
};
point vec(line L){
  return L.B - L.A;
}
double point_line_distance(point P, line L){
  return abs(cross(P - L.A, vec(L))) / abs(vec(L));
}
point line_intersection(line L1, line L2){
  return L1.A + vec(L1) * cross(L2.A - L1.A, vec(L2)) / cross(vec(L1), vec(L2));
}
long double get_dist(long double a, long double b, long double one_minus_cost){
  return sqrt(max(pow(a - b, 2) + 2 * a * b * one_minus_cost, (long double) 0));
}
int main(){
  cout << fixed << setprecision(20);
  int n;
  cin >> n;
  vector<point> P(n + 1);
  for (int i = 0; i < n; i++){
    cin >> P[i].x >> P[i].y;
  }
  P[n] = P[0];
  long double S = 0;
  for (int i = 0; i < n; i++){
    S += cross(P[i], P[i + 1]) / 2;
  }
  long double mn = INF, mx = 0;
  for (int i = 0; i < n; i++){
    long double S1 = 0, S2 = S;
    for (int p = 1; p < n; p++){
      int j = (i + p) % n;
      S2 -= area(P[i], P[j], P[j + 1]);
      if (p >= 3){
        int j2 = (j - 1 + n) % n;
        S1 += area(P[j], P[j2], P[i + 1]);
      }
      if (S1 <= S / 2 + eps && S2 <= S / 2 + eps){
        long double C = cross(P[j + 1] - P[j], P[i] - P[i + 1]);
        if (C > eps){
          long double sint = angle_sine(P[j + 1] - P[j], P[i] - P[i + 1]);
          long double one_minus_cost = angle_one_minus_cosine(P[j + 1] - P[j], P[i] - P[i + 1]);
          point A = line_intersection(line(P[i], P[i + 1]), line(P[j], P[j + 1]));
          long double a1 = dist(A, P[i + 1]);
          long double a2 = dist(A, P[i]);
          long double b1 = dist(A, P[j]);
          long double b2 = dist(A, P[j + 1]);
          long double pr = a1 * b1 + (S - S1 * 2) / sint;
          long double xmn = max(a1, pr / max(b2, eps));
          long double xmx = min(a2, pr / max(b1, eps));
          if (xmn <= xmx + eps){
            long double d1 = get_dist(xmn, pr / xmn, one_minus_cost);
            mn = min(mn, d1);
            mx = max(mx, d1);
            long double d2 = get_dist(xmx, pr / xmx, one_minus_cost);
            mn = min(mn, d2);
            mx = max(mx, d2);
            long double x = sqrt(pr);
            if (xmn - eps <= x && x <= xmx + eps){
              long double d3 = get_dist(x, x, one_minus_cost);
              mn = min(mn, d3);
            }
          }
        } else if (C > -eps){
          long double a = dist(P[i], P[i + 1]);
          long double b = dist(P[j], P[j + 1]);
          long double h = point_line_distance(P[i + 1], line(P[j], P[j + 1]));
          long double s = (S - S1 * 2) / h;
          if (-eps <= s && s <= a + b + eps){
            long double xmn = max((long double) 0, s - b);
            long double xmx = min(a, s);
            point P1 = P[i + 1] + unit(P[i] - P[i + 1]) * xmn;
            point Q1 = P[j] + unit(P[j + 1] - P[j]) * (s - xmn);
            long double d1 = dist(P1, Q1);
            mn = min(mn, d1);
            mx = max(mx, d1);
            point P2 = P[i + 1] + unit(P[i] - P[i + 1]) * xmx;
            point Q2 = P[j] + unit(P[j + 1] - P[j]) * (s - xmx);
            long double d2 = dist(P2, Q2);
            mn = min(mn, d2);
            mx = max(mx, d2);
            if (dot(P[i] - P1, Q1 - P1) > -eps && dot(P[i] - P2, Q2 - P2) < eps){
              mn = min(mn, h);
            }
          }
        }
      }
    }
  }
  cout << mn << endl;
  cout << mx << endl;
}