$$\newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\N}{\mathbb{N}}\newcommand{\C}{\mathbb{C}} \newcommand{\oiv}[1]{\left] #1 \right[} \newcommand{\civ}[1]{\left[ #1 \right]} \newcommand{\ad}[1]{\text{ad}(#1)} \newcommand{\acc}[1]{\text{acc}(#1)} \newcommand{\Setcond}[2]{ \left\{\, #1 \mid #2 \, \right\}} \newcommand{\Set}[1]{ \left\{ #1 \right\}} \newcommand{\abs}[1]{ \left\lvert #1 \right\rvert}\newcommand{\norm}[1]{ \left\| #1 \right\|}\newcommand{\prt}{\mathcal{P}}\newcommand{\st}{\text{ such that }}\newcommand{\for}{\text{ for }} \newcommand{\cl}[1]{\text{cl}(#1)}\newcommand{\oiv}[1]{\left] #1 \right[}\newcommand{\interior}[1]{\text{int}(#1)}$$

BOJ 14347 / 14346 Radioactive Islands::::Gratus' Blog

BOJ 14347 / 14346 Radioactive Islands

알고리즘 문제풀이/BOJ 2020. 2. 4. 04:14

난이도 : Solved AC 기준 루비 3 (Large), 루비 5 (Small)

출처 : Google CodeJam World Final 2016 

 

 

내가 백준에 제출한 문제중에 가장 어려웠던 것 같다. 미적분학 II 책을 끼고 열심히 띵킹해서 풀었는데, 푸는 과정에서 잠시 백준 채점 서버가 멈추고 (내 코드가 5분 동안 돌아서 그런게 아닐까 싶은 정황이 있고, 완전히 그것때문은 아니라도 상당부분 기여한건 맞는거 같다. :cry:) 여러 일이 있었지만 아무튼 252초의 시간을 쓰면서 맞았다. 우웩;;;

 

먼저 14346 Small을 기준으로 풀이를 작성하고, 14347을 위해 필요한 관찰만 더 소개하는 것으로 하자.

 

1. 문제 설명

$(-10, a)$ 에서 보트가 출발하고, 목적지 $(10, b)$ 로 가려고 한다.

이때, 다음과 같은 조건에 따라 매 초마다 데미지를 입는다.

  • $(0, c)$ 에 있는 섬까지의 거리의 제곱에 반비례하는 데미지, 섬까지의 거리가 $r$ 인 점에서는 $\frac{1}{r^2}$ 의 데미지를 받는다.
  • 어디에 있는, Background radiation damage 를 1씩 받는다.

최소의 데미지를 입으면서 도착할 수 있는 경로를 찾아야 한다.

 

2. 풀이 설명

 

2-1 : 문제 관찰

상황에 대해서 간단하게 몇가지를 먼저 관찰하고 넘어가자.

 

먼저, 섬에 지나치게 가까이 다가가는 것은 엄청나게 큰 손해이다. 섬에 가까울수록 제곱에 비례하는 데미지를 입는 상황이므로, 가급적 섬을 피해 멀리 돌아가고 싶다. 특히 섬 옆을 바로 스쳐 지나가는 상황은 거의 무한대의 데미지를 일시에 받을 수 있으므로 절대 답이 되지 않을 것이다.

 

그러나, 섬으로부터 멀어짐으로써 경로가 지나치게 길어진다면, 그동안 받는 Background radiation 1 에 의한 데미지가 생각보다 커질 수 있다. 그러므로 지나치게 경로가 길어지는 것 역시 피하고 싶다. 필요하다면 섬에 조금 가깝더라도 빠르게 지나가서 배경 방사선 영향을 줄이고 싶을 수도 있다. 

문제에서 제공된 예시. 섬을 피해 최대한 돌면서, 너무 지나치게 멀어지지 않도록 주의한 느낌을 받을 수 있다.

섬이 하나뿐이므로 계산의 편의를 위해 점들을 밀어서, 섬이 원점 $(0, 0)$ 에 있도록 하자. 일반성을 잃지 않음은 자명.

 

2-2 : 적분을 이용한 최적화

주어진 문제의 상황을 어떤 적분식의 최적화로 표현하려고 한다고 생각해 보자. 먼저, 경로가 $y = f(x)$ 의 함수를 따라 가는 경로라고 생각하고, 이 경로를 가는 동안 받는 방사선을 생각해 보기로 하자.

 

배경으로부터 받는 데미지는 시간당 1의 속도로 움직이고 있으므로 전체 경로의 길이만큼이다. 경로 곡선의 길이는 다음과 같은 식으로 구할 수 있음이 잘 알려져 있다.

$$I_1 = \int_{-10}^{10} \sqrt{1 + f'(x)^2} dx$$

여기에 섬으로부터 받는 양을 계산하려면, 각 구간 dx마다 받는 방사선에 그 구간을 지나는 시간을 고려하여 적분한다.

$$I_2 = \int_{-10}^{10} \frac{1}{x^2 + f(x)^2} {\sqrt{1 + f'(x)^2}} dx$$

이 두 함수의 합을 최적화하고 싶다.

 

다양한 접근이 있지만, 나는 변분법을 이용하여 풀었다.

변분법 (Calculus of Variations) 은 어떤 함수들의 공간에서 뭔가를 최적화하는 함수 $f$ 를 찾는 문제를 푸는 방법이다.  대표적인 예시는, 어떤 점에서 다른 점으로 '중력의 가속력으로만 강하하는' 최단 시간 경로가 사이클로이드 형태라는 잘 알려진 결과를 들 수 있다. 참고 : The Brachistochrone

구체적으로, 이 문제에서 쓰는 것은 오일러-라그랑주 방정식이다. 다음과 같은 적분의 값을 최적화하고 싶다고 하자.

$$I[f] = \int_{x_0}^{x_1} \mathcal{L}(x, f(x), f'(x) \cdots, f^{k}) \text{ d}x$$

단, 여기서 $f^{k}$ 는 $f$를 $x$에 대해 $k$번 미분한 함수를 의미한다.

오일러-라그랑주 방정식은 이 식의 최적해(Functional을 최소화하는 함수 $f$) 가 다음 조건을 만족한다는 것이다.

$$\sum_{i = 0}^{k} \frac{\text{d}^i}{\text{d} x^i}\frac{\partial \mathcal{L}}{\partial f^{i}} = 0$$

 

이 문제에서는 이계도함수 이상이 필요하지 않으므로, $f$와 $f'$에 대해서만 한정해서 쓰면, 

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

 

이제, 이 식을 우리가 가지고 있는 앞서의 함수에 적용하면 다음과 같은 식을 얻는다. 

$$\frac{f''}{(1+f'^2)^{3/2}}\left(\frac{1}{x^2+f^2} +1 \right) - \frac{2x}{(x^2+f^2)^2}\frac{f'}{(1+f'^2)^{1/2}} =- \frac{2f}{(x^2+f^2)^2}\frac{1}{(1+f'^2)^{1/2}}$$

 

조금 정리하면, $f''$ 을 $f'$ 과 $f$, $x$ 로 풀어 쓸 수 있다. 

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

 

 

2-3 : 초기값 문제 해결

결국, 식정리와 미적분 매드무비가 알려주는 것은 어떤 $x, f(x), f'(x)$ 의 초기값이 구해지면 그것으로부터 $f''(x)$ 의 값을 구할 수 있다는 것이고, 이를 다시 $f'$ 에 먹이고 $f$에 먹여서 다음 값을 구하는 식으로 구간에 따른 수치적분을 수행할 수 있다는 사실이다. 초기 $x$와 $f(x)$ 의 값은 주어져 있으므로, $f'$의 초기값만 있으면 (최초 방향) 이를 이용하여 전체 경로를 그려볼 수 있다.

 

그러나, $f'$ 의 초기값에 대해서는 아무런 정보가 주어져 있지 않다. 추가로, 우리는 도착점이 $(10, b)$ 여야 한다는 강제적인 조건이 있는데 이 부분은 전혀 반영되어 있지 않다.

 

이를 해결하기 위해, 먼저 $f'$ 의 초기값이 증가한다면, 즉 다른 모든 조건이 그대로 유지되는 와중에 최초 방향만 위로 맞추고 출발한다면 적어도 아래로 맞춰져 있을 때보다는 높은 위치에 도착할 것이라는 추측을 해볼 수 있다. (수식에 대입해서 이를 계산해도 어떻게 잘 나올 것 같지만, 굳이 이를 전부 해보진 않았다) 

 

그렇다면, 역으로 일단 적당히 아무렇게나 내보내 본 다음, $(10, h)$ 에 도착한다면 이 도착점을 보고 $f'$의 초기값을 조정하는 방법으로 이분 탐색을 시도할 수 있을 것이다. 특히, 대략적으로 '섬 위를 지나는 경로' 와 '섬 아래를 지나는 경로' 를 하나씩 잡아 놓고, 위쪽 경로 중 최적이면서 원하는 위치에 도착하는 경로와, 아래쪽 경로 중 최적이면서 원하는 위치에 도착하는 경로를 찾아 비교할 수 있다. 이 두 경로 중 더 낮은 dose 를 받는 경로가 최적일 것이다.

 

 

2-4. 두 개의 섬

Large 문제는 섬 두 개가 있는 상황에서 이 문제를 똑같이 풀어야 한다. 크게 달라지는 점은 없고, 2-2에서 유도한 적분식에서 더이상 $(0, 0)$ 에만 섬이 있는 편리한 상황이 아니라, 어딘가에 섬이 또 있다는 사실 때문에 여기저기에 $f$ 대신 $f-c_0$ 나 $f-c_1$ 로 치환되어야 한다.

 

또한, 초기값 문제를 해결하는 과정에서 경로를 둘이 아닌 셋으로 나누어야 한다. 위쪽에 있는 섬보다 더 위로 돌아 가는 Upper path, 아래 있는 섬보다 더 아래로 돌아 가는 Lower path, 두 섬 사이를 과감하게 통과하는 Middle path 를 생각해 볼 수 있는데, Upper path 나 Lower path 가 답이라고 생각하기 쉽지만 생각보다 많은 경우에, 섬들이 서로 적당히 멀리 떨어져 있다면 가운데로 파고들어서 나가는 편이 낫고 섬을 피하려고 돌다 보면 굉장히 오래 걸려서 그동안 받는 Background가 더 커진다는 것을 알 수 있다. 

 

 

사람마다 숫자를 보고 느끼는 것은 많이 다르겠지만, 나는 처음에 코딩을 실제로 해보고 Small 을 풀면서 랜덤 데이터를 먹여보고 느꼈던 점은 생각보다는 섬에 가까이 붙는다는 것이다. 사실, 섬에서 2km 떨어진 곳까지 가도 0.25의 추가 데미지를 받는 셈인데, 0.25를 받으면서 3의 거리를 달리는 것이 경로가 0.75만큼 길어지고 섬에서 멀어지는 것보다 이득이다. 

 

3. 코드

코드는 상당히 더러운데, 우선 디버깅용으로 추가된 부분도 꽤 많고, 섬이 하나일 때와 두 개일 때를 나누어서 풀었다. 백준에서 제출했을 때는 제한시간이 5분임에도 불구하고 시간초과를 받았는데, 어차피 최종 결과에서 0.001이라는 넉넉한 실수오차를 허용하기 때문에 구간을 크게 잡고 대-충 적분해도 된다. 적분 스텝의 크기 0.00001 (-10부터 10까지 200만 스텝) 은 4분 10초로 간신히 들어오고, 0.000005 (400만 스텝) 은 시간 초과를 받았다. 

 

더보기
#include <bits/stdc++.h>
#pragma GCC optimize("O3")
#pragma GCC optimize("Ofast")
#pragma GCC target("avx,avx2,fma")
#define ll long long
#define eps 1e-7
#define all(x) ((x).begin()),((x).end())
#define usecppio ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
using namespace std;
using pii = pair<int, int>;
#define int ll
double twosolve();
inline double val_ypp(double &y, double &yp, double &x)
{
    double u = 2*x*yp-2*y;
    double uu = (x*x+y*y);
    uu *= (uu+1);
    double lu = (1+yp*yp);
    return u*lu/uu;
}
inline double recieve_dose(double &xpos, double &ypos, double &interlen)
{
    double dd = 1 + 1/(xpos*xpos+ypos*ypos);
    return dd*interlen;
}
pair <double, double> trial(double init_yp, double start_y)
{
    double yp = init_yp;
    double y = start_y;
    double x = -10;
    double xp = 0.00001;
    double dose = 0;
    while(x<10)
    {
        double ilen = sqrt(1+yp*yp)*xp;
        dose += recieve_dose(x,y,ilen);
        x += xp;
        double ypp = val_ypp(y, yp, x);
        y += yp*xp;
        yp += ypp*xp;
    }
    return {y, dose};
}

double solve()
{
    double a, b, c;
    int n;
    scanf("%lld",&n);
    if (n == 2)
        return twosolve();
    scanf("%lf %lf %lf",&a,&b,&c);
    a -= c;
    b -= c;
    double up_ans, down_ans;

    double lo = -100;
    double hi = -a/10.0;
    while(abs(lo-hi)>eps)
    {
        double mid = (lo+hi)/2;
        auto guess = trial(mid,a);
        double endpos = guess.first;
        if (endpos<b)
            lo = mid;
        else hi = mid;
    }
    down_ans = trial(lo,a).second;

    double lo2 = -a/10;
    double hi2 = 100;
    while(abs(lo2-hi2)>eps)
    {
        double mid = (lo2+hi2)/2;
        auto guess = trial(mid,a);
        double endpos = guess.first;
        if (endpos<b)
            lo2 = mid;
        else hi2 = mid;
    }
    up_ans = trial(lo2,a).second;

    return min(up_ans, down_ans);
}

int32_t main()
{
    int TC;
    cin >> TC;
    for (int tt = 1; tt<=TC; tt++)
    {
        printf("Case #%lld: ",tt);
        printf("%.10lf\n",solve());
    }
}

inline double val_ypp_2(double &y, double &yp, double &x, double &c0, double &c1)
{
    double lu = (1+yp*yp);
    double tmp_1 = (x*x+(y-c0)*(y-c0));
    double tmp_2 = (x*x+(y-c1)*(y-c1));
    double up1 = -2*(y-c0)/(tmp_1*tmp_1);
    double up2 = -2*(y-c1)/(tmp_2*tmp_2);
    double down1 = 1/tmp_1;
    double down2 = 1/tmp_2;

    return lu*(up1+up2)/(down1+down2+1);
}
inline double recieve_dose_2(double &xpos, double &ypos, double &interlen, double &c0, double &c1)
{
    double t = abs(ypos-c0), t2 = abs(ypos-c1);
    return (1/(xpos*xpos+t*t)+1/(xpos*xpos+t2*t2)+1)*interlen;
}
pair <double, double> trial_2(double init_yp, double start_y, double &c0, double &c1)
{
    double yp = init_yp;
    double y = start_y;
    double x = -10;
    double xp = 0.00001;
    double dose = 0;
    int cnt = 0;
    while(x<10)
    {
        cnt++;
        double ilen = sqrt(1+yp*yp)*xp;
        dose += recieve_dose_2(x,y,ilen,c0,c1);
        x += xp;
        double ypp = val_ypp_2(y, yp, x,c0,c1);
        y += yp*xp;
        yp += ypp*xp;
    }
    return {y, dose};
}


double twosolve()
{
    double a, b, c0, c1;
    scanf("%lf %lf %lf %lf",&a,&b,&c0,&c1);
    if (c0 > c1) swap(c0, c1);
    double up_ans, mid_ans, down_ans;

    double lo = -100;
    double hi = -(c0-a)/10.0;

    while(abs(lo-hi)>eps)
    {
        double mid = (lo+hi)/2;
        auto guess = trial_2(mid,a,c0,c1);
        double endpos = guess.first;
        if (abs(endpos-b)<=eps)
            break;
        if (endpos<b)
            lo = mid;
        else hi = mid;
    }
    double mid1 = (lo+hi)/2;
    down_ans = trial_2(mid1,a,c0,c1).second;

    double lo2 = -(c1-a)/10.0;
    double hi2 = 100;
    while(abs(lo2-hi2)>eps)
    {
        double mid = (lo2+hi2)/2;
        auto guess = trial_2(mid,a,c0,c1);
        double endpos = guess.first;
        if (abs(endpos-b)<=eps)
            break;
        if (endpos<b)
            lo2 = mid;
        else hi2 = mid;
    }
    double mid2 = (lo2+hi2)/2;
    up_ans = trial_2(mid2,a,c0,c1).second;
    double lo3 = -(c0-a)/10.0;
    double hi3 = -(c1-a)/10.0;
    while(abs(lo3-hi3)>eps)
    {
        double mid = (lo3+hi3)/2;
        auto guess = trial_2(mid,a,c0,c1);
        double endpos = guess.first;
        if (abs(endpos-b)<=eps)
            break;
        if (endpos<b)
            lo3 = mid;
        else hi3 = mid;
    }
    double mid3 = (lo3+hi3)/2;
    mid_ans = trial_2(mid3,a,c0,c1).second;
    return min(min(mid_ans,up_ans), down_ans);
}

'알고리즘 문제풀이 > BOJ' 카테고리의 다른 글

BOJ 15310 아티스트  (2) 2020.03.07
BOJ 13558 등차수열의 개수  (3) 2020.02.29
BOJ 16709 Congruence Equation  (2) 2020.02.01
BOJ 12728 n제곱 계산 / 12925 Numbers  (2) 2020.01.30
BOJ 15745 Snow Boots  (0) 2020.01.23
admin