関係者からのリクエストで書いた記事です．おそらくほとんどの方は「誰得……」と思うことでしょう．
なお，推測で書いている部分もあるので，以下には間違いが含まれている可能性があります．記述に間違いを見つけた場合はご一報ください．メールアドレスは home にあります．
「Linux の cosf は精度が悪い．Windows の MinGW の方が良い結果が出る」と言われたので，その真相を調べた結果，どうやらこれは glibc の開発チームの意図通りの結果で，確かにそれが問題になる人もいるかもしれないが，まあ個人的にはこれぐらい許容範囲かなあ，という結論が出ました．
2018年12月の頭のこと，某研究会での数値計算についての講演で「最新の Linux (Ubuntu) では cosf の結果と cos を float にキャストした結果が合わないことがある．この講演の提案法では cosf を多用するため，この差が蓄積して大きな精度悪化を引き起こす．Windows の MinGW では両者の値は一致して精度悪化は起こらない．」という報告がありました．普段から Gentoo Linux を使っている自分としては「本当かなあ．もしそんな問題があるのならばバグ報告するべきなのでは」ということで調べ始めました．
ちなみに，cos ではなくわざわざ cosf を使う理由は（私の理解では）倍精度演算が単精度演算よりはるかに遅い GPU との比較実験を将来的にはしたいから，という数値計算屋の都合以外の何物でもない話で，これをもって「Linux は信用できない」と言われるのは難癖もいいところなのですが，あまり気にしないことにします．ただし，あまりにも誤差が大きいようだと問題なわけで，標準的な数学ライブラリの誤差について知っておくことは意義があると言えます．
主張を確かめるために，次の簡単なコードを走らせました．
#include <stdio.h>
#include <math.h>
int main(void){
float d;
for(d = 1e-7; d < 4; d = nextafterf(d, 4))
if((float)cos(d) != cosf(d))
printf("NG: %e, err = %e\n", d, (float)cos(d)-cosf(d));
return 0;
}
要は 1e-7 から 4 まで，nextafterf を使って 1 ulp ずつ加え，この範囲で float で表現可能な全ての値について cosf と cos を float にキャストした値を比較して，一致しない場合はその差を出力します．実験に用いた PC は次の通りです：
% uname -a Linux localhost 4.19.2-ck #1 SMP Sat Nov 24 11:15:45 JST 2018 x86_64 Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz GenuineIntel GNU/Linux % gcc --version gcc (Gentoo 8.2.0-r5 p1.6) 8.2.0 ...
調べたところ，後述するように glibc の cosf 関数は最近頻繁に更新されているということがわかりました．そこで，2.25, 2.28, devel (執筆段階で最新の 7c9a7c) の3つのバージョンを比較することにしました．2.28 は現在リリースされている最新版で Gentoo の Portage でビルドしたもの，2.25 と devel はソースをとってきて（--prefix 以外は）特に何もオプションをつけずに ./configure; make -j8 したものに /lib64/libm.so.6 からシンボリックリンクを張ることで実験しました．なお，2.28 は念のため同様にソースからビルドしたものも試して，結果が Gentoo の色々パッチをあてたものと同じであることを確認しています．この他のバージョンについて，調べた限りでは，2.27 は 2.28 と同じ，2.24，2.26 は 2.25 と同じ結果になります．
結果のファイルは大きいので圧縮して置きました (result_err_cosf_glibc.zip)．出力の行数（一致しなかった引数の個数）を wc -l で出した結果は次の通りでした：
なお，どの引数でも誤差は 1 ulp のようです．過去のバージョンに比べて，最近の cosf は 1 ulp とはいえ確かに精度が悪化していることがわかります．2.28 から devel では改善しているように見えますが，よく見ると 2.28 では 3.125000e-02 より小さい引数に集中しているのに対し，devel では誤差の出る引数が 4 までに満遍なく広がっていることがわかります．分布の密度は下がっているとはいえ，使用法によっては精度がより悪化する原因になる可能性があると言えそうです．
これに加えて，Windows 機で Cygwin, MinGW, Visual Studio 2017 を使った同様の実験もしてみたところ，Cygwin は（出力はファイル化しておらず，ちゃんと見てもいませんが）最近の glibc とほぼ同じ結果，MinGW と Visual Studio 2017 は一切出力が出ない（！）という結果が得られました．よく知られている話のようですが，Cygwin の C ライブラリは cygwin1.dll に含まれる独自のもので，MinGW は単に msvcrt.dll にリンクするようです．このことは ldd コマンドで簡単に確認できます．したがって，MinGW と Visual Studio では常に同じ結果になるといえます．冒頭の講演をした人にこの結果を伝えたところ，「おそらく msvcrt の cosf は単に cos の結果を float にキャストしているだけなのではないか」ということでした（ソースがないので確かめようもありませんが）．
単精度での 1 ulp ぐらい別にええやん……，とも思いますが，まじめに調査して考えてみます．
これも冒頭の講演者経由で伝え聞いた話ですが，元々C言語の単精度演算というのは倍精度にキャストしてから演算をして，結果を単精度に丸めるという実装が主流だったそうです．実際，古い資料と思われますが「Sun WorkShop 6.0 C プログラミング言語コンパイラ」のマニュアルの付録「K&R Sun C と Sun ANSI/ISO C の違い」には次の記述があります：
- 単精度計算と倍精度計算
- Sun C (K&R): 浮動小数点式のオペランドを double に拡張する。float を返すように宣言された関数の戻り値は、常に double に拡張される。
SUN ANSI/ISO C: float の演算を単精度計算で行うことができる。このような関数に float の戻り型を使用できる。
この古い時代の慣習を引き継いでいれば，msvcrt の cosf の結果が cos を float にキャストした結果と一致するのも自然であるといえます．
一方，技術の進歩は著しく，近年の CPU では SIMD 演算が標準装備されていて，特に最近の Intel の製品には AVX2 や FMA といった最新技術が投入されています．どうも最近の glibc の cosf の精度が落ちたのは，この積和演算を利用した高速化を狙ったコードの変更のためのようで，実際に精度が落ちた 2.27 では FMA による高速化が導入されたという記事が見つかります．ちなみに，この 1 ulp の誤差というのは，glibc の方針によるとバグではないようです．詳しくは次を読んでください．
下のリンクは実装ミスで cosf(1.57079697) の結果が 3 ulp ずれているというバグ報告ですが，修正のコミットメッセージで次のように明言されています：
(To be clear, this is a quality-of-implementation issue relating to the apparent intent of those particular cosf and tanf implementations; 3ulp is within the general glibc accuracy goals, so not inherently a bug.)
3 ulp の誤差というのは，glibc の開発方針としてはバグではないのだが，意図通りの実装になっていないので修正したということです．さらに，2.28 以降の commit ea6c662 には次のようにはっきりと書かれています：
This patch is a complete rewrite of sincosf. The new version is significantly faster, as well as simple and accurate. The worst-case ULP is 0.5607, maximum relative error is 0.5303 * 2^-23 over all 4 billion inputs. In non-nearest rounding modes the error is 1ULP.
真値からの 0.5607 ulp までの誤差を許容するので，丸めると cos をキャストした結果とは 1 ulp ずれることがありえるというわけです．この記述から，今回の cosf の結果は計算速度の改善を狙ったために生じた意図的なものだということがわかります．実は 2.27, 2.28 の cosf の実装は，ソース ([glibc.git]/sysdeps/ieee754/flt-32/s_cosf.c) を読むと 0x1p-5 = 0.03125 未満の場合はチェビシェフ近似式の次数を落としてサボるようになっていることがわかるので，これによる誤差が意図していないものである可能性があるならば報告しようかとも考えていたのですが，commit ea6c662 を見て devel の誤差も十分にテストをした意図的な結果だということがわかったので，報告はやめました．
まあ，私にとっては単精度の 1 ulp とかどうでもいいんですがね（ひどい）．単精度で計算する時点で……，というわけですが，1 ulp ずれるのがどうしても困る場合は cos で計算してキャストしましょう，というのが結論です．ネタとしては面白かったので，このような機会を与えていただいた冒頭の講演者の方に深く感謝します．
最後に，あまり「どうでもいい」とか書いてしまうと怒る人がいるかもしれないので補足します．今回調べたのは float ですが，おそらく double でも glibc の数学関数の結果は 1 ulp 程度ずれることは十分にありえると思われます．したがって，特に数値計算など精度が必要な状況で数学関数を利用する場合は，その精度で一番下の桁が1ずれることもありえると覚悟して，それで困る場合はより精度の高い関数を自分で組むか，4倍精度や多倍長を利用することを検討する必要があるでしょう（ただし，自分で数学関数を組むのは信頼性の面でハードルが高いと思われます）．
追記：柏木雅英さんが過去に glibc の double の数学関数の精度を調べた記事を見つけました．Ubuntu 14.04 なので，glibc 2.18 か 2.19 でしょうか．「1 ulp より大きい誤差を検出」という基準だと，sin, cos などは現行版でもパスでしょう．0.5 ulp より大きい誤差は十分にありえることに注意する必要があります．
従って、C言語などから使うexp, sinなどの数学関数の返す値は信用することは出来ません。kvライブラリの区間演算でも数学関数は全て自力で実装しています。
ということで，ガチ勢はやっぱり自分で実装するんだなあ，と納得しました．