Category: Mathematics

실수의 실수: 스페셜 저지만 붙인다고 끝나는 것이 아니다

$% TeX defines$ $\newcommand{\eps}{\varepsilon_{\mathrm{m}}}$ $\newcommand{\F}{\mathbb{F}}$

프로그래밍 문제를 출제하다 보면 가끔 실수를 다루는 문제를 출제하고 싶어질 때가 있습니다. 웬만하면 모든 걸 정수로 하면 좋겠지만 문제 상황이 그렇지 못할 경우도 존재하죠. 그러나 우리가 다루는 실수는 사실 정확히 말하면 실수가 아니기 때문에 조심해야 할 점들이 꽤나 있습니다.

이 글은 floating point 실수 연산의 오류 분석과, 이를 기반으로 프로그래밍 문제에서 좀 더 옳은 정답 데이터를 만드는 방법에 대해 소개합니다. 굳이 문제나 데이터를 만들지 않더라도 일반적으로 프로그래밍 문제에서 실수를 다루는 데에 참고하면 도움이 될 것입니다.

IEEE 754

By Vectorization: Stannered – Own work based on: Float example.PNG, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=3357169

아마 float이나 double 등을 처음 배울 때 한 번쯤은 봤을 법한 그림입니다. 지금은 잊어버렸을 수도 있겠습니다만, 아무튼 컴퓨터가 실수를 그렇게 정확하게 나타내지는 못한다는 것쯤은 알고 계시리라 생각합니다.

사실 float과 double은 아무 실수나 나타낼 수 있는 것이 아니며, 분모가 $2^n$인 유리수만을 나타낼 수 있습니다. 심지어 $0.1$과 같은 우리가 보기에 간단해 보이는 유리수도 분모를 $2^n$으로 나타낼 수 없기에, double x = 0.1을 하면 $0.1$이 아니라 $0.1$을 최대한 근사한 값¹이 들어가게 됩니다.

양수 float 값 $q$는 일반적으로 다음과 같이 저장됩니다²:

\[ q = 2^{\left(A-\mathit{bias}\right)} \times \left(1+B\right). \]

여기서 $A$는 음이 아닌 정수이고, $B$는 소숫점 아래에 쓸 수 있는 $0 \le B < 1$의 이진수 소수입니다. 수의 크기에 따라 $A$를 조절해서 $\left[1, 2\right)$보다 크거나 작은 수를 표현할 수 있습니다.

int, long long 등에 제한이 있듯이 float, double도 한정된 비트를 사용하는 만큼 제한이 있습니다. 아래 정보는 일반적인 온라인 저지 환경에서의 GCC를 가정합니다.

	$A$	$\mathit{bias}$	$B$
`float` ($32$비트)	$8$비트	$127$	$23$비트
`double` ($64$비트)	$11$비트	$1\,023$	$52$비트
`long double` ($80$비트)	$15$비트	$16\,383$	$64$비트
`__float128` ($128$비트)	$15$비트	$16\,383$	$112$비트

자릿수가 제한되어 있는 만큼 이산적입니다. 예를 들어 $3$자리만 사용할 수 있는 상황이고 $A = \mathit{bias}$, 즉 $q = 2^0 \times \left(1+B\right) = 1+B$라고 한다면 $1.\underline{000}_2=1$과 $1.\underline{001}_2=1.125$ 사이의 수는 표현할 수 없게 되는 것입니다.

사실 $A$나 $\mathit{bias}$보다 중요하게 봐야 할 건 $B$의 유효숫자 수입니다. 이게 전부 $2$진수라서 생각하기 어려울 수 있으니 $\log_2 10$로 나눠서 생각해 보면 이렇습니다.

	$B$	$10$진 자릿수
`float` ($32$비트)	$23$비트	$6.92$자리
`double` ($64$비트)	$52$비트	$15.65$자리
`long double` ($80$비트)	$64$비트	$19.26$자리
`__float128` ($128$비트)	$112$비트	$33.71$자리

이 이후부터는 참값과 근삿값을 구분해야 할 일이 있으면 다음과 같은 표현법을 쓰도록 하겠습니다. $\F$는 ECMAScript 스펙에서 따왔습니다.

\[\begin{aligned} x_\R &= \textrm{참값 }x\\ x_\F &= \textrm{부동소수점 값으로 표현된 }x\textrm{의 근삿값} \end{aligned}\]

그리고 편의를 위해 $\mathit{bias}$는 $A$에 이미 더해져 있다고 생각하겠습니다.

정확도의 소실

유효숫자가 정해져 있다 보니 불가피하게 정확도의 소실이 일어납니다. 그래서 일반적으로 실수를 다루는 문제들은 절대/상대 오차 얼마 이하는 정답으로 간주합니다. 하지만 문제를 만드는 입장에서는 정해의 오차가 얼마 정도인지를 정확히 알고 있어야 참가자에게 ‘이 정도의 정확성을 갖는 풀이를 제출하십사’라고 요구할 수 있습니다. 반대로는 ‘이 정도의 정확도를 요구하려면 정해가 어느 정도로 정확해야겠구나’를 생각할 수 있겠죠.

여기서 ‘절대 오차’와 ‘상대 오차’에 대해 먼저 설명하도록 하겠습니다.

$x$의 근삿값 $y$가 $\delta_a$만큼의 ‘절대 오차’를 갖고 있다는 것은 $|x-y| \le \delta_a$임을 의미합니다.
반면 $y$가 $\delta_r$만큼의 ‘상대 오차’를 갖고 있다는 것은 $\frac{|x-y|}{x} \le \delta_r$임을 의미합니다. 실생활에서 ‘$1$% 정도의 차이가 발생할 수 있다’같은 표현이 상대 오차 $0.01$이 발생할 수 있다는 것과 같은 표현이겠네요.

이제 자주 사용하는 연산들에서 오차가 얼마나 발생할 수 있는지 하나하나 뜯어보도록 합시다.

입력

조금 전에 참값을 어떻게 float이 다룰 수 있는 값으로 근사하는지 설명했습니다. 참값 $x_\R$를 유효숫자 $p$비트의 float 값 $x_\F$로 나타내려고 한다면 정확히 얼마만큼의 오차가 발생할까요?

일단 $x_\R = 2^A \times \left(1+B_\R\right)$라고 한다면, $B$의 자릿수 제한을 맞춰 주기 위해 $B_\R$을 최대한 가까운 $p$자리 이진 소수로 반올림합니다. 이렇게 하면 $B$에는 최대 $0.5 \times 2^{-p} = 2^{-p-1}$만큼의 절대 오차가 발생하게 됩니다.

$x_\R$과 $x_\F$의 입장에서는 어떨까요? 앞서 언급한 것은 $B$에서의 오차였으므로, $x$ 입장에서 생각하려면 $2^A$를 곱해 줘야 합니다. 그러면 $x$ 입장에서의 절대 오차는 $2^A \times 2^{-p-1}$입니다.

앞서 $A$는 수의 크기에 따라 조정된다고 설명했습니다. 음? 입력받는 값의 스케일에 따라 결정된다구요? 그러면 오차를 절대 오차보다는 상대 오차로 나타내는 것이 괜찮은 옵션인 것 같아 보입니다.

상대 오차는 최대

\[ |x_\R-x_\F| \le 2^A \times 2^{-p-1} \quad \textrm{이므로} \quad \frac{|x_\R-x_\F|}{x_\R} \le \frac{2^A \times 2^{-p-1}}{2^A \times \left(1+B\right)} = \frac{2^{-p-1}}{1+B} \]

이고, $0 \le B < 1$이므로

\[\frac{2^{-p-1}}{1+B} \le \frac{2^{-p-1}}{1+0} = 2^{-p-1} = \frac{2^{-p}}{2}\]

가 됩니다. 따라서 반올림 오차는 상대 오차로 최대 $\frac{2^{-p}}{2}$입니다. 이를 간단하게 나타내기 위해, 부동소수점 수에서 다음 수와의 간격을 나타내는 기계 입실론^{machine epsilon}을 도입해 표기해 봅시다. 기계 입실론은 $\eps$과 같이 나타내고, $\eps = 2^{-p}$이므로 $\frac{2^{-p}}{2} = 0.5 \eps$로 표현할 수 있겠습니다.

결론: 실수를 입력받으면 최대 $0.5 \eps$만큼의 상대 오차가 발생합니다.

double에서 $p=52$임을 생각하면 유효숫자가 $15$자리쯤 된다고 볼 수 있습니다. 따라서 double을 입력받았다면 앞의 $15$자리까지는 믿을 수 있다고 할 수 있겠습니다.

부호가 같은 수의 덧셈

유효숫자 세 자리 float에서, $a_\F=1.\underline{101}_2$와 $b_\F=1\underline{10}.\underline{1}_2$를 더하는 상황을 가정해 봅시다.

그냥 더하면 $\left(a_\F + b_\F\right)_\R = 1\,\underline{000}.001_2$이 되지만, 우리의 float은 유효숫자 $3$개까지만을 허용하므로 $.001_2$ 부분을 날려야 합니다. 따라서 계산된 값을 다시 float에 저장하면 $\left(a_\F + b_\F\right)_\F = 1\,\underline{000}_2$ 부분만이 남습니다.

컴퓨터는 계산 과정에서는 저장하는 float의 두 배 정도의 공간을 사용해 정확하게 계산하고 이후 반올림하기 때문에, 덧셈의 경우에도 반올림 오차 한 번만이 발생합니다. 이는 상대 오차 $0.5 \eps$입니다.

부호가 같은 두 float의 덧셈에서는 $0.5 \eps$만큼의 상대 오차가 발생합니다.

그러나 여기에는 함정이 있는데, $a_\F \ne a_\R$이라는 점입니다. $a_\F$와 $a_\R$는 이미 $0.5 \eps$만큼의 상대 오차를 갖고 있습니다. 따라서 $\left(a_\F + b_\F\right)_\R$와 $a_\R + b_\R$의 상대 오차는 이미 $0.5 \eps$입니다. 여기서 반올림도 수행해야 하므로, 발생할 수 있는 오차는 최대

\[\left(1+0.5 \eps\right)\left(1+0.5 \eps\right)-1 = \eps + 0.25 \eps^2\]

입니다. $0.25 \eps^2$는 너무 작아서 무시할 수 있습니다(후술). 따라서 다시 결론을 내리면,

결론: 부호가 같은 두 실수를 입력받아 더하면 $\eps$만큼의 상대 오차가 발생합니다.

여담으로 큰 차이가 나는 두 수를 더하면 더하는 것에 의미가 없어질 수도 있습니다. 예를 들어 $a_\F=1\underline{001}_2$과 $b_\F=0.001\underline{101}_2$을 더하는 경우 계산 결과는 $\left(a_\F + b_\F\right)_\R=1\underline{001}.001101_2$이 되지만, 반올림하면 $1\underline{001}_2$이 되어 $a_\F$와 같아집니다. $\eps^2$ 자체는 float으로 정확히 표현할 수 있을지 몰라도, $\eps \times 2$ 이상의 값과 더하면 이와 같은 현상이 나타나기 때문에 생각하지 않아도 괜찮습니다.

곱셈과 나눗셈

이번에는 $a_\F=1.\underline{110}_2$와 $b_\F=1\underline{10}.\underline{1}_2$를 곱하거나 나누는 상황을 생각해 봅시다.

\[ \begin{aligned} \left(a_\F \times b_\F\right)_\R &= 1\underline{0}.\underline{11}0110_2 \\ \left(a_\F / b_\F\right)_\R &= 1.\underline{000}100111011000\cdots_2 \end{aligned} \]

반올림하면 이렇게 됩니다.

\[ \begin{aligned} \left(a_\F \times b_\F\right)_\F &= 1\underline{0}.\underline{11}_2 \\ \left(a_\F / b_\F\right)_\F &= 1.\underline{001}_2 \end{aligned} \]

계산 과정에서의 정확도 손실은 발생하지 않으며(나눗셈의 경우 무시할 수 있을 정도로 작습니다), 반올림 오차 $0.5 \eps$가 발생합니다.

이번에도 $a_\F \ne a_\R$임을 고려해 봅시다. $\left(a_\F\times b_\F\right)_\R$와 $a_\R\times b_\R$의 상대 오차는 $\eps$가 되고, 반올림을 수행해 $\left(a_\F\times b_\F\right)_\F$로 만들면 $1.5 \eps$가 됩니다. 나눗셈도 같은 방법으로 계산 가능합니다.

결론: 두 float의 곱셈/나눗셈에서는 $0.5 \eps$만큼의 상대 오차가 발생합니다.
두 실수를 입력받아 곱하거나 나누면 $1.5 \eps$만큼의 상대 오차가 발생합니다.

부호가 같은 수의 뺄셈

뺄셈, 혹은 부호가 다른 수의 덧셈은 위에서 언급한 덧셈과 비슷하지만 상황이 조금 다릅니다. 특히 $x_\F \ne x_\R$이면서 차이가 얼마 안 나는 값을 계산하려고 하는 경우입니다.

예를 들어 $a_\F=1.\underline{111}_2$에서 $b_\F=1.\underline{110}_2$을 빼는 경우를 생각해 봅시다. 결과는 $\left(a_\F-b_\F\right)_\R=0.001\underline{0\,00}_2$입니다. 이는 별 문제가 없어 보이지만, 계산 결과에서의 밑줄 친 $\underline{000}$은 $a_\F$에도 $b_\F$에도 없었던 값입니다. 이 새롭게 생긴 유효 비트 $3$개를 그냥 믿기에는 조금 꺼림칙합니다. 다행히도 두 값이 근삿값이 아니라 참값이라면, 즉 $a_\F=a_\R$고 $b_\F=b_\R$라면 계산 결과는 정확하고 믿을 수 있습니다.

하지만 $a_\F$가 정확한 값이 아니라면, 즉 $a_\F$과 $a_\R$의 차이가 크다면 문제가 됩니다. 예를 들어 $1.93$에서 $1.69$를 빼는 상황을 생각해 본다면 정확한 결과는 $0.24$여야 하지만, $1.93_\F=1.\underline{111}_2$, $1.69_\F=1.\underline{110}_2$이므로 $\left(1.93_\F-1.69_\F\right)_\F=0.001\underline{0\,00}_2 = 0.125$가 됩니다. 무려 두 배나 차이가 나네요!

조금 더 극단적인 예로 $|\varepsilon_\F| < \eps$에 대해 다음과 같은 계산을 한다고 가정해 봅시다.

\[ \left[ (1_\F+\varepsilon_\F)_\F-(1_\F-\varepsilon_\F)_\F \right]_\F \]

$2 \varepsilon_\F$이 나와야 할 것처럼 보입니다. 그러나 $\varepsilon_\F$이 반올림 오차보다 작으므로 실제로는 아래와 같이 계산됩니다.

\[ (1_\F+\varepsilon_\F)_\F-(1_\F-\varepsilon_\F)_\F= 1_\F-1_\F = 0_\F \]

이는 무려 $100$%의 상대 오차입니다. 이런 현상을 재앙적인 소거^{catastrophic cancellation}라고도 합니다.³

따라서 뺄셈을 정확하게 하고자 할 경우 비슷한 두 값의 뺄셈을 최대한 피하도록 해야 합니다. 예를 들어 앞서 언급한 $(1+\varepsilon_\F)-(1-\varepsilon_\F)$의 경우, $(1-1)-(-\varepsilon_\F-\varepsilon_\F)$로 계산 순서를 약간 바꾸면 조금 더 정확한 값을 얻을 수 있을 것입니다.

결론: 부호가 같은 두 실수의 뺄셈을 하는 상황은 가급적 피하는 것이 좋습니다.

덧셈, 곱셈, 나눗셈은 상대적으로 훨씬 ‘안전한’ 연산이므로 가능한 경우 덧셈, 곱셈, 나눗셈만을 사용하도록 수정합니다. 일례로 실수 값의 누적 합을 구해야 한다면 대신 펜윅 혹은 세그먼트 트리를 사용할 수 있습니다.

여담으로, 같은 부호이더라도 두 배 이상 차이가 나는 값의 뺄셈은 $\eps$ 이하의 오차만이 발생함을 덧셈의 경우와 비슷한 방법으로 보일 수 있습니다.

뺄셈을 포기할 수 없다면

뺄셈의 오차를 다른 상황과 같이 생각할 수 없는 이유는 없는 유효숫자를 만들어야 하는 상황 때문입니다. 따라서 상대 오차를 계산하기는 어렵습니다. 하지만 절대 오차를 계산해 볼 수는 있습니다.

$a_\R$에서 $b_\R$을 빼는 연산을 생각해 봅시다.

우선 이 값들을 입력받는 과정에서 반올림 오차 $0.5 \eps$가 생깁니다. 이 오차는 상대 오차이므로, 절대 오차로 바꾸면 각각 $0.5 \eps a_\R$과 $0.5 \eps b_\R$이 됩니다. 따라서 $a_\R-b_\R$과 $\left(a_\F-b_\F\right)_\R$은 최대 $0.5 \eps \left(|a|_\R+|b|_\R\right)$의 절대 오차를 가집니다.
이제 $\left(a_\F-b_\F\right)_\R$을 $\left(a_\F-b_\F\right)_\F$로 근사합니다. 이 과정에서 반올림 오차 $0.5 \eps$가 한 번 더 생깁니다. $0.5 \eps \left(|a|_\R+|b|_\R\right)$의 절대 오차에 $0.5 \eps$의 상대 오차가 한 번 더 생기므로, 총 절대 오차는 $\left[0.5 \eps \left(|a|_\R+|b|_\R\right)\right]\left(1+0.5\eps\right)$입니다.
앞서 언급했던 것처럼 $0.25 \eps^2$는 $0.5 \eps$에 비하면 무시할 수 있는 수준이므로, $0.5 \eps \left(|a|_\R+|b|_\R\right)$만큼의 절대 오차가 발생한다고 할 수 있습니다.

따라서 다루는 수들의 스케일을 알고 있다면, 다행히도 발생하는 오차를 절대 오차로나마 가늠할 수 있습니다.

결론: 부호가 같은 두 실수 $a$와 $b$가 있을 때,
$a$에서 $b$를 빼면 최대 $0.5 \eps \left(|a|+|b|\right)$만큼의 절대 오차가 발생합니다.

정리해 보면 이렇습니다.

연산	부호가 같은 경우 오차	부호가 다른 경우 오차
입력/반올림	상대; $0.5 \eps$	–
$a_\R+b_\R$	상대; $\eps$	절대; $0.5 \eps \left(\|a\|+\|b\|\right)$
$a_\R-b_\R$	절대; $0.5 \eps \left(\|a\|+\|b\|\right)$	상대; $\eps$
$a_\R\times b_\R$	상대; $1.5 \eps$	상대; $1.5 \eps$
$a_\R / b_\R$	상대; $1.5 \eps$	상대; $1.5 \eps$

그러나 이 표에서 언급한 덧셈과 뺄셈에서의 상대 오차는 덧뺄셈 한 번에서 발생하는 오차이고, 덧셈을 계속 누적해나간다면 덧셈을 할 때마다 반올림 오차가 새로 발생하고, 이는 절대 오차로서 더해짐을 간과하면 안 됩니다. 프로그래밍 문제에서는 일반적으로 덧셈을 여러 번 하는 상황이 많기 때문에, 덧뺄셈은 아예 절대 오차로만 생각하는 편이 좋겠습니다. 이렇게 하면 부호가 같은 경우와 다른 경우에 차이가 없어지며, 다음과 같이 다시 한 번 정리할 수 있습니다.

연산	오차
입력/반올림	상대; $0.5 \eps$
$a_\R+b_\R$	절대; $0.5 \eps \left(\|a\|+\|b\|\right)$
$a_\R-b_\R$	절대; $0.5 \eps \left(\|a\|+\|b\|\right)$
$a_\R\times b_\R$	상대; $1.5 \eps$
$a_\R / b_\R$	상대; $1.5 \eps$

정해 코드의 오차 계산 및 문제 설계

문제에서 얼마만큼의 오차를 허용할지는 정해 코드의 최대 오차를 계산함으로서 정할 수 있습니다.

BOJ #1546 평균과 똑같은 문제를 조금 더 어렵게 세팅한다고 가정해 봅시다. 정해 코드를 이렇게 짰습니다.

double s = 0;

int n;
cin >> n;
while (n--) {
    double x;
    cin >> x;
    s += x;
}

cout << fixed << setprecision(12) << s / n << '\n';

일단 $s$의 오차는 덧셈 $N$번이 누적되므로 높게 잡아도 $N \times 0.5 \eps \left(s+s\right) = Ns\eps$ 이하의 절대 오차가 생깁니다. 마지막에 나눗셈 한 번을 하므로 상대 오차 $1.5\eps$가 한 번 더 생깁니다.

그러면 총 절대 오차는 $Ns\eps$, 총 상대 오차는

\[\begin{aligned} \left(1+\frac{Ns\eps}{s}\right)\left(1+1.5\eps\right)-1 &= \left(1+N\eps\right)\left(1+1.5\eps\right)-1 \\ &= N\eps+1.5\eps+1.5N\eps^2 \\ &= \left(N + 1.5\right)\eps \end{aligned}\]

가 됩니다.

허용 오차 설정하기

double로 풀리게 하고 싶다면 어떻게 해야 할까요? double에서는 $10^{-16} < \eps < 10^{-15}$이므로 넉넉하게 $\eps = 10^{-15}$로 두고 계산하면, 위 코드의 절대 오차는 $Ns \times 10^{-15}$, 상대 오차는 $\left(N + 1.5\right) \times 10^{-15}$가 나옵니다. 따라서 $N$의 값에 따라 오차 허용 범위가 달라지는데, $N = 10^5$로 설정한다면 상대 오차가 약 $10^{5-15}=10^{-10}$이 되므로 $10^{-9}$ 혹은 그보다 넉넉한 허용 범위를 잡으면 무리 없이 통과할 것입니다. (통상적으로 상대 오차 혹은 절대 오차 중 하나만 범위에 들어오면 맞도록 하기 때문에 그렇습니다.)

반대로 $N = 10^5$인 상황에서 float을 저격하고 싶다면 어떻게 해야 할까요? float에서는 $10^{-7} < \eps < 10^{-6}$이므로, $\eps = 10^{-7}$로 보수적이게 잡으면 상대 오차 $10^{5-7}=10^{-2}$라는 계산이 나옵니다. 따라서 허용 오차 범위를 $10^{-2}$로 잡으면 float을 사용한 코드가 통과하지 못할 가능성이 있고, $10^{-3}$ 혹은 그 이하로 잡으면 float이 통과하기 어려운 문제로 세팅할 수 있을 것입니다.

위의 예시에서도 그렇듯이 float의 경우 덧셈을 $100\,000$번만 하는데도 $10^{-2}$ 혹은 그 이상의 오차가 발생하기 때문에, float 코드는 부정확한 값을 도출하기 상당히 쉽습니다. 그래서 float을 굳이 통과시켜 주려고 의도하는 상황이 아니고서야 float으로 작성한 코드는 그냥 틀렸다고 생각하고, double으로 해결하는 경우만 신경써도 충분합니다.

In Practice

그러나 문제에서 요구하는 계산이 많아지고 복잡해지면 코드마다 사칙연산을 수행하는 횟수가 천차만별일 수 있습니다. 따라서 실수 오차를 신경쓰는 게 문제가 물어보는 주요 포인트가 아닌 이상 실제 세팅할 때는 위에서 언급한 표의 모든 상수를 $10$ 정도로 놓고 오차 허용 범위를 설계하게 됩니다. 앞의 상수는 계산 횟수에 비례하기 때문입니다.

연산	일반적으로 허용하는 오차
$a_\R+b_\R$	절대; $10 \eps \left(\|a\|+\|b\|\right)$
$a_\R-b_\R$	절대; $10 \eps \left(\|a\|+\|b\|\right)$
$a_\R\times b_\R$	상대; $10 \eps$
$a_\R / b_\R$	상대; $10 \eps$

굳이 $10$으로 두는 이유는 앞서 말했듯이 상수와 계산 횟수는 비례하기 때문에, 상수가 너무 커질 정도로 계산 횟수가 많아지면 차라리 시간 초과가 나는 편이 낫기 때문입니다.

데이터 제작

특히 데이터를 제작할 때 주의해야 하는 점이 있습니다. 바로 데이터 자체에 오차가 있는 경우입니다. 데이터도 컴퓨터로 계산하는 것이기에 일반적인 방법으로 만든다면 데이터에도 그만큼의 오차가 생길 수밖에 없고, 이렇게 오차가 있는 데이터를 기준으로 채점하면 오히려 맞아야 하는 답을 틀리다고 채점하게 될 수도 있습니다. 실수를 많이 다뤄보지 않은 세터가 꽤 흔하게 하는 실수입니다.

이를 방지할 수 있는 방법은 여러 가지가 있는데, 몇 가지를 소개하자면:

long double, __float128, Python Decimal, Java BigDecimal 등 double보다 훨씬 정확한 자료형을 통해 계산합니다.
가능한 경우 모든 계산을 정수로 수행합니다. 분모와 분자를 관리하는 분수 구조체를 만들거나, Python Fraction을 사용하는 것도 고려해 볼 수 있습니다.

정해 소스의 오차가 최대 얼마인지를 계산하는 것의 중요성을 강조하고 싶습니다. 개인적으로 double로 풀리게 하고 싶다면 데이터는 적어도 최대 $10^{-15}$만큼의 오차만을 갖도록 만드는 것을 권장합니다.

여담

사칙연산 이외의 함수들의 오차는 어떻게 구하나요?

IEEE 754는 사실 아래의 연산들에 대해 가장 가까운 실수를 계산할 수 있도록 하고 있습니다⁴. 의외로 제곱근이 들어가 있습니다.

사칙연산
제곱근
단일 곱셈 누산
나머지 연산
최대 및 최솟값

C는 이를 따라서, sqrt와 fmod는 반올림 오차, 즉 상대 오차 $0.5\eps$만을 가지도록 합니다. 단, 이는 인자들이 참값임을 가정한 오차임에 주의해야 합니다.

제곱근을 제외하고, GNU C에서 실수와 관련된 함수들의 구현과 오류의 대략적인 양은 문서에 명시되어 있습니다. 문서에서 명시된 ULP^{unit in the last place}는 이 글에서 언급한 기계 입실론과 같은 단위인데, 상대 오차입니다.

보통 채점 환경으로 주로 사용되는 x86_64에서의 주요 함수의 오차는 다음과 같습니다. 모두 계산 결과에 대한 상대 오차이므로, 입력값의 오차를 따로 고려해 줘야 함에 유의합니다.

sin, cos, tan, asin, acos, atan: $\eps$
pow, exp, expm1, log, log1p: $\eps$
hypot: $\eps$
sinh, cosh, tanh, asinh, acosh, atanh: $2\eps$
cbrt: $4\eps$

일부 함수들은 정의역이 제한되어 있으므로(예를 들어 acos는 $-1$과 $+1$ 사이의 입력만 허용합니다) 계산 과정에서 실수로 정의역을 벗어나는 값을 넣지 않도록 주의해야 합니다. 또한 아래 함수들은 특정 범위의 입출력을 피하거나, 입력 범위에 각별히 주의하는 것을 요합니다.

log: 출력값이 $0$에 가까울 것으로 예상하는 경우.
- $\log \left(x+1\right)$이 필요한 상황이라면 $x$의 값을 조금 더 정확히 입력할 수 있는 log1p 함수를 사용합니다: $\operatorname{log1p} x=\log \left(x+1\right)$.
- $x$가 작다면, $x+1$을 계산하는 데에서 이미 정확도의 손실이 일어납니다.
exp, pow: 출력값이 $1$에 가까울 것으로 예상되는 경우.
- $e^x-1$의 정확한 값이 필요한 상황이라면 마찬가지로 expm1 함수를 사용합니다: $\operatorname{expm1} x=e^x-1$.
fmod: 나눠지는 값이 나누는 값보다 너무 큰 경우.
- 인자가 참값이 아니라 반올림된 값일 경우 뺄셈과 비슷하게 재앙적 소거가 쉽게 발생할 수 있습니다.

여담으로, tan 함수의 경우 $\left(\frac{1}{2}+n\right)\pi$ 왼쪽에서는 $+\infty$로 발산하고 오른쪽에서는 $-\infty$로 발산하지만, 어떤 부동소수점 값도 $\left(\frac{1}{2}+n\right)\pi$를 정확히 표현할 수 없어서 우려하는 오류가 발생하지 않습니다.

일부 연산은 CPU가 직접 지원하고, CPU의 구현에 따라 오차가 다르게 발생하는 경우가 존재합니다. 예를 들어 exp10 함수는 CPU가 10의 $N$제곱 연산을 지원한다면 코드로 계산하려고 하지 않고 하드웨어로 __builtin_exp10을 통해 더 빠르게 계산하려고 할 것입니다.

glibc 소스에서도 상세 구현을 확인하고 오차를 가늠할 수 있습니다.

다른 방식으로 오차를 구하는 글도 있던데 여기서는 왜 이렇게 구하나요?

아마 $a \pm b$의 오차 $\delta \left(a \pm b\right)$는 $\sqrt{\left(\delta a\right)^2+\left(\delta b\right)^2}$라고 설명하는 글을 보셨을지도 모르겠습니다. 이는 물리학 등에서 $a$와 $b$가 여러 번의 측정을 통해 얻어낸 평균값이고, 오차가 확률분포(특히 정규분포)에서 추출될 것이라고 가정할 수 있을 때 사용하는 방법입니다. 즉, 우리가 다루는 오차와는 다른 성격의 오차입니다. 여기에서는 수치해석적 관점에서의 오차 전파를 다루고 있습니다.

마무리하면서

부동소수점 오차는 참 다루기 힘든 주제입니다. 이 글이 부동소수점 오차에 대한 막연한 가늠과 추측을 해소하고 좀 더 좋은 문제들을 만드는 데 기여해, 제가 백준에서 맞왜틀을 덜 당할 수 있기를 희망합니다.

이 글 작성에 도움을 주신 김영현^kipa00님과 김준원^junie님께 감사의 말씀을 전합니다.

각주

double의 경우 $0.1_\F$은 정확히 \[0.1000\,0000\,0000\,0000\,0555\,1115\,1231\,2578\,2702\,1181\,5834\,0454\,1015\,625\]입니다. ↑
$q=0$, $q=+\infty$, $q=\,$NaN, 그리고 $A = 0$인 경우는 예외입니다. ↑
반면 두 값이 참값인 경우에서 일어나는 현상은 온화한 소거^{benign cancellation}라고 합니다. ↑
https://en.wikipedia.org/wiki/IEEE_754#Required_operations ↑

(23/7/19) 사칙연산 이외 함수들에 대한 설명을 추가했습니다. kiwiyou님께 감사의 말씀을 드립니다.

모든 캐릭터를 모으기 위해서는 얼마를 써야 하는가

가챠(뽑기)는 나쁜 문명입니다.

뜬금없이 웬 가챠냐고 하면, 요즘 게임을 ‘확률성 아이템’을 빼놓고는 이야기할 수 없을 정도이기 때문입니다. 물론 필자가 하고 있는 게임 대부분에도 확률성 아이템이 들어가 있습니다. 하지만 저는 확률성 아이템의 윤리성 같은 것을 논하는 현명한 소비자와는 꽤 거리가 있으므로, 전혀 반대의 논지, 즉 몇 가지 게임에서 ‘모든 캐릭터를 모으려면 얼마를 써야 할까?’(Shut up and take my money)를 중심으로 이야기해보겠습니다.

확률 계산의 기본

통계학에서 확률을 계산할 때 수학적인 표현을 위해 사용되는 개념들이 있습니다. 확률변수, 확률질량함수, 기댓값이 그들 중 하나입니다.

확률변수는 값과 확률을 갖는 수입니다. 보통 대문자로 표시합니다. 예를 들어 주사위를 던져서 나오는 수를 확률변수 $X$라고 할 수 있습니다.

확률질량함수는 확률변수와 그 값을 대응시켜 주는 함수입니다. $P\left(X=x\right)$과 같은 방식으로 표기합니다. 쉽게 말하자면, $P\left(X=x\right)$는 $X$가 $x$일 확률을 말합니다. 위에서 언급한 주사위를 던져서 나오는 수를 예로 들면,
\[P\left(X=x\right) = \dfrac{1}{6} \left(x = 1, 2, 3, 4, 5, 6\right)\]
이라고 표시할 수 있겠습니다.

그리고 기댓값은 말 그대로 확률변수에 기대하는 값입니다. 평균이라고도 하는 기댓값은 $E\left(X\right)$로 표기합니다. 주사위의 눈 같은 경우엔 평균 $E\left(X\right) = 3.5$입니다.

그러나 주사위 안쪽에 쇠 같은 걸 붙여 둬서 한 면만 나올 확률이 비정상적으로 높은 주사위가 있다고 생각해 봅시다. 예를 들어, 6은 나올 확률이 $\dfrac{95}{100}$고, 나머지는 나올 확률이 $\dfrac{1}{100}$밖에 안 된다고 생각해 봅시다. 이 때도 주사위를 많이 던졌을 때 평균 $E\left(X\right) = 3.5$라고 생각할 수 있을까요?

아닙니다. 이 때는 거의 6이 나올 것이므로 평균이 6에 가까우면 가까웠지 3.5에 가깝진 않을 겁니다. 그래서 평균, 즉 기댓값을 구할 때는 평균과 확률을 곱해 줘야 합니다. 이 때 기댓값은 이렇게 구할 수 있습니다.

\[E\left(X\right) = 1\times\dfrac{1}{100} + 2\times\dfrac{1}{100} + 3\times\dfrac{1}{100} + 4\times\dfrac{1}{100} + 5\times\dfrac{1}{100} + 6\times\dfrac{95}{100} = 5.85\]

5.85는 3.5보다는 조금 믿을 수 있는 숫자 같네요!

아무튼 이걸 일반화하면, $P\left(X=x_i\right) = p_i \left(i=1, 2, \ldots, n\right)$일 때 아래와 같이 적을 수 있습니다.
\[E\left(X\right) =\sum_{i=1}^{n} {x_i}{p_i}\]
또한 보이다시피 $E\left(X\right)$는 ${x_i}{p_i}$들의 합으로 이루어져 있기 때문에, $E\left(aX+b\right)=aE\left(X\right)+b$, $\sum E\left(X\right)=E\left(\sum_{}^{} X\right)$ 등도 성립합니다.

모든 카드를 얻으려면 – 등급별 가중치가 없을 때

여러 가지 게임들의 예를 들어 보겠습니다만, 일단 ‘사운드 볼텍스’의 예를 들어 보겠습니다. 사운드 볼텍스는 리듬게임이지만 스토리를 진행함에 있어 높은 등급의 카드를 뽑는 것이 거의 필수적입니다. 등급이 나누어져 있긴 하고 실제로 중복 카드를 일부러 더 잘 나오게 설정했다는 게 체감된다고 하지만, 일단은 실제 확률이 고지되어 있지 않으므로 등급별 확률이 정해져 있지 않고 모든 카드가 나올 확률이 같다고 가정하겠습니다.

‘사운드 볼텍스 IV’의 PUR 등급 카드 중 하나의 모습. 필자도 갖고 있습니다.

일단 몇 가지를 가정해 봅시다.

카드는 한 번 뽑는 데 1,000원입니다.
카드는 한 장 씩 뽑고, 다음 카드를 뽑기 전에 카드를 확인합니다.
카드는 전체 $n$종류입니다. 각 카드 종류가 등장할 확률은 $\dfrac{1}{n}$입니다.
$n$종류의 카드를 전부 뽑으면 더 이상 뽑지 않습니다.
게임에서는 이미 나온 카드를 10번 다시 뽑으면 다음 한 장은 중복이 아님이 보장되지만, 여기서는 무시합니다.

이 때 우리는 평균적으로 카드를 총 몇 장 뽑아야 모든 종류의 카드를 뽑을 수 있을까요? 즉 뽑아야 하는 카드 수를 확률변수 $X$라고 둘 때 $\left(X\right)$는 얼마일까요?

그냥 구하고자 하니 막막합니다. 조금 더 쉬운 방법을 생각해 봅시다. 이미 $i-1$종류의 카드를 모았을 때 새로운 종류의 카드를 모으는 데 뽑는 횟수를 확률변수 $Y_i$라고 두면 $X = Y_1 + Y_2 + Y_3 + \ldots + Y_n$으로 생각할 수 있습니다.

이렇게 하면 $P\left(Y_i = k\right)$는 $i-1$종류를 뽑은 상황에서 $k$번째에 새로운 종류의 카드를 뽑은 셈이 되므로, $k-1$번은 원래 뽑았던 카드 $i-1$종류 중 하나를, 마지막에는 남은 카드 종류 $n-\left(i-1\right)$개 중 하나를 뽑는 셈이 됩니다. 즉,

\[P\left(Y_i = k\right) = {\left(\dfrac{i-1}{n}\right)}^{k-1}{\left(1-\dfrac{i-1}{n}\right)}\]

가 됩니다. 식이 나왔으니 카드 $i-1$종류가 있을 때 새로운 종류의 카드를 뽑으려면 몇 장을 뽑아야 하는지에 대한 기댓값을 구할 수 있는데요, 이론적으로는 무한 장 뽑아도 새 종류가 안 나올 수도 있으니 $k$에는 사실상 1부터 무한대까지 들어갈 수 있겠습니다. 그러므로,

\[E\left(Y_i\right) =\sum_{k=1}^{\infty} kP\left(Y_i = k\right)\]
\[=\sum_{k=1}^{\infty} k{\left(\dfrac{i-1}{n}\right)}^{k-1}{\left(1-\dfrac{i-1}{n}\right)}\]

식이 복잡하므로 $\dfrac{i-1}{n} = r$으로 치환하고 계속 계산해 보면,

\[=\sum_{k=1}^{\infty} k{r^{k-1}{\left(1-r\right)}}\]
\[=\left(1-r\right)\sum_{k=1}^{\infty} kr^{k-1}\]

여기서 $r$이 1보다 작으므로 위의 멱급수는 무한등비급수의 합을 이용해 계산할 수 있습니다. 항등식 $\sum_{k=1}^{\infty} r^{k} = \dfrac{1}{1-r}$에서 양변을 $r$에 대해 미분하면 $\sum_{k=1}^{\infty} kr^{k-1} = \dfrac{1}{{\left(1-r\right)}^2}$이므로, 이를 이용해서 정리하면

\[E\left(Y_i\right) = \left(1-r\right)\dfrac{1}{{\left(1-r\right)}^2}\]
\[= \dfrac{1}{1-r}\]

라는 꽤 간단한 식이 됩니다. 최종적으로 $r = \dfrac{i-1}{n}$을 다시 대입하고 정리해 주면

\[E\left(Y_i\right) = \dfrac{n}{n-i+1}\]

입니다.

처음에 $X$를 $n$개의 $Y_i$로 쪼갰으니까 $X = \sum_{i=1}^{n} Y_i$라고 쓸 수 있겠습니다. 앞서 언급한 기댓값의 정의에 의해

\[E\left(X\right) = E\left(\sum_{i=1}^{n} Y_i\right) = \sum_{i=1}^{n} E\left(Y_i\right)\]

이므로, $E\left(Y_i\right) = \dfrac{n}{n-i+1}$를 대입하면

\[E\left(X\right) = \sum_{i=1}^{n} \dfrac{n}{n-i+1}\]

입니다. 여기서 $u = n-i+1$, 즉 $i = n-u+1$로 치환해 버린다면 1에서부터 $n$까지의 $i$의 범위는 $n$에서부터 1까지의 $u$의 범위가 되고, 시그마는 범위가 반대여도 결과는 같으므로 식은

\[E\left(X\right) = \sum_{u=n}^{1} \dfrac{n}{u} = n\sum_{u=1}^{n} \dfrac{1}{u}\]

이 나옵니다. $\sum_{u=1}^{n} \dfrac{1}{u}$는 달리 간단히 할 방법이 없으므로 여기서 만족하고 Wolfram|Alpha에 맡기도록 합시다.

적용

이제 식에 대입하는 일만 남았습니다! 만약 카드 한 세트가 총 31장이라고 하고 31장 안에서만 카드가 나온다고 한다면, $n = 31$일 때 $E\left(X\right)$는

\[E\left(X\right) =31\sum_{u=1}^{31} \dfrac{1}{u} \approx 124.84460\]

이므로, 약 124,845원을 넣으면 모든 카드를 수집할 수 있다는 뜻이 되겠습니다! 물론 실제 게임에서는 이미 나온 카드를 10번 다시 뽑으면 다음 한 장은 중복이 아님이 보장되지만 이를 무시하고 계산한 것이기에 숫자가 현저히 크게 나오는 감도 없지 않습니다.

모든 아이돌을 프로듀스하려면 – 등급별 가중치가 있을 때

그러나 많은 게임들은 그렇게 호락호락하지 않습니다. 바로 등급에도 확률을 부여하기 때문입니다.

대표적으로 ‘아이돌 마스터 신데렐라 걸즈 스타라이트 스테이지’는 노멀, 레어, S레어, SS레어의 4개 등급의 카드를 만들어 놓았습니다. 5명의 아이돌을 유닛으로 편성해 즐기는 리듬게임인데, SS레어 아이돌들이 점수와 직결되는 합계 어필 수치가 현저히 높고 부가 능력치도 훨씬 좋기 때문에 S레어 5인으로 덱을 구성해 게임을 하는 것보다 SS레어 5인으로 구성하는 것이 훨씬 점수가 잘 나와서 높은 점수를 얻으려면 SS레어 아이돌의 존재는 필수적입니다. 결정적으로 일러스트가 상당히 예뻐서 수집욕구를 불러일으키기도 합니다.

하지만 이 게임은 SSR를 얻을 수 있는 확률은 1.5%라고 공지해 두고 있습니다. 결국 98.5%는 S레어 이하의 등급이라는 뜻이죠. 그러면 모든 종류의 SSR 아이돌을 얻기 위해서는 얼마만큼의 돈이 필요할까요?

몇 가지를 가정하겠습니다.

아이돌은 한 번 스카우트하는 데 스타 쥬얼 250개가 필요합니다. 스타 쥬얼 0.85개는 1엔입니다.
아이돌은 한 명 씩 스카우트하고, 다음 아이돌을 스카우트하기 전에 아이돌의 이름을 기억해 둡니다.
SS레어 아이돌은 전체 $n$명입니다. 각 아이돌이 등장할 확률은 등급 내에서 $\dfrac{1}{n}$입니다. 한정 아이돌은 한정으로 생각하지 않습니다.
SS레어 아이돌이 등장할 확률은 $p$입니다.
$n$명의 SS레어 아이돌을 전부 스카우트하면 더 이상 스카우트하지 않습니다.

이번에도 사운드 볼텍스에서 계산했던 것과 같이 확률변수 $X$를 $n$개의 $Y_i$로 쪼개서 계산합시다. 하지만 이번에는 SSR가 아닌 아이돌은 고려하지 않습니다. 저번처럼 $P\left(Y_i = k\right)$는 $i-1$명을 스카우트한 상황에서 $k$번째에 새로운 아이돌을 스카우트한 것인데, $k-1$번은 원래 스카우트했던 SSR 아이돌 $i-1$명 중 하나 또는 SR 이하 아이돌을, 마지막에는 SSR 등급에 아직 스카우트하지 못한 아이돌 $n-\left(i-1\right)$명 중 하나를 뽑는 셈이 됩니다.

여기서 SSR 등급에 남은 아이돌 $n-\left(i-1\right)$명 중 한 명을 스카우트하는 확률은 $p\left(1-\dfrac{i-1}{n}\right)$입니다. 원래 스카우트했던 SSR 아이돌 $i-1$명 중 하나 또는 SR 이하 아이돌을 스카우트하는 경우는 전체 카드 중 SSR 등급에 아직 스카우트하지 못한 아이돌 $n-\left(i-1\right)$명 중 한 명을 제외하고 누구를 스카우트하든 상관없는 경우이므로, 확률은 $1 – p\left(1-\dfrac{i-1}{n}\right)$이라고 할 수 있겠습니다. 따라서,

\[P\left(Y_i = k\right) = {\left(1 – p\left(1-\dfrac{i-1}{n}\right)\right)}^{k-1}{\left(p\left(1-\dfrac{i-1}{n}\right)\right)}\]

계산하기 쉽게 $p\left(1-\dfrac{i-1}{n}\right) = s$로 놓으면

\[P\left(Y_i = k\right) = {\left(1 – s\right)}^{k-1}{\left(s\right)}\]

$1-s=r$로 놓으면

\[P\left(Y_i = k\right) = r^{k-1}{\left(1-r\right)}\]

으로 지난번과 똑같은 형태가 됩니다. 다만 $r$ 안의 식이 지난번과 다를 뿐입니다.

역시 이론적으로는 무한 번 뽑아도 안 나올 수 있으므로 지난번과 같은 계산을 하면

\[E\left(Y_i\right) = \dfrac{1}{1-r} = \dfrac{1}{s}\]
\[= \dfrac{n}{p\left(n-i+1\right)}\]

이 됩니다. $E\left(X\right)$도 계산해 보면
\[E\left(X\right) = \sum_{i=1}^{n}\dfrac{n}{p\left(n-i+1\right)}\]
\[= \dfrac{1}{p}n\sum_{i=1}^{n}\dfrac{1}{\left(n-i+1\right)}\]
이므로, 앞서 $\sum_{i=1}^{n}\dfrac{1}{\left(n-i+1\right)} = \sum_{u=1}^{n}\dfrac{1}{u}$로 계산한 것과 같이

\[E\left(X\right) = \dfrac{n}{p}\sum_{u=1}^{n}\dfrac{1}{u}\]

가 됩니다. 따라서 기댓값은 등장 확률에 반비례함을 알 수 있습니다.

적용

데레스테 SSR 소유율이라는 사이트에 의하면 현재(2017년 8월 13일) SSR 아이돌은 총 127명입니다. (통상 70명, 한정 57명) 또한, SSR 확률은 1.5%이므로 $p = \dfrac{3}{200}$입니다.

한정 아이돌을 고려하지 않고 127명이 가챠에서 모두 나올 수 있다고 가정하면, $n = 127$일 때

\[E\left(X\right) = \dfrac{200\times 127}{3}\sum_{u=1}^{127}\dfrac{1}{u}\]

\[\approx\dfrac{200\times 127}{3}\times 5.4253346\approx45934.5\]

이므로, 평균적으로 스타 쥬얼을 11,483,625(1148만 3625)개, 즉 구매량에 따른 효율에 따라 평균적으로 13,353,052엔(1335만 3052엔, $\approx$ 1.35억원)을 과금하면 SSR 아이돌을 전부 모을 수 있습니다.

반면 통상 아이돌 70명만 고려한다면, $n = 70$일 때

\[E\left(X\right) = \dfrac{200\times 70}{3}\sum_{u=1}^{70}\dfrac{1}{u}\]

\[\approx\dfrac{200\times 70}{3}\times 4.8328368\approx22553.2\]

이므로, 이 경우엔 스타 쥬얼 5,638,300개가 필요합니다. 그러므로 평균적으로 6,633,294엔(663만 3294엔, $\approx$ 0.670억원)을 과금하면 통상 SSR 아이돌을 전부 모을 수 있습니다.

신데렐라 페스 한정 아이돌, SSR 마에카와 미쿠. 한 아이돌은 보통 여러 카드에 등장합니다.

6월 말에 신데렐라 페스를 진행했을 때 SSR 확률이 3%로 고정되었습니다. 이 때의 $p = \dfrac{3}{100}$이므로 페스 기간 동안에만 가챠를 돌린다면, 소모되는 비용은 등장 확률에 반비례하기 때문에 현재는 위에 적힌 가격의 반값으로 전부 스카우트하는 것을 기대해 볼 수 있습니다.

스타 쥬얼을 얻을 수 있는 경로가 필요한 스타 쥬얼의 양에 비해 꽤 적은 비율이지만 과금 외에도 존재하므로 이를 고려한다면 위에 적힌 가격보다는 적을 수도 있습니다.

모든 용사를 얻으려면 – 등급별 가중치와 결제 한정 캐릭터를 고려할 때

용사에게 등급이 있지만 4성 이후로는 등급을 올리는 게 가능해 오히려 도감을 채우기 위해 5성 이상보다는 4성이 선호되는 게임 ‘크루세이더 퀘스트’

‘크루세이더 퀘스트’는 특이한 게임입니다. 뽑기로 얻을 수 없는 ‘전설 용사’를 제외하고 1~3성 용사는 ‘일반 용사’, 4~6성은 ‘고급 용사’로 분류됩니다. 용사에 등급이 매겨져 있는 건 다른 게임과 다를 게 없으나, 용사들은 조금의 노력으로 더 높은 단계로 승급할 수 있습니다. 2~4성으로의 승급은 아예 랜덤이지만 5~6성으로의 승급은 고정됩니다. 예를 들어 3성 가챠레인저B는 4성 바이올렛으로 승급이 가능하지만, 4성 바이올렛이 승급하면 5성 인형사 바이올렛, 6성 마법인형사 바이올렛이 고정되는 식입니다.

하지만 고급 용사는 또 ‘승급 용사’와 ‘계약 전용 용사’로 나뉩니다. 계약 전용 용사는 보석을 5~6개 소모하는 ‘황금 계약서’에서 확률적으로 나옵니다. 그러면 모든 ‘계약 전용 용사’를 얻기 위해서는 얼마가 필요할까요?

이번의 조건은 다음과 같습니다.

용사는 연속으로 계약한다고 가정합니다. 첫 번째 계약 가격을 무시하고, 한 번 계약하는 데 보석 $5$개를 소모합니다. 보석 1개는 550원입니다.
용사는 한 명 씩 계약하고, 다음 용사를 계약하기 전에 용사의 이름을 확인합니다.
4성 이상의 용사는 전체 $n$명입니다. 각 용사가 등장할 확률은 $\dfrac{1}{n}$입니다. 이 중 한정 용사는 $m$명입니다.
$n$종류의 용사와 모두 계약하면 더 이상 계약하지 않습니다.
게임에서는 10명씩 계약하면 마지막 1명은 4성 이상이 보장되지만, 여기에서는 고려하지 않습니다.
게임에서는 기간 한정으로 특정 용사의 계약확률이 증가하지만, 여기에서는 고려하지 않습니다.

앞서 확률이 $p$인 등급의 아이돌 $n$명을 전부 스카우트할 때의 기댓값은

\[E\left(X\right) = \dfrac{n}{p}\sum_{u=1}^{n}\dfrac{1}{u}\]

이었습니다. 여기서 ‘모든 계약 전용 용사’를 4~6성 내의 또 하나의 등급으로 생각한다면, 4~6성 내에서 계약 전용 용사와 계약할 확률은 $\dfrac{m}{n}$이므로 더 복잡한 계산을 할 필요 없이 모든 계약 전용 용사와 계약했을 때 계약 횟수에 대한 기댓값은

\[E\left(X\right) = \dfrac{n}{\dfrac{m}{n}p}\sum_{u=1}^{m}\dfrac{1}{u}\]

\[= \dfrac{n^2}{mp}\sum_{u=1}^{m}\dfrac{1}{u}\]

이 됩니다.

적용

황금 계약서에서 4성 이상의 용사가 나올 확률은 19%입니다. 또한 용사등급표에 따르면 현재(2017년 8월 13일) 크루세이더 퀘스트에는 계약서로 획득할 수 없는 이벤트, 전설 용사를 제외하고 승급 용사 43명, 계약 용사 58명이 있으므로 $n = 43 + 58 = 101$, $m = 58$, $p = \dfrac{19}{100}$으로 두면

\[E\left(X\right) = \dfrac{101^2 \times 100}{58 \times 19}\sum_{u=1}^{58}\dfrac{1}{u}\]

\[\approx \dfrac{101^2 \times 100}{58 \times 19} \times 4.646254 \approx 4300.95\]

이므로 평균적으로 보석을 21,505개, 즉 11,827,750원(1182만 7750원)을 과금하면 모든 계약 용사를 얻을 수 있습니다. 물론 이 게임은 보석을 획득할 수 있는 경로가 과금 이외에도 굉장히 다양하고 계약 용사 확정을 고려하지 않았으므로 실제 과금액은 이보다 현저히 적으리라 예상됩니다.

정리

확률형 아이템이 게임을 진행하는 데 중요한 요소가 되는 대표적인 게임 ‘아이돌 마스터 신데렐라 걸즈 스타라이트 스테이지’. 물론 여기다 쓸 돈을 다른 데다 쓰는 게 나을 수도 있겠지만, 적어도 자동차에 쓰면 자동차가 춤추고 노래하진 않을 것 같습니다.

계산하면서 큰 값이 나오리라고는 예상하고 있었지만, 예상을 훨씬 웃도는 숫자들이 나오면서 많이 놀랐습니다.

특히 ‘아이돌 마스터 신데렐라 걸즈 스타라이트 스테이지’의 경우 모든 아이돌을 모으는 데에 평균적으로 몇 천만원대의 과금이 필요하다는 것을 계산하고 나서는 차라리 자동차를 사는 게 낫다고 생각될 수도 있겠지만, 역시 자동차가 춤추고 노래하진 않기에 나름의 의미가 있다고 생각합니다.

마지막으로, 언급된 모든 게임들이 캐릭터를 모두 모으지는 않아도 진행하는 데에 큰 무리는 없으므로, 언급된 가격들은 전부 모으실 게 아니라면 ‘이렇게 큰 숫자가 나오는구나’ 라고 생각해 주시길 부탁드리겠습니다.

컴퓨터는 삼각함수를 어떻게 계산하는가

주의: 이 포스트는 약간의 미적분학 지식이 있어야 이해하기 쉽습니다. 그래도 설명을 위해 약간의 증명은 들어가 있으니, 아는 부분은 스킵하셔도 됩니다.

우리는 언제나 직각과 직사각형이 편합니다. 인류는 직교좌표계를 발명했습니다. 그리고 컴퓨터 모니터의 픽셀 배열은 대부분 가로세로 수백~수천 개 픽셀로 이루어진 격자로 이루어져 있고, 이는 직교좌표계로 접근할 수 있습니다. 직교좌표계로 화면을 그리면 픽셀 하나하나를 관리하기가 제일 쉽고 자유로워서가 아닐까 싶습니다.

만약에 누군가가 사칙연산만 제공되는 엔진으로 게임을 개발한다고 생각합시다. 직교좌표계 위에 그려진 어떤 도형이 있는데 이 도형을 회전시켜야 한다고 합니다. 어떻게 하면 될까요?

사실 수학은 이미 답을 알고 있습니다. 어떤 픽셀의 좌표 $\left ( x, y\right )$에 대해 회전된 좌표 $\left ( x′, y′\right )$는
\[\begin{bmatrix}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}x′\\ y′\end{bmatrix}\]
다른 말로 \[\begin{cases}x′=x \cos \theta-y\sin \theta\\ y′=x \sin \theta+y \cos \theta \end{cases}\] 임이 자명합니다. 이걸 회전하기 전의 이미지의 모든 픽셀의 $\left ( x, y\right )$에 대해 한 번씩 해 주면 회전한 이미지가 나옵니다. 참 쉽죠?

사칙연산만 갖고 삼각함수를 계산하라고요?

컴퓨터가 어떻게 삼각함수를 계산하는지를 논하기로 했으니까, 컴퓨터가 할 수 있는 연산만을 이용해야겠습니다. 아주 기본적인 연산만 사용 가능한 프로세서에서 말입니다.

그럼 대체 삼각함수는 어떻게 근사할 수 있을까요? $\sin$과 $\cos$를 구하면 $\tan$ 등은 나눗셈으로도 구할 수 있으니까 $\sin$과 $\cos$를 구하는 데 집중해 봅시다.

삼각함수를 다항함수로 근사하면 되지 않을까요?

다항함수는 사칙연산만으로 계산할 수 있으니 삼각함수의 그래프와 비슷한 다항함수를 만들어서 거기다 집어넣으면 되겠네요!

근데 그 다항함수는 대체 어떻게 구할까요? 여기서 테일러 급수가 등장합니다.

테일러 급수

테일러 급수를 이용하면 어떤 미분 가능한 함수라도 다항함수로 근사해 버릴 수 있습니다. 세상에 그런 흑마법이 존재하냐고요? 네, 놀랍게도 존재합니다!

증명

우선 어떤 함수 $f^\prime\left(x\right )$를 $a$부터 $x$까지 정적분해 봅시다. 미적분의 정의에 의해 아래 식과 같이 표현할 수 있습니다.
\[\int_{a}^{x}f^\prime\left(t\right )\mathrm{d}t=f\left(x\right )-f\left(a\right )\]
그리고 위 식을 변형해 봅시다.
\[\int_{a}^{x}f^\prime\left(t\right )\mathrm{d}t=\int_{a}^{x}\left(-1\right )\left(-f^\prime\left(t\right )\right )\mathrm{d}t\]
이제 $-1$을 적분하고 $-f^\prime\left(t\right )$를 미분해 부분적분법을 씁시다. 이 때 $f\left(x\right )$가 무한히 미분 가능하다면, 부분적분법도 무한히 써 버릴 수 있습니다. 그러니까
\[\cdots = \left.\begin{matrix}\left (-\left ( x-t \right )f^\prime\left ( t \right )- \dfrac{\left( x-t \right )^2}{2}f^{\prime\prime}\left ( t \right )- \dfrac{\left( x-t \right )^3}{6}f^{\prime\prime\prime}\left ( t \right )- \cdots\right )\end{matrix}\right|_{a}^{x}\]
가 됩니다. 이 때 적분변수 $t$와 관계없는 $x$는 상수취급됩니다. 이 식을 전개해 보면
\[\cdots =\left ( x-a \right )f^\prime\left ( a \right )+ \dfrac{\left( x-a \right )^2}{2!}f^{\prime\prime}\left ( a \right )+ \dfrac{\left( x-a \right)^3}{3!}f^{\prime\prime\prime}\left ( a \right )+ \cdots\\ \]
\[= f\left ( x \right )-f\left ( a \right )\]
마지막으로 $f\left ( a \right )$를 이항하면
\[f\left ( x \right ) =f\left ( a \right )+ \left ( x-a \right )f^\prime\left ( a \right )+ \dfrac{\left( x-a \right )^2}{2!}f^{\prime\prime}\left ( a \right )+ \dfrac{\left( x-a \right )^3}{3!}f^{\prime\prime\prime}\left ( a \right )+ \cdots\]
즉 $f^{\left( n \right)}$이 $f$를 $n$번 미분한 함수라고 할 때
\[f\left ( x \right ) = \sum_{n=0}^{\infty}\dfrac{f^{\left ( n \right )}\left ( a \right )}{n!}\left( x-a \right )^n\]
가 유도됩니다. 이 때 $a$에 어떤 수를 넣어도 근사됩니다.

그래서 이게 삼각함수랑 무슨 상관이에요!

삼각함수는 무한히 미분 가능합니다. 그러니까, $f$ 대신 삼각함수를 넣으면 삼각함수를 다항함수로 근사할 수 있습니다. 어떤 삼각함수 $f\left ( x \right )$에 대해 $a$에 $0$을 대입해 보면
\[f\left ( x \right ) =f\left ( 0 \right )+ \left ( x-0 \right )f^\prime\left ( 0 \right )+ \dfrac{\left( x-0 \right )^2}{2!}f^{\prime\prime}\left ( 0 \right )+ \dfrac{\left( x-0 \right )^3}{3!}f^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
\[=f\left ( 0 \right )+ x f^\prime\left ( 0 \right )+ \dfrac{x^2}{2!}f^{\prime\prime}\left ( 0 \right )+ \dfrac{x^3}{3!}f^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
인데요, $f\left ( x \right ) = \sin x$라고 하고 대입해 보면
\[\sin x =\sin\left ( 0 \right )+ x \sin^\prime\left ( 0 \right )+ \dfrac{x^2}{2!}\sin^{\prime\prime}\left ( 0 \right )+ \dfrac{x^3}{3!}\sin^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
\[=0+ \left ( x \cdot 1 \right )+ \left ( \dfrac{x^2}{2!}\cdot 0 \right )+ \left ( \dfrac{x^3}{3!}\cdot -1 \right )+ \left ( \dfrac{x^4}{4!}\cdot 0 \right )+ \left ( \dfrac{x^5}{5!}\cdot 1 \right )+ \cdots\]
정리하면
\[\therefore \sin x =x- \dfrac{x^3}{3!}+ \dfrac{x^5}{5!}- \dfrac{x^7}{7!}+ \dfrac{x^9}{9!}+ \cdots\]
마찬가지로
\[\cos x =1-\dfrac{x^2}{2!}+\dfrac{x^4}{4!}-\dfrac{x^6}{6!}+\dfrac{x^8}{8!}+ \cdots\]
이 됩니다. 중간 과정에 비해서 정말 간단한 식이 되었습니다. 더구나, 우리의 궁극적인 목표인 사칙연산만으로 삼각함수 표현하기에 성공했습니다! (팩토리얼은 곱셈의 연속일 뿐이니까요!)

근데, 식이 참… 무한합니다. 이건 어떻게 하면 좋을까요? 우리는 어차피 floating-point 변수들은 굉장히 정확하지 않다는 걸 잘 알고 있습니다. 그러니까 저 식을 어디까지만 계산하고 그 이후의 오차는 무시해 버려도 됩니다.

위에서 구한 다항식을 최고차항이 $n$인 항까지만 계산하고, 그 이후는 무시해버립시다. 그리고 이걸 $n$차 근사식이라고 합시다. 이제부터 차수를 올려나가면서 $n$차 근사식의 그래프를 그려보겠습니다.

그래프

▶ 1차 근사식: $f\left ( x \right ) = x$

… 하나도 비슷해보이진 않습니다. 다음!

▶ 3차 근사식: $f\left ( x \right ) = x-\dfrac{x^3}{3!}$

앞에서는 감을 좀 잠은 거 같기도 한데, 겨우 $\dfrac{\pi}{2}$도 가기 전에부터 조금씩 이상하더니 아래로 곤두박질쳐버리는군요. 다음!

▶ 5차 근사식: $f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}$

이제 $\dfrac{\pi}{2}$에서도 비슷해졌네요! 다음!

▶ 7차 근사식: $f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!}$

$pi$ 근처까지도 꽤 비슷해졌습니다. 몇 개만 더 해 봅시다.

▶ 9차 근사식: $f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}-\dfrac{x^7}{7!} + \dfrac{x^9}{9!}$

거의 원본 함수와 같아졌습니다. 사실 여기부터는 이제 $0 < x < \pi$일 때는 그냥 가져다 써도 오차는 크지 않을 것 같습니다. 어차피 $sin$은 주기함수이고 대칭함수이니까 $0 < x < \dfrac{\pi}{2}$에서만 정확해도 다른 범위에서는 함수의 특징을 이용해 계산해내면 됩니다.

▶ 11차 근사식: $f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!}$

차이가 얼마나 나는지 점점 이렇게 봐서는 확인하기가 어렵습니다. 마지막으로 13차 근사식까지만 그려보겠습니다.

▶ 13차 근사식: $f\left ( x \right ) = x- \dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!} + \dfrac{x^{13}}{13!}$

그래프 범위 내에서는 거의 똑같은 곡선이 그려집니다.

프로그래밍 언어에서의 구현 사례

OpenJDK의 StrictMath 네이티브 코드가 13차 근사식을 통해 삼각함수를 구현하고 있습니다. 정확성을 위해서 13차 근사식을 그대로 사용하진 않고, 이렇게 조금 변형해서 사용합니다.
\[\sin x \sim x- \dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!} + \dfrac{x^{13}}{13!}\]
\[\frac{\sin x}{x} \sim 1- \dfrac{x^2}{3!} + \dfrac{x^4}{5!}- \dfrac{x^6}{7!} + \dfrac{x^8}{9!}- \dfrac{x^{10}}{11!} + \dfrac{x^{12}}{13!}\]
그러므로 $r$을 이렇게 정의할 때
\[r = x^3 \left(\dfrac{1}{5!} + x^2 \left(-\dfrac{1}{7!} + x^2 \left(\dfrac{1}{9!} + x^2 \left(-\dfrac{1}{11!} + {x^{2}} \cdot \dfrac{1}{13!}\right)\right)\right)\right)\]
$\sin x$는 이렇게 표현할 수 있습니다.
\[\sin x \sim x + x^3 \cdot \left(-\dfrac{1}{3!} + x^2 r \right)\]
예를 봅시다. sin 함수에서 __kernel_sin 함수를 호출할 때 코드에서 $|x| \prec \pi/4$라면 y, iy의 값은 모두 0입니다. 이외의 범위에서는 적절한 범위로 평행이동시켜 __kernel_sin 혹은 __kernel_cos 함수값을 구합니다.

#include "fdlibm.h"
#ifdef __STDC__
    static const double
#else
    static double
#endif
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */

#ifdef __STDC__
    double __kernel_sin(double x, double y, int iy)
#else
    double __kernel_sin(x, y, iy) double x,y; int iy; /* iy=0 if y is zero */
#endif
{
    double z,r,v; int ix;
    ix = __HI(x)&0x7fffffff; /* high word of x */
    if(ix < 0x3e400000) /* |x| < 2**-27 */ {
        if((int)x==0) return x;
    } /* generate inexact */
    z = x*x;
    v = z*x;
    r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
    if(iy==0) return x+v*(S1+z*r);
    else return x-((z*(half*y-v*r)-y)-v*S1);
}

소스에서
$S1 = -\dfrac{1}{3!} \approx -1.66666667 \times {10}^{-1}$
$S2 = \dfrac{1}{5!} \approx 8.33333333 \times {10}^{-2}$
$S3 = -\dfrac{1}{7!} \approx -1.98412698 \times {10}^{-4}$
$S4 = \dfrac{1}{9!} \approx 2.75573137 \times {10}^{-6}$
$S5 = -\dfrac{1}{11!} \approx -2.50521084 \times {10}^{-8}$
$S6 = \dfrac{1}{13!} \approx 1.58969099 \times {10}^{-10}$
으로 근사식의 계수들이 근사되어 하드코딩되어 있는 것을 알 수 있습니다. 또한 $z=x^2$, $v=zx=x^3$으로 계산하고 있습니다. 13차 근사식을 그대로 쓴다면 항마다 1번, 3번, 5번, …, 13번의 곱셈을 해야 하지만 $x^2$를 새 상수로 두면 곱셈을 줄일 수 있어서 이런 방법을 사용하지 않았을까 하는 생각입니다.

다른 방법으로도 계산할 수 있나요?

물론 다른 방법들도 있습니다. 예를 들어 볼더가 1956년에 고안한 CORDIC(COordinate Rotation DIgital Computer) 알고리즘을 이용해 삼각함수의 값을 근사하는 것도 가능합니다.

간단하게 설명하자면, 맨 위에서 회전을 하기 위해 곱하는 행렬을
\[\begin{bmatrix}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta\end{bmatrix}\]
으로 소개했습니다만, CORDIC 알고리즘은 이를 $\tan$ 함수만으로 표현해
\[R_i = \dfrac{1}{\sqrt{1 + \tan^2 \left(\gamma_i\right)}}\begin{bmatrix}1 & -\tan \left(\gamma_i\right) \\\tan \left(\gamma_i\right) & 1\end{bmatrix}\]
로 변형하고, 컴퓨터가 빠르게 계산할 수 있도록 $\tan \left(\gamma_i\right) = \pm2^{-i}$가 되는 $\gamma_i$, 즉 $\arctan\left(2^{-i}\right)$의 값들과 이 때의 $\dfrac{1}{\sqrt{1 + \tan^2 \left(\gamma_i\right)}}$의 값들을 하드코딩해 두고 원하는 각도에 가까워질 때까지 $i$를 하나씩 증가시키면서
$\arctan\left(2^{-i}\right)$를 더하거나 빼면서 행렬 계산을 해나가는 알고리즘입니다.

이 알고리즘은 계산기에 주로 쓰이고 있고, 인텔의 8087 ~ 80486 CPU에도 채용되어 왔습니다.

컴퓨터는 이렇게 삼각함수를 계산합니다.