[C++] 34041번 회전체와 쿼리

https://www.acmicpc.net/problem/34041

---

 

25/10/18

 

 

미적분학 + 기하학 내용으로 풀 수 있는 다이아문제다. 내용이 그리 어렵진 않은 것 같다.

---

 

문제 접근 방식:

 

 

기본적으로 Pappus Centroid Theorem을 이용한다.

 

해당 정리에 대한 증명은 다음과 같다.

# Centroid of a Region
The centroid $(\bar{x}, \bar{y})$ of a region bounded by the graphs of the continuous functions $f$ and $g$ such that $f(x) \geq g(x)$ on the interval $[a, b]$, $a \leq x \leq b$ is given by $$ \bar{x} = \frac{1}{A} \int_{a}^{b} x(f(x) -g(x))dx $$

- 직관적인 설명
영역의 전체 넓이 $A$는 $$ \int_{a}^{b}[f(x) - g(x)]dx $$ 즉, 세로($f(x)-g(x)$)$\times$가로($dx$)의 모든 합임.
$x$축 방향으로의 얇은 세로 띠들의 면적합이 전체 면적임.
무게중심의 $x$좌표는 모든 얇은 띠들이 가진 $x$좌표의 가중 평균과 같음.
각 띠의 대표 $x$좌표는 띠의 $x$좌표이고, 가중치는 해당 띠의 면적($[f(x)-g(x)]\cdot dx$)임.
따라서, $$ \bar{x} = \frac{1}{A}\int_{a}^{b}x[f(x)-g(x)]dx $$

# Pappus Centroid Theorem
The volume $V$ of a solid of revolution generated by rotating a plane figure $F$ about an external axis $=$ (The area $A$ of $F$)$\cdot$(The distance $d$ traveled by the geometric centroid of $F$)

Proof)
Let the area $A$ is bounded by the two functions $y = f(x), y \geq 0$ and $y = g(x), f(x) \geq g(x)$, and bounded by the two lines $x=a\geq {0}$ and $x = b \geq a$.
Note that $$ A = \int_{a}^{b} dA = \int_{a}^{b} [f(x)-g(x)]dx $$
If this area is rotated about the $y$-axis, the volume generated can be calculated using the shell method,
Hence, $$ V = 2\pi \int_{a}^{b}x[f(x)-g(x)]dx $$ Hence, $V = 2\pi \bar{x}A$.

- 직관적인 이해
면적 $A$를 이루는 얇은 세로 띠를 보자. 이 얇은 세로 띠를 회전시켰다고 생각하면, 그 모양은 얇은 원통 껍질모양이 된다.
미소 부피 $dV$ = (해당 껍질의 넓이)$\times$(해당 껍질의 두께)가 되므로, $dV = 2\pi x[f(x)-g(x)] \times dx$

 

먼저, 해당 정리를 이용하기 위해서는 3개의 정보가 필요하다.

 

1. 분할된 다각형의 넓이 ($A$)

2. 분할된 다각형에 대한 무게중심 좌표 $(C_{x}, C_{y})$

3. 해당 좌표와 분할 직선 사이의 거리 $d$

 

회전체의 부피는 해당 정리에 의하여 $2\pi dA$가 된다.

 

먼저, 분할된 두 다각형의 넓이는 아래 공식(신발끈 공식)을 통해 $\mathcal{O}(N)$의 시간복잡도로 구할 수 있다.

 

$$A = \frac{1}{2}\sum_{i=0}^{n-1}(x_{i}y_{i+1} - x_{i+1}y_{i})$$

 

문제는 각 쿼리마다 이를 구하려고 한다면, 쿼리가 총 $100 \ 000$개가 주어지기 때문에 시간초과를 받는다는 점이다.

 

다각형의 넓이 공식을 보면 현재 점과 다음 점 사이의 관계식으로 모두 쪼개질 수 있기 때문에, 이 부분을 "누적합"으로 처리한다면 각 쿼리를 $\mathcal{O}(1)$의 시간 복잡도로 빠르게 처리할 수 있을 것이다.

 

예를 들어, 0번 점부터 3번 점까지의 관계는 기존 누적합으로 처리하고, 3번 점과 0번 점 사이의 관계는 따로 더해줌으로써 파란색 영역의 넓이를 구할 수 있다.

 

두번째로, 분할된 다각형에 대한 무게 중심 좌표를 구해야 한다. 해당 무게 중심 좌표는 다음과 같은 공식을 통해 구할 수 있다.

 

$$C_{x} = \frac{1}{6A}\sum_{i=0}^{n-1}(x_{i}+x_{i+1})(x_{i}y_{i+1} - x_{i+1}y_{i}), \quad C_{y}=\frac{1}{6A}\sum_{i=0}^{n-1}(y_{i} + y_{i+1})(x_{i}y_{i+1} - x_{i+1}y_{i})$$

 

해당 무게 중심 좌표가 왜 저렇게 나오는 지에 대한 직관적인 설명을 하면, 다음과 같다.

 

먼저, $N$-다각형을 $N-2$개의 삼각형으로 분할한다.

 

각각의 삼각형의 면적은 $\frac{1}{2}(x_{i}y_{i+1} - x_{i+1}y_{i})$이다.

 

각각의 삼각형의 무게중심 좌표는 $\frac{1}{3}(x_{i} + x_{i+1})$이다.

 

전체 다각형의 무게중심 좌표는 각각의 삼각형의 무게중심 좌표의 "가중 평균"인데, 면적에 대한 가중평균이므로 해당 공식이 나오게 된다.

 

이것도 넓이를 구하는 과정과 마찬가지로 각 쿼리마다 해당 값을 일일히 구하려고 한다면 무리이므로, 누적합을 통해서 $\mathcal{O}(1)$의 시간 복잡도로 구해준다.

 

마지막으로 분할된 다각형에 대한 무게 중심 좌표를 구했다면 점과 직선 사이의 거리 공식을 이용해 $d$를 구한다.

 

점과 직선 사이의 거리 공식은 다음과 같다.

$$d = \frac{|ax+by+c|}{\sqrt{ a^{2}+b^{2} }}$$

 

만약, 두 점 $(p_{x}, p_{y}), (q_{x}, q_{y})$를 이은 직선이 있다면, 직선의 계수에 해당하는 $a, b, c$의 값은 다음과 같다.

 

$$\begin{align} y &= \frac{\Delta y}{\Delta x}(x-p_{x}) + p_{y} \\ \rightarrow \Delta x \cdot y &= \Delta y(x-p_{x}) + p_{y}\cdot \Delta x \\ \rightarrow (q_{x}-p_{x})x &+ (p_{y}-q_{y})y + p_{y}(q_{x}-p_{x}) -p_{x}(q_{y}-p_{y}) = 0 \\ \rightarrow &\begin{cases} a = q_{x} - p_{x} \\ b = p_{y} - q_{y} \\ c = p_{y}(q_{x}-p_{x}) -p_{x}(q_{y}-p_{y}) \end{cases}  \end{align} $$

 

그렇게 각각의 다각형에 대해서 부피를 구해준 다음에 비교해서 각 쿼리를 $\mathcal{O}(1)$에 처리하여 구하면 된다.

 

유의할 점은, 넓이의 정확성을 위해 long double을 사용해야한다는 점이다.

---

아래는 내가 위의 접근 방식과 같이 작성한 C++ 코드이다. 더보기를 누르면 확인할 수 있다.

// 34041번 회전체와 쿼리
// 기하학, 수학, 누적 합
/*
접근 방법:
다각형의 영역의 넓이
무게중심 좌표의 x값, y값
모두 구할 때 현재 점과 다음 점 사이의 관계로 표현된다.
누적합으로 일정한 점까지 빠르게 구할 수 있다.
그리고나서 두 점은 따로 처리하면 될듯.
*/
#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

using namespace std;
#define fastio ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
#define endl '\n'
#define PI 3.14159265358979323846
#define EPS 1e-12L
typedef long double ld;


struct Point{
    ld x, y;
};

ld cross(Point& p, Point& q){
    return (p.x*q.y - q.x*p.y);
}
ld x_cross(Point& p, Point& q){
    return (p.x+q.x)*(p.x*q.y - q.x*p.y);
}
ld y_cross(Point& p, Point& q){
    return (p.y+q.y)*(p.x*q.y - q.x*p.y);
}
// return coefficient of ax+by+c = 0. (a, b, c)
vector<ld> make_line(Point& p, Point& q){
    return {q.y-p.y, p.x-q.x, p.y*(q.x-p.x)-p.x*(q.y-p.y)};
}
// return dist between line and point
ld dist_line_point(Point& LP1, Point& LP2, Point& p){
    vector<ld> coeff = make_line(LP1, LP2);
    ld a = coeff[0], b = coeff[1], c = coeff[2];
    ld ret = abs(a*p.x+b*p.y+c) / sqrt(a*a + b*b);
    return ret;
}
// return accumed summation vector
vector<ld> accum(vector<ld>& arr){
    vector<ld> ret(arr.size()+1, 0);
    for (size_t i = 0; i < arr.size(); ++i){
        ret[i+1] = (ret[i] + arr[i]);
    }
    return ret;
}

int main(void){
    fastio

    int N; cin >> N;
    vector<Point> points(N);
    for (int i = 0; i < N; ++i){
        cin >> points[i].x >> points[i].y;
    }
    
    vector<ld> CROSS(N), X_CROSS(N), Y_CROSS(N);
    for (int i = 0; i < N; ++i){
        CROSS[i] = cross(points[i], points[(i+1)%N]);
        X_CROSS[i] = x_cross(points[i], points[(i+1)%N]);
        Y_CROSS[i] = y_cross(points[i], points[(i+1)%N]);
    }
    vector<ld> aCROSS = accum(CROSS), aX_CROSS = accum(X_CROSS), aY_CROSS = accum(Y_CROSS);

    int Q; cin >> Q;
    int i, j;
    cout << fixed << setprecision(20);
    while (Q--){
        cin >> i >> j; i--; j--;
        if (i < j) swap(i, j); // i가 큰거, j가 작은거.
        ld dA1 = (aCROSS[i] - aCROSS[j] + cross(points[i], points[j]));
        ld dA2 = aCROSS[N] - dA1;
        if (fabsl(dA1) < EPS || fabsl(dA2) < EPS){
            cout << 0 << endl;
            continue;
        }
        ld x1 = (aX_CROSS[i] - aX_CROSS[j] + x_cross(points[i], points[j]));
        ld x2 = (aX_CROSS[N] - aX_CROSS[i] + aX_CROSS[j] + x_cross(points[j], points[i]));
        ld y1 = (aY_CROSS[i] - aY_CROSS[j] + y_cross(points[i], points[j]));
        ld y2 = (aY_CROSS[N] - aY_CROSS[i] + aY_CROSS[j] + y_cross(points[j], points[i]));
        Point p1 = {x1/(3*dA1), y1/(3*dA1)};
        Point p2 = {x2/(3*dA2), y2/(3*dA2)};

        if (dist_line_point(points[i], points[j], p1)*dA1 < dist_line_point(points[i], points[j], p2)*dA2){
            cout << PI*dist_line_point(points[i], points[j], p1)*dA1 << endl;
        } else {
            cout << PI*dist_line_point(points[i], points[j], p2)*dA2 << endl;
        }
    }
    return 0;
}