🐡

ちょっとだけ怖い三角関数のはなし

2020/11/04に公開

math

tech

TL;DR

sin()とかcos()の計算時、つい計算誤差を「引数に」入れがち

よくある課題

競技プログラミングとかでよくある

x,y平面上の単位円をn分割する点の座標をdouble精度(16桁)で表示する関数を作りなさい。ただし一つ目の点の座標は(1,0)とする

みたいなやつをなーんも考えずに、たとえばC(C99)で書くと以下のような感じになるかと思います

#include <stdio.h>
#include <math.h>

void pointsOnCircle(int n)
{
    double pi = 3.14159265358979323;
    for (int i = 0; i < n; i++) {
        printf("(%1.16f, %1.16f)\n", cos(2 * pi / n * i), sin(2 * pi / n * i));
    }
}

じゃあこれ実行するとうなるでしょうか。int main() { pointsOnCircle(2); return 0;} みたいな関数とくっつけてコンパイル・実行すると

(1.0000000000000000, 0.0000000000000000)
(-1.0000000000000000, 0.0000000000000001)

なんということでしょう。二つ目の頂点のy座標が予想外のところに来てしまいました!
n=4とかでも似た現象を確認できます。

もちろんこうなるのはコンパイラのバグとかじゃないです。規格に準拠した、正しい動作をするコンパイラならこうなるのが普通です。また、詳細は省きますが、floatやlong doubleとかを使ったとしても同様の問題が発生します。

なんでそうなるのか

C言語(に限らず、多くの計算機上の言語)のsin(), cos()といった関数は、数学的な同名の関数を、できるだけ近似したモノになります。なのですが、今回の誤差についてはこの近似が原因ではありません。

原因は sin() に与えている引数にあります。

最初に上げた例にあった、n=2 における二つ目の頂点のy座標はsin(2 * pi / 2), つまりsin(pi)です。piには適当に長い桁数まで $\pi$ の値を入れていますが、本物の $\pi$ ではなく、適当な桁で丸められた近似値になります。この近似が問題の原因となっています。この近似は一般的に浮動小数点値を使う場合に常に発生してしまうため、多くの計算機言語において、 $\sin(\pi)$ を計算しようとしてストレートに引数に $\pi$ ぽい値を使ったのでは、どう頑張っても期待する値(正確に0)を得ることはできません。

数学的に $\sin(x-\pi) = -\sin(x)$ なので、C言語において sin(pi) を計算すると、それはつまり $\sin(3.1415926535....) = \sin(3.1415926535... - \pi)$ になります。精度の高い変数を使って、sin(x)のxを限りなく $\pi$ に近い値にしていっても、ご存じの通り $\pi$ は無限の桁数がある値なので、x - π はとても小さい値にはなりますが、0になりません。

そして、これまた数学的に $\sin(x) = x - x^3/6 + x^5/120 - ...$ なので、 $x$ がとても小さい値の時、 $\sin(x) \simeq x$ です。x をどんなに限りなく $\pi$ に近づけていっても、両者の誤差と同じぐらいの誤差が $\sin(\pi)=0$ と sin(x) の間に発生することが分かります。

じゃあどうすんだよ

この例では競技プログラミングっぽい課題を出しましたが、FFTとか、円や球面をポリゴンで分割近似するとかいった実アプリケーションとかでもごく普通に使うコードです。場合によってはこの計算誤差が無視できないことがあります(筆者は実際にこれによっていくつかのバグにはまりました…)。

解決するには、数学的に頑張るのが良いと思われます。つまり、最初のコード例において、 sin(2 * pi / n * i) を計算する際に、i と n/4, n*2/4, n*3/4 に対する大小関係を使って、sin(2*pi/n*i), cos(2*pi/n*(i-n/4.0)), -sin(2*pi/n*(i-n/2.0)), -cos(2*pi/n*(i-(3*n)/4.0)) を使い分ける、というものです。

この方法で最初に例に上げた問題は解決します。ただ、詳細は略しますが、さらに精度にこだわるのであれば、(i-n/4.0), (i-n/2.0),... のような部分の「絶対値」が最小になるような選び方をしたほうが、より良い結果が得られることがあります。

別の解として、より高い精度の変数を使って計算するというお手軽方法もありますが、これだと今回問題にした部分以外でいわゆる「二重丸め」の問題が起きてしまうほか、高精度浮動小数点フォーマットに対応してない環境では解決しないという根本的な問題が起きてしまいます。

まとめ

こんなんで結果が変わるようなアルゴリズム、めんどくせーよ！あーもー

と言いたくなるわけですが、たかがsin, cosを求めるだけでこんな罠がある、ということで。