サイン、コサインをインテルの CPU で計算すると少しバグっているらしい
こんにちは。
ちょっと前にどこかの偉い人が「女性にサイン、コサイン教えても仕方ない」ような発言をして炎上したのは記憶に新しいところです。面白いですね。
「女の子にサインコサイン教えて何になる」 鹿児島県知事の発言がネットで非難 発言撤回の会見動画も公開 – ねとらぼ
http://nlab.itmedia.co.jp/nl/articles/1508/28/news145.html
「いやいやサイン、コサイン使うでしょ」という反論には完全に同意です。それはよいのですが、プログラマー的な視点で見れば「サイン、コサインはコンピュータでどうやって求めているか」を説明できてやっと安心（？）ではないでしょうか。
そこでちょっと復習と思って調べてみたら、なんかインテルのCPUのバグ？の話とか出てきてわけが分からなくなってきたので、まとめておきました。
sin(x) の求め方（基本）
理系の方は学校で習ったと思いますが、テーラー展開した多項式で求めるのが常套手段です。
例えばサインの値は、以下の多項式で表現できます。
自然科学のための数学I/II（2014年度）- 琉球大学
http://www.phys.u-ryukyu.ac.jp/~maeno/sizensuugaku/lec11_sin.html
これを実際に計算するサンプルプログラムとかはネット上にごろごろしてるので、ここでは挑戦しません。ここでは、この計算が最近のコンピュータのどの部分で行われるかに注目してみます。
コンピュータのどの部分で計算しているのか
コンピュータで sin(x) の値を計算させるときは、昔も今も上記の多項式が使われているようです。（※他のアルゴリズムもありますがあまり使われていないらしい）
最近（と言っても 8087 にすでにあったようなので20年以上前から）はハードウエアレベルでその機能が提供されています。別の言い方をすると、sin(x) の値を求めるアセンブラ命令があるということです。インテルの x87 では FSIN, FCOS などです。
なので当然 sin(x) の実装にはこのアセンブラ命令 FSIN が使われていると思いますよね？
しかし念のため確認してみると意外な事実が判明しました。確認したのは GNU の C 標準ライブラリ glibc のソースで、バージョンは 2.22.90 です。
GNU C Library (glibc)
おそらく sin(x), cos(x), tan(x) とかはそのハードウエアのアセンブラ命令を叩くだけのはず。
sin(x) とかはアーキテクチャ依存の実装になるので、double __sin(double) の実装があちこちにあって悩ましいですが、とりあえず x86_64 アーキテクチャでは、
- sysdeps/x86_64/fpu/multiarch/s_sin.c
- sysdeps/ieee754/dbl-64/s_sin.c
の実装が有効のようです。すると、
<sysdeps/ieee754/dbl-64/s_sin.c>
https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_sin.c;hb=HEAD#l279
あ、あれ？
なんか思っていたのと違いますね・・・と言うか、完全にソフトウェア的に計算しています。
なんでハードウエアサポート使わないの？
せっかくハードウェアサポートがあるのに何で使わないんだろうか・・・。って思っていろいろ調べてみましたが、以下ページを見つけました。
Intel overstates FPU accuracy – 06/01/2013
http://notabs.org/fpuaccuracy/index.htm
リンク先を見てみると「インテルの FPU (x87) のコサイン、タンジェント命令は、スペック通りの精度がでない場合がある。」と言った趣旨のことが書かれています。
また以下のページにも興味深いことが書いてあります。
How does C compute sin() and other math functions? – Stackoverflow
http://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions
この質問に対する回答によると、「2011年の10月以降 x86-64の Linuxではこの（ソフトウェア）コードが使われており、明らかにFSINより速い（Since October 2011, this is the code that actually runs when you call sin() on a typical x86-64 Linux system. It is apparently faster than the fsin assembly instruction.）」とのことです。
というわけで、現在かなり普及している x86_64 アーキテクチャー（Windowsも Mac も大部分が採用している）の CPU の三角関数命令 FSIN は、速度、精度などが問題視されており（少なくとも glibc では）、別途ソフトウェア的に実現されているようです。
他の実装ではどうなっているかは不明ですが、glibc での実装と似たようなことをやっているのではないかと思います。
インテルも大変なんだな・・・。
まとめ
- サインなどの三角関数の値はテーラー展開を使用した有名な式を昔も今も使って計算していることが多い。
- 20年以上昔から浮動小数点数演算用のFPU(コプロセッサー) には三角関数の値を計算する命令（ハードウェアサポート）がある。
- しかしインテルの x86_64 系アーキテクチャーのこの命令は精度及び速度に問題があるらしく、少なくとも glibc 2.22 ではFPUのサイン命令は使っておらず、ソフトウェア的に計算しているらしい。
という感じでした。
何はともあれ glibc の sin 関数のソースは勉強になるので、数値計算に興味のある方には一度じっくり見てみることをおすすめします。
Pingback: （インテルのCPUの）三角関数は不要 | 配電盤