[C++] 24558번 Downsizing

 

 

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

---

 

25/10/21

 

 

반전 기하학의 내용을 담은 문제로, 나는 그린 정리의 내용을 활용하여 문제를 해결했다.

---

 

문제 접근 방식:

 

 

문제를 요약하면, 어떤 원과 그 원 내부에 있지 않은 다각형의 좌표들이 주어질 때, 다각형을 원에 대해서 반전시킨 넓이를 구하는 것이 목적이다.

 

편의 상 원의 중심의 좌표를 $(0, 0)$으로 옮기고, 다각형의 모든 점의 좌표 또한, 해당 이동에 맞게 전부 다 옮겨진 상태라고 가정하자.

 

원의 반지름을 $R$이라고 하자.

 

점 $P = (x, y)$의 Inversion $P' = (x', y')$은 닮음에 의하여 $(kx, ky)$가 된다.

 

이때, $k$의 값을 구하기 위해 $\text{length}(\overline{OP})\cdot\text{length}(\overline{OP'}) = R^2$을 이용하자.

 

$$\begin{align} k\sqrt{ x^{2} + y^{2}}\cdot \sqrt{ x^{2}+y^{2} } &= R^{2} \\ \rightarrow k &= \frac{R^{2}}{x^{2}+y^{2}} \end{align}$$

 

즉, 점 $(x, y)$가 $(x', y')$으로 이동하는 변환 $f(x, y)$는 다음과 같이 정의된다.

$$f(x, y) = \frac{R^{2}}{x^{2}+y^{2}}(x,y)$$

 

따라서, 변환 후의 미소 면적 $dA'$를 구하기 위해 자코비안 $\mathbf{J}$를 구해보자.

$$x' = \frac{R^{2}x}{x^{2}+y^{2}}, \quad y' = \frac{R^{2}y}{x^{2}+y^{2}}$$

 

$$\begin{align} \frac{ \partial x' }{ \partial x } &= R^{2} \frac{ \partial }{ \partial x } \left( x\cdot(x^{2}+y^{2})^{-1} \right) \\ &= R^{2}\left( (x^{2}+y^{2})^{-1}-2x^{2}\cdot(x^{2}+y^{2})^{-2} \right) \\ &= \frac{R^{2}}{(x^{2}+y^{2})^{2}} (y^{2}-x^{2}) \end{align}$$

 

$$\begin{align} \frac{ \partial x' }{ \partial y } &= R^{2} \frac{ \partial }{ \partial y } \left( x\cdot(x^{2}+y^{2})^{-1} \right) \\ &= -\frac{2R^{2}xy}{(x^{2}+y^{2})^{2}} \\ &= \frac{R^{2}}{(x^{2}+y^{2})^{2}}\cdot(-2xy) \end{align}$$

 

$$\frac{ \partial y' }{ \partial x } =  \frac{R^{2}}{(x^{2}+y^{2})^{2}}\cdot(-2xy)$$

 

$$\frac{ \partial y' }{ \partial y } =  \frac{R^{2}}{(x^{2}+y^{2})^{2}} (x^{2}-y^{2})$$

 

따라서 $\det(\mathbf{J})$의 값은,

$$\begin{align} \det(\mathbf{J}) &= \begin{vmatrix} \frac{ \partial x' }{ \partial x } & \frac{ \partial x' }{ \partial y } \\ \frac{ \partial y' }{ \partial x } & \frac{ \partial y' }{ \partial y } \\ \end{vmatrix} \\ &= \left(\frac{R^{2}}{(x^{2}+y^{2})^{2}}\right)^{2}\begin{vmatrix} y^{2}-x^{2} & -2xy \\ -2xy & x^{2}-y^{2} \end{vmatrix} \\ &= \left(\frac{R^{2}}{(x^{2}+y^{2})^{2}}\right)^{2}(-4x^{2}y^{2}-(x^{2}-y^{2})^{2}) \\ &= \left(\frac{R^{2}}{(x^{2}+y^{2})^{2}}\right)^{2}(-x^{4}-2x^{2}y^{2}-y^{4}) \\ &= -\left(\frac{R^{2}}{(x^{2}+y^{2})^{2}}\right)^{2}(x^{2}+y^{2})^{2} \\ &= -\frac{R^{4}}{(x^{2}+y^{2})^{2}} = -\frac{R^{4}}{r^4} \end{align}$$

 

즉, 미소 면적은 $dA'$은 $$\begin{align} dA' &= |\det(\mathbf{J})|dA \\ &= \frac{R^{4}}{r^{4}} dA \end{align}$$가 된다.

 

이를 직관적으로 이해해보자.

 

Inversion은 각 점을 원의 중심으로부터 뻗어나가는 방향(방사형 방향)으로만 늘리거나 줄인다.

 

이때의 Scale Factor $\lambda_{r}=\frac{R^{2}}{r^{2}}$이다.

 

평면에서의 면적은 방사형 방향과 수직 방향의 곱으로 이루어져있고, 둘 다 같은 Scale Factor를 가지므로 $$\lvert \lambda_{r} \rvert\times \lvert \lambda_{\perp} \rvert = \frac{R^{4}}{r^{4}}$$가 된다.

 

따라서, 원 밖에 있는 어떤 영역 $D$가 Inversion을 거치면 다음과 같이 변화한다.

 

$$dA' = \frac{R^{4}}{r^{4}}dA \rightarrow A' = \iint_{D} \frac{R^{4}}{(x^{2}+y^{2})^{2}}dA$$

 

이제 다각형에 대해서 면적분을 진행해야하는데, 면적분을 하기 힘드므로, 그린 정리를 사용하여 면적분을 선적분으로 바꾸자.

 

Green's Theorem)
Let $C$ be a positively oriented (counterclockwise), piecewise smooth, simple closed curve in the plane, and let $D$ be the region enclosed by $C$.
If $P(x, y)$ and $Q(x, y)$ have continuous partial derivatives on an open region containing $D$, then $$ \oint_{C}(Pdx+Qdy) = \iint_{D} \left( \frac{ \partial Q }{ \partial x } -\frac{ \partial P }{ \partial y } \right)dxdy $$

 

결국 그린 정리를 활용하기 위해서는 $$\frac{ \partial Q }{ \partial x } -\frac{ \partial P }{ \partial y } =\frac{R^{4}}{(x^{2}+y^{2})^{2}}$$가 되도록 $Q$와 $P$를 적절하게 잡아야 한다.

 

이때, 두 함수 $P$와 $Q$가 같은 함수에서 나온 형태라면 보기가 좋을 것이다. 즉, 라플라시안의 형태로 관찰해보자.

 

즉, 어떤 함수 $\psi(x, y)$가 존재해서, $Q = \psi_{x}, P=-\psi_{y}$꼴이라면 $$\frac{ \partial Q }{ \partial x } -\frac{ \partial P }{ \partial y } =\psi_{xx}+\psi_{yy} = \Delta \psi$$로 작성할 수 있으므로, 해당하는 $\psi$를 찾아보자.

 

즉, $\Delta \psi = R^{4} / (x^{2}+y^{2})^{2}$가 되는 $\psi$를 찾자. 편의 상 $x^{2}+y^{2} = r^{2}$라고 놓고(극좌표 변환) $\psi = \psi(r)$을 구해보자. 이때, $\theta$에 의존하지 않으므로, 라플라시안을 다음과 같이 나타낼 수 있다.

 

$$\Delta \psi(r) = \psi''(r) + \frac{1}{r}\psi'(r) = \frac{R^{4}}{r^{4}}$$

 

이유)

Polar Coordinate의 Laplacian에 대한 유도과정은 아래 링크의 동영상을 참고하자.

https://www.youtube.com/watch?v=wlbZxa66ncg

Polar coordinate의 Laplacian은 $$f_{xx} + f_{yy}   = f_{rr} + \frac{1}{r}f_r + \frac{1}{r^2} f_{\theta\theta}$$가 된다.

 

$f=f(r)$인 중심대칭 함수라면 $\partial f/\partial \theta = 0$이므로 $$\boxed{\Delta f = f''(r) + \frac{1}{r} f'(r)}$$.

이제 양변에 $r$을 곱하면 $$\begin{align} r\psi''(r) + \psi'(r) &= \frac{R^{4}}{r^{3}} \\ \rightarrow \frac{d}{dr}(r\psi'(r)) &= \frac{R^{4}}{r^{3}} \\ \rightarrow r\psi'(r)&=-\frac{R^{4}}{2r^{2}} + C_{1} \\ \rightarrow \psi'(r) &= -\frac{R^{4}}{2r^{3}} + \frac{C_{1}}{r} \\ \rightarrow \psi(r) &=\frac{R^{4}}{4r^{2}} + C_{1}\ln r + C_{2} \end{align}$$

 

이제 $C_{1} = 0, C_{2} = 0$이라고 잡으면, $\psi(r) = \frac{R^{4}}{4r^{2}}$이 된다.

 

따라서, $P, Q$는 다음과 같다.
$$\begin{align} \psi(x, y) &= \frac{R^{4}}{4(x^{2}+y^{2})} \\ P = -\psi_{y} &= \frac{R^{4}\cdot y}{2(x^{2}+y^{2})^{2}} \\ Q = \psi_{x} &= \frac{R^{4}\cdot(-x)}{2(x^{2}+y^{2})^{2}} \end{align}$$

 

따라서, 그린 정리에 의해 $$\begin{align} A' &= \iint_{D} \frac{R^{4}}{(x^{2}+y^{2})^{2}}dA \\ &= \oint_{C} \frac{R^{4}y}{2(x^{2}+y^{2})^{2}}dx - \frac{R^{4}x}{2(x^{2}+y^{2})^{2}}dy \\ &= \boxed{\frac{R^{4}}{2} \oint_{C} \frac{ydx - xdy}{(x^{2}+y^{2})^{2}}} \end{align}$$

 

 

선 적분은 다각형의 경계들의 합으로 표현할 수 있으므로, 이제 $p = (x_{0}, y_{0})$에서 $q = (x_{1}, y_{1})$로 이동하는 하나의 선적분에 대해서 생각해보자.

 

즉, 다음 내용을 생각하자.

 

$$\int_{p \to q} \frac{ydx-xdy}{(x^{2}+y^{2})^{2}}$$

 

이제, $x(t) = x_{0} + (x_{1}-x_{0})t, \quad y(t) = y_{0}+(y_{1}-y_{0})t$라고 두자. 이때, $dx = (x_{1}-x_{0})dt, dy = (y_{1}-y_{0})dt$가 된다.

 

그러면 적분 식은 다음과 같이 쓸 수 있다.
$$\begin{align}&\int_{p \to q} \frac{ydx-xdy}{(x^{2}+y^{2})^{2}}  \\=& \int_{0}^1 \frac{(y_{0}+(y_{1}-y_{0})t)(x_{1}-x_{0})-(x_{0}+(x_{1}-x_{0})t)(y_{1}-y_{0})}{(((x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2})t^{2} + (2x_{0}(x_{1}-x_{0}) + 2y_{0}(y_{1}-y_{0})t + x_{0}^{2}+y_{0}^{2}))^{2}}dt \\=& \int_{0}^{1} \frac{y_{0}x_{1}-y_{0}x_{0}-x_{0}y_{1}+y_{0}x_{0}}{(((x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2})t^{2} + (2x_{0}(x_{1}-x_{0}) + 2y_{0}(y_{1}-y_{0})t + x_{0}^{2}+y_{0}^{2}))^{2}}dt \\=& \int_{0}^{1} \frac{y_{0}x_{1}-x_{0}y_{1}}{(((x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2})t^{2} + (2x_{0}(x_{1}-x_{0}) + 2y_{0}(y_{1}-y_{0})t + x_{0}^{2}+y_{0}^{2}))^{2}}dt \\\rightarrow =& (y_{0}x_{1}-x_{0}y_{1})\int_{0}^{1} \frac{1}{(at^{2}+bt+c)^{2}}dt \\ \end{align}$$

 

이때 $a, b, c$는 각각 $(x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2}, \ 2x_{0}(x_{1}-x_{0}) + 2y_{0}(y_{1}-y_{0}), \ x_{0}^{2}+y_{0}^{2}$이다.

 

울프람 알파로 저걸 적분하면, 그 결과는 다음과 같다.

 

$$\begin{align} \int_{0}^{1} \frac{1}{(at^{2}+bt+c)^{2}}dt &= \frac{2a+b}{\Delta \cdot(a+b+c)}+\frac{4a}{\Delta^{3/2}}\tan ^{-1}\left( \frac{2a+b}{\sqrt{ \Delta }} \right) \\ &- \frac{b}{\Delta \cdot c}-\frac{4a}{\Delta^{3/2}}\tan ^{-1}\left( \frac{b}{\sqrt{ \Delta }} \right) \end{align}$$

 

이때, $\Delta = 4ac - b^{2}$이다.

이제 이걸 모든 변에 대해서 수행해주면 된다!!

 

 

유의할 점은 3개 이상의 점이 하나의 직선 위에 놓여 있을 수 있기 때문에, $\Delta = 0$인 경우 선적분의 값을 $0$으로 리턴하도록 예외처리해줘야 된다.

---

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

// 24558번 Downsizing
// 미적분학
#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>

using namespace std;
#define fastio ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
#define endl '\n'
typedef long double ld;

struct Point {
    ld x, y;
};

ld integrand(ld x0, ld y0, ld x1, ld y1){
    ld a, b, c;
    a = (x1-x0)*(x1-x0) + (y1-y0)*(y1-y0);
    b = 2*x0*(x1-x0) + 2*y0*(y1-y0);
    c = x0*x0 + y0*y0;
    ld d = 4*a*c - b*b;
    const ld eps = 1e-18L;
    if (d < 0 || fabsl(d) < eps) d = 0;
    if (d == 0) return 0.0L;
    
    ld ret;
    ret = (2*a + b)/(d*(a+b+c)) + 4*a/(d*sqrtl(d))*atanl((2*a + b)/sqrtl(d)) - b/(d*c) - 4*a/(d*sqrtl(d))*atanl(b/sqrtl(d));
    ret *= (y0*x1 - x0*y1);
    return ret;
}

int main(void){
    fastio

    ld x, y, r; cin >> x >> y >> r;
    int N; cin >> N;
    vector<Point> points(N);
    for (int i = 0; i < N; ++i){
        cin >> points[i].x >> points[i].y;
    }
    for (int i = 0; i < N; ++i){
        points[i].x -= x;
        points[i].y -= y;
    }
    ld ans = 0;
    for (int i = 0; i < N; ++i){
        ans += integrand(points[i].x, points[i].y, points[(i+1)%N].x, points[(i+1)%N].y);
    }
    ans = ans*r*r*r*r*0.5L;
    
    cout << fixed << setprecision(30);
    cout << ans;

    return 0;
}