[백준 9206] 나무 말고 꽃 (Platinum II)

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

9206번: 나무 말고 꽃

수치해석 문제들이 다 그렇듯이 수많은 시도를 해야 풀 수 있다... 일단 constexpr로 각종 상수를 정의해 주자. constexpr value들은 컴파일 시점에 평가되기 때문에, #define 같은 매크로의 기능을 거의 모두 쓸 수 있으면서 그것보다 더 안전하다고 할 수 있겠다.

constexpr long double PI = 3.1415926535897932384626l;
constexpr int STEP = 20000;

PI의 맨 끝에 l을 안 붙이면 long double이 아니라 double로 인식이 돼서 WA를 받는다. (난 그것 때문에 개고생했다) 정확히는, PI 자체는 long double형 변수가 되는데 PI를 init하는 값은 double이라 정확도를 자기 스스로 깎아먹는 짓을 하는 것이다.

주어진 식은 부정적분을 초등함수 내에서 나타낼 수 없으므로 수치해석 기법을 써야 한다. 떠오르는 접근 방법은 테일러 급수와 심슨 공식의 2가지가 있을 수 있는데 난 첫 번째로는 WA를 받고 두 번째로 AC를 받았다. 그 이유를 살펴보자.

테일러 급수

주어진 수식

$$ \pi {\left( ae^{-x^2} + b\sqrt{x} \right)^2} $$

는 곱셈공식에 의해 기초 함수들의 덧셈으로 쪼개질 수 있다. 이것들을 일일이 테일러 전개해 부정적분을 무한급수 꼴로 나타낸 뒤, h를 대입하면 끝. 그러나 다음과 같이 주어지는 테일러 급수의 오차항을 생각해 보자.

$$ \left| \triangle V \right| \leq {{\sup{\left| f^{(n+1)}(t) \right|} \times \left| x-a \right|^{n + 1}} \over {(n + 1)!}} $$

(n은 n차항까지 사용했음을 의미한다.) 이 식은 오차의 대략적인 모습을 제공하는데, 주요하게 보아야 할 항은 \( \left| x - a \right|^n \). 이 값이 최대 10까지 들어올 수 있다. 따라서 n을 작게 잡으면 10^n / n!의 값은 2000을 넘어서는 끔찍한 오차를 보일 수도 있다... 물론 테일러 급수는 결국 n의 값을 무한대로 늘리면 오차가 0으로 수렴한다는 것이 보장되므로 n을 겁나 크게 잡으면 되는 거 아닌가 하고 생각할 수가 있는데, 안타깝게도 long double의 계산오차 때문에 그렇게 할 수 없다. 이공계 대학생이라면 소수들의 덧셈에서 유효숫자는 가장 작은 것을 따라간다는 사실을 알고 있을 것이다. 비슷한 논리가 long double의 계산오차에 그대로 적용된다.

심슨 공식

심슨 공식은 라그랑주 보간법을 이용해 임의의 곡선을 2차함수 또는 3차함수의 연속으로 근사하면서 적분한 결과물이다. 유도 과정은 생략하고, 우리는 3차함수의 근사를 사용하도록 하자.

$$ \int_0^{3h}{f(x)dx} \sim {1 \over 8} \left( f(0) + 3f(h) + 3f(2h) + f(3h) \right) $$

물론 이런 하찮은 근사로 10^-4 수준 정확도가 나올 리가 없다. 그래서 마치 구분구적법 하는 것처럼 구간을 잘게 쪼개서 각 구간마다 심슨 공식을 따로 사용할 것이다.

long double simpson_a_part(const Bottle &b, long double s, long double e) {
    long double m1 = (2 * s + e) / 3, m2 = (s + 2 * e) / 3;
    return (e - s) / 8 * (b(s) + 3 * b(m1) + 3 * b(m2) + b(e));
}

long double simpson_integration(const Bottle &b, long double s, long double e) {
    long double step = (e - s) / STEP, ret = 0;
    for (int i = 0; i < STEP; ++i) {
        ret += simpson_a_part(b, s + i * step, s + (i + 1) * step);
    }
    return ret;
}

 

클래스 Bottle은 함수 객체다. 첫 번째 함수에서 볼 수 있듯이 Bottle의 인스턴스를 마치 함수처럼 취급할 수 있다. Bottle 클래스도 이해를 돕기 위해 삽입해야겠다.

class Bottle {
public:
    long double a, b;
    long double operator()(long double x) const {
        long double R = a * std::exp(-x * x) + b * std::sqrt(x);
        return R * R * PI;
    }
};

 

아니면 람다 함수를 사용해서 함수 객체를 만드는 것도 가능하다.

using Bottle = std::function<long double(long double)>;

Bottle bottle_factory(long double a, long double b) {
    return [a, b](long double x) -> long double {
        long double R = a * std::exp(-x * x) + b * std::sqrt(x);
        return R * R * PI;
    };
}

 

여기서 std::exp 대신 그냥 <cmath>에 정의된 exp를 쓰면 long double이 강제로 double로 바뀌기 때문에 오차가 난다. sqrt도 마찬가지.

오차항의 상한은 다음과 같다.

$$ \left| {(b - a)^5 \over {6480}} \min_{a\leq x\leq b} f^{(4)}(x) \right| $$

테일러 공식 때와 달리, (b - a) 값은 우리가 조절 가능한 값이다. 그래서 오차항을 관리하는 데 그렇게 큰 신경을 쓸 필요가 없다.