BOJ 14346 - Radioactive Islands (Small) (Python3)

여러모로 의미가 있는 문제이다. 벼르고 있긴 했지만 그래도 루비라서 손 댈 생각이 없었는데, 정확하게 루비 한 문제로 내 solved.ac 티어를 다이아까지 올려 놓을 수 있는 상황이 되어 동기가 생겼다. 저녁부터 씨름하다 딴짓하다 씨름하다 딴짓하다를 반복한 끝에 오늘 새벽 3시 20분이 되어서야 마무리지을 수 있었다. 영 다른 함수를 적분하는 대참사가 벌어져서...

역시 영어로 된 문제다. 이제는 역사의 뒤안길로 사라진 Google Code Jam의 2016년 World Finals 문제란다. 명성에 걸맞게 끔찍하게 어려운 난이도를 자랑한다. 나는 PS하면서 변분법을 쓰게 될 줄은 몰랐다...

원문 뒤에 번역을 붙이는 방식으로 문제 소개를 하려고 했는데 풀이과정도 상당히 길어질 것 같고 그냥 100% 내 번역으로 집어넣으려고 한다. O-ey 번역을 날먹하는 건 영 아니다.

 

그럼 본격적으로 풀이에 돌입해 보자.

---

문제

당신은 \((-10, A)\)에 위치한 보트를 \((10, B)\)까지 몰고 가려고 한다. 좌표의 단위는 km 단위이고, 당신의 보트는 일정한 속력 1km/h로 이동한다. 당신은 당신 마음대로 보트를 움직일 수 있다. 보트를 좌표평면상의 한 점이라고 가정하자.

이 구역에는 \(N\)개의 섬이 있다. 섬들 또한 좌표평면상의 한 점이라고 가정하자. \(i\)번째 섬의 좌표는 \((0,C_i)\)이다.

이 구역은 방사능을 띄고 있으며, 당신은 구역 내 어느 위치에 있건 간에 한 시간당 1mSv씩 일정하게 방사능에 피폭된다. 또한 구역 내 섬들도 방사능을 띄고 있으며, 한 시간에 \(i\)번째 섬과 떨어진 거리 \(D_i\)km의 제곱에 반비례하는 양인 \(D_i^{-2}\)mSv씩 일정하게 방사능에 피폭된다.

정확하게 말하자면, \(i\)번째 섬과 당신 간의 거리가 함수 \(D_i(t)\)로 나타내어지고 \(0\)초부터 \(X\)초까지 여행한다고 하자. 그러면 당신이 \(i\)번째 섬의 방사능에 피폭된 양은 \(\int_0^XD_i(t)^{-2}dt\)mSv이다.

당신은 섬과 정확하게 같은 위치에 있지 않는 한, 섬에 원하는 만큼 가까이 접근할 수 있다.

당신이 피폭량을 최소화하는 경로를 따라서 보트를 운행할 시 그 때 받게 되는 방사능의 총량을 구하라.

 

입력

첫 번째 줄에 테스트케이스의 수 \(T\)이 주어진다. \((1\le T\le 20)\)

각 테스트케이스는 두 줄로 구성된다. 첫 번째 줄에 정수 \(N\)과 두 실수 \(A, B\)가 공백으로 구분되어 주어진다. \(N\)은 구역 내 섬의 수이고. \(A,B\)는 각각 출발점과 도착점의 \(y\)좌표이다. \((N=1;\) \(-10.00\le A,B\le 10.00)\)

두 번째 줄에는 \(N\)개의 실수 \(C_i\)가 공백으로 구분되어 주어진다. \(i\)번째 수 \(C_i\)는 \(i\)번째 섬의 \(y\)좌표이다. \((-10.00\le C_i\le 10.00)\)

입력으로 주어지는 실수는 정확하게 소수점 아래 두 자리까지만 주어진다.

 

예제 입력)

1
1 1.00 -2.00
0.00

 

출력

각 테스트케이스마다 Case #x: y의 행식으로 출력하라. x는 1부터 시작하는 테스트케이스의 번호이고, y는 mSv 단위로 주어진 최소 방사성 피폭량이다.

y와 절대 혹은 상대 오차가 \(10^{-3}\) 이하인 수를 출력하더라도 정답으로 인정된다.

 

예제 출력)

Case #1: 21.81

 

---

내 코드

tol = 0.001
step = 20000
h = 20/step
def f(x,y,u): return ((2*x*u-2*y)*(1+u**2))/((x**2+y**2)*(x**2+y**2+1))
def RK4(y0,u0):
    x0 = -10
    liy,liu = [y0],[u0]
    for _ in range(step):
        m1 = u0; k1 = f(x0,y0,u0)
        m2 = u0 + h*k1/2; k2 = f(x0+h/2,y0+h*m1/2,u0+h*k1/2)
        m3 = u0 + h*k2/2; k3 = f(x0+h/2,y0+h*m2/2,u0+h*k2/2)
        m4 = u0+h*k3; k4 = f(x0+h,y0+h*m3,u0+h*k3)
        x0 += h
        y0 += h*(m1+2*m2+2*m3+m4)/6
        u0 += h*(k1+2*k2+2*k3+k4)/6
        liy.append(y0); liu.append(u0)
    return (y0,u0,liy,liu)
def length(liy,liu):
    def F(x,y,u): return pow(1+u**2,0.5)*(1+1/(x**2+y**2))
    ans = F(-10,liy[0],liu[0]) + F(10,liy[-1],liu[-1])
    for i in range(1,step):
        if i%2: ans += 4*F(-10+h*i,liy[i],liu[i])
        else: ans += 2*F(-10+h*i,liy[i],liu[i])
    return ans*(h/3)
for _ in range(int(input())):
    N,A,B = map(float,input().split())
    C = float(input())
    slope = (C-A)/10
    A -= C; B -= C
    ans = []
    l,r = -2,slope
    while r-l > tol:
        mid = (l+r)/2
        y,u,liy,liu = RK4(A,mid)
        if y > B: r = mid
        else: l = mid
    ans.append(length(liy,liu))
    l,r = slope,2
    while r-l > tol:
        mid = (l+r)/2
        y,u,liy,liu = RK4(A,mid)
        if y > B: r = mid
        else: l = mid
    ans.append(length(liy,liu))
    print(f'Case #{_+1}: {min(ans)}')

 

실수 오차가 널널하지만 않았더라도 난이도는 저세상으로 뛰었을 문제다. 물론 지금도 결코 쉬운 문제는 아니다. 농담으로라도 쉽다고 해서는 안 되는 문제.

 

계산의 편의성을 위해서, 영역 내 좌표계를 유일한 섬이 원점에 있도록 구성하자. 그러면 여정은 \((-10,A-C)\)에서 시작해서 \((10,B-C)\)에서 끝난다. 물론 이것도 계산의 편의성을 위해서 \(A,B\)를 \(A-C, B-C\) 대신 사용하자. 그냥 처음부터 원점에 섬이 위치했다고 보자는 이야기다.

상식적으로, 여정 내 같은 \(x\)좌표는 단 한 번만 지날 것이다. 오른쪽으로 가다가 왼쪽으로 되돌아가고, 다시 오른쪽으로 나아가는 상황은 발생하지 않을 것이라는 말이다. 왜? 시간이 낭비되고, 그만큼 방사선 피폭량은 늘어날 것이니까. 그러므로, 경로는 \(y=f(x)\)의 형태를 띄고 있다고 가정해도 무방하다.

그리고 하나 더. 컴퓨터 프로그래밍이니만큼 경로는 매끄럼지 않게 그려지겠으나 실제로 배가 그리는 궤적은 매끄러울 것이다. 이게 굉장히 중요하다. 매끄럽다는 것은 미분 가능이면서 연속을 의미하는 것이니까. \(y=|x|\)처럼 첨점이 존재하는 경우는 무시해도 된다는 것이다. 배가 텔레포트하는 일은 없을 테니 당연히 불연속함수일 수도 없을 것이고.

 

그러면 배경 방사선 양은 경로 곡선의 길이와 같은 \(I_1=\int_{-10}^{10}\sqrt{1+f'(x)^2}dx\)일 것이고, 섬으로부터 받는 방사선의 양을 고려한다면 미소구간 \(dx\)마다 받는 방사선, 즉 거리 제곱의 역수 \(\frac1{x^2+f(x)^2}\)에 해당 구간을 지나는 미소 시간을 곱해 준 다음 적분하면 된다. 그런데 배의 속력은 1km/h이므로, 미소 시간 대신 미소 거리를 써도 무방하다. 미소 거리는 \(ds=\sqrt{1+f'(x)^2}\)과 같다. 따라서 섬에서 받는 방사선의 양은 \(I_2=\int_{-10}^{10}\frac1{x^2+f(x)^2}\sqrt{1+f'(x)^2}dx\)이 된다.

즉, 우리가 최소화하고 싶은 것은 다음 범함수의 값이다. (이제부터 편의상 \(f(x),f'(x)\) 대신 \(f,f'\)을 쓰겠다.)

$$J=\int_{-10}^{10}\left(1+\frac1{x^2+f^2}\right)\sqrt{1+f'^2}dx$$

 

이제 변분법을 배워야 한다. 유도과정은 조금 길고 복잡하고 (나야 좋지만 이 글을 보는 사람은 별로 좋아하지 않을) 수학 투성이이기 때문에 중간과정을 살짝 생략시켜 주자면, 범함수 \(J=\int_i^o\mathcal L(f,f';x)dx\)이 극값을 가지게 하는 함수 \(\mathcal L\)은 다음과 같은 미분방정식을 만족시킨다. (이를 오일러-라그랑주 방정식이라고 한다.)

$$\frac{\partial\mathcal L}{\partial f}=\frac{d}{dx}\frac{\partial\mathcal L}{\partial f'}$$

이 문제에서 \(\mathcal L(f,f';x)=\left(1+\frac{1}{x^2+f^2}\right)\sqrt{1+f'^2}\)이다. 좌변은 \(\frac{\partial\mathcal L}{\partial f}=-\frac{2f}{x^2+f^2}\sqrt{1+f'^2}\)이고, 우변은 \(\frac{d}{dx}\frac{\partial\mathcal L}{\partial f'}=\frac{d}{dx}\left(1+\frac{1}{x^2+f^2}\right)\frac{f'}{\sqrt{1+f'^2}}=\left(1+\frac1{x^2+f^2}\right)\frac{f''(x)}{(1+f'(x))^{\frac{3}{2}}}-\frac{2x}{(x^2+f^2)^2}\frac{f'}{\sqrt{1+f'^2}}\)이 된다. 이 미분의 향연 뒤로 식을 예쁘게 정리해 보면, 다음과 같이 \(f''=F(x,f,f')\)의 형태로 나타낼 수 있다.

$$f''=\frac{2xf'-2f}{x^2+f^2}\frac{1+f'^2}{1+x^2+f^2}$$

 

이 2계미분방정식은 다음과 같은 연립 1계 미분방정식으로 치환할 수 있다.

$$\left\{\begin{matrix}u=f'\\ u'=F(x,f,u)=\frac{2xu-2f}{x^2+f^2}\frac{1+u^2}{1+x^2+f^2} \end{matrix}\right.$$

미분방정식의 수치해법은 Euler's method, improved Euler's method, Runge-Kutta method 등 여러 가지가 있지만, 아무래도 Euler's method는 전체절단오차가 \(\mathcal O(h)\), improved Euler's method라고 할지라도 전체절단오차가 \(\mathcal O(h^2)\)이다. 아무리 실수오차가 널널하다 한들 그대로 적용하기에는 조금 무리가 있다. 4차 Runge-Kutta method는 전체절단오차가 \(\mathcal O(h^4)\)밖에 되지 않으며, 이는 더 적은 횟수의 연산만으로도 실수오차를 충분히 줄일 수 있음을 의미한다.

같은 시간복잡도를 가지면서도 메모리 사용량이 1/4로 줄어든 4차 Adams-Bashforth-Moulton method를 사용할 수도 있지만, 처음 세 값을 RK4로 구해야 한다. 메모리도 0.5GB로 널널하니 그냥 RK4를 쓰자.

 

RK4를 사용해서 이 미분방정식을 풀려고 하는데 맙소사, \(f'(-10)\)이 주어져 있지 않다. 초기값 문제를 풀려고 하는데 초기값 대신 경계값 조건 \(f(10)=B\)이 생겨 버렸고, 슬프게도 이건 RK4를 쓰는 데 그렇게 큰 도움이 되지 못한다.

그렇다면 임의로 초기값 \(f'(-10)=u(-10)\)을 지정하고 RK4를 열심히 돌린 다음 결과값으로 나온 \(f(10)\)이 \(B\)보다 크면 \(u(-10)\)을 조금 줄이고, 반대로 \(f(10)\)이 \(B\)보다 작으면 \(u(-10)\)을 조금 키우면 된다. 이 방법에 가장 알맞은 방법은? 이분탐색이다.

 

이 문제는 최적홰를 두 개 가진다. 섬 위로 돌아가느냐, 섬 밑으로 돌아가느냐. 그리고 \(A, B, C\)의 값 때문에 \(u(-10)\)이 그렇게 크가나 작은 값이 아닐 테고. 그러면 섬 위로 돌아가는 경로는 \(l=\frac{-A}{10}, r=2\)로 두고, 섬 아래로 돌아가는 경로는 \(l=-2,r=\frac{-A}{10}\)로 두고 이분탐색을 돌리면 된다.

굉장히 작위적으로 보이는 이 \(\frac{-A}{10}\)이라는 값은 시작점에서 섬을 향해 배를 일직선으로 모는 케이스를 의미한다. 위에서 오일러-라그랑주 방정식을 소개할 때 범함수 \(J\)가 최솟값을 가지는 게 아니라 극값을 가진다고 이야기했었다. 그렇다. 섬으로 그대로 직진하는 경로는 방사능 피폭량을 최대한으로 만들어버리기 때문에 극값이다. 이 경우를 원천봉쇄하기 위해서 일부러 이분탐색의 시작점을 이 값으로 두었다.

 

이제 정말로 열심히 RK4만 돌리면 된다. 위의 연립미분방정식을 푸는 RK4 iteration은, 하나의 step을 \(h\)라고 둘 시 다음과 같다. 물론 초기값 \(x_0=-10, f_0=f(-10)=A, u_0=mid\)이다.

$$x_{n+1}=x_n+h,\,\,f_{n+1}=f_n+\frac{h}{6}(m_1+2m_2+2m_3+m_4),\,\,u_{n+1}=u_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$$

$$m_1=u_n,\,\,m_2=u_n+\frac12hk_1,\,\,m_3=u_n+\frac12kh_2,\,\,m_4=u_n+kh_3$$

$$k_1=F\left(x_n,f_n,u_n\right),\,\,k_2=F\left(x_n+\frac12h,f_n+\frac12hm_1,u_n+\frac12hk_1\right)$$

$$k_3=F\left(x_n+\frac12h,f_n+\frac12hm_2,u_n+\frac12hk_2\right),\,\,k_4=F\left(x_n+h,f_n+hm_3,u_n+hk_3\right)$$

 

마지막으로 RK4에서 나오는 각 \(f_n,u_n\)값은 버리지 말고 차곡차곡 모아두었다가 방사능 피폭량 \(J=\int_{-10}^{10}\left(1+\frac1{x^2+f^2}\right)\sqrt{1+f'^2}dx\)을 심프슨 공식으로 계산하는 데 사용하자. 심프슨 공식에 익숙해져 있지 않다면 9206번 문제로 연슴을 해 보는 것을 권한다. 조만간 해당 문제 풀이도 여기에 올리겠다.

여하튼, 어떻게 잘 버무려 보면 그닥 어렵지 않게 AC를 받을 수 있다. 상대 오차가 0.1%씩이나 나도 통과시켜주는 관대함 덕분에 이분탐색의 허용오차를 0.001이라는 터무니없이 큰 값으로 두어도 통과할 수 있었다. 시간은 걱정 안 해도 된다. 30초나 준다.

 

강화판 문제로 섬이 두 개 있는 14347번 문제가 있는데, 그 문제를 풀려면 평행이동으로 섬들 위치를 잘 지정해 줘야 하고, 섬들 위로 가는 경로와 아래로 가는 경로 말고도 두 섬 사이를 지나가는 경로를 감안해 줘야 하며, 변분법에 사용해 줄 범함수가 많이 복잡해진다. 대신 제한시간이 300초로 늘어난다. 파이썬으로 풀면 15분이다.

 

이로서 14346번의 풀이를 마친다.