🐷

GPUの単精度 sin / cos 関数の誤差を引数 32bit 全探索で調べた

に公開

GPUの単精度 sin / cos には「標準版」と「近似版」がある

NVIDIA / AMD など多くの GPU には、sincos といった超越関数を高速(1サイクルあたり1命令のスループット[1])に 近似計算する 専用ハードウェアが搭載されています。NVIDIA ではこれを SFU (Special Function Unit) [2]と呼びます。本記事では便宜上、SFU で計算されるものを 近似版 と呼び、非SFUのほうを 標準版 と呼びます。

Unity をはじめとするゲームエンジンで何気なく sin / cos を Shader(HLSL/GLSL)から呼ぶと、HLSL → 中間表現 → GPU アセンブリ となってSFU 命令に置き換えられ ます[3] [4]
このときの誤差特性を把握していないと、回転行列FFT 系など三角関数を多用するアルゴリズムで想定外のブレが出ることがあります。

本稿ではNVIDIA/AMD/Intel GPUで SFU 版の sin / cos 命令がどの程度ズレるのか を明らかにします。せっかく GPU を使うので引数32bit(≈ 42 億通り)を全探索する方針で臨みました。

まず仕様から

NVIDIA SFU 版 sin / cos の保証範囲は -100pi ~ +100pi

NVIDIA GPU の SFUが提供する三角関数命令は、入力範囲に応じて誤差の上限が定義されています。
Parallel Thread Execution ISA Version 9.0によると、

  • 入力範囲が -2pi ~ +2pi の場合
    最大誤差は 2^(-20.5) 程度
  • 入力範囲が -100pi ~ +100pi の場合
    最大誤差は 2^(-14.7) 程度

と明示されています。これを超える範囲の入力については 「best effort」 とだけ記載されており、誤差に関する保証はありません。
そんなことを知らず、自作Shaderで sin(time * x) のように書きまくってましたけどね・・
入力が大きくなると精度も劣化するのですが、その詳細は以下の検証で。

※誤差範囲や数値はアーキテクチャやベンダで異なります。該当 GPU/ドライバのリファレンスで確認を、と言いたいところですが実質NVIDIA以外公開してないので、https://zenn.dev/toropippi/articles/d8def1c994e0a9#検証4.1:nvidia保証範囲内(±100pi)でamd%2Cintelのsfu版絶対誤差最大値 このセクションを参考にしてください。

地味に重要ですが、"absolute error" = 絶対誤差 のみが保証されている点に注意です。後に扱いますがULPは相対誤差のほうを見ています。


標準版 sin / cos は「最大 2 ULP」まで(単精度)

NVIDIA の CUDA C++ Programming Guide では、単精度の標準版 sinf / cosf(※SFU を使わずソフトウェア近似で計算)は最大誤差 2 ULPとされています。速度は遅い代わりに精度は高いという位置づけです。

この「標準版」は OpenCL / CUDA などの GPGPU API で普通に使える一方、Unity の Shader(Compute Shader 含む)、つまり HLSL / GLSL からは原則として直接選べません(対応する命令や切り替え手段がないため)。少なくとも私の調べた限りそのようです。
どのような書き方が 近似版/標準版 関数になるか、という点についてはもう1本記事をかく予定です。

ULP とは?

ULP(Units in the Last Place)とは、「浮動小数点数で隣り合う表現可能な値の差」を 1 と数える単位です。
たとえば、0.25f (0x3E800000) と 0.2499999850988388f (0x3E7FFFFF) の差はちょうど 1 ULP です。


標準版より高精度な「正解」をどう作るか

誤差検証には、「1 ULP のズレも無視できない基準値」が必要です。単精度の標準版でも最大 2 ULPの誤差があるので基準としては不十分。そこで本稿では、

  • doublesin / cos を計算してから float に丸め(キャスト)

で得られた値を正解とました。
(この時点で最大0.5 ULPの誤差がありますが他にいい方法が思いつきません)

方法

検証に使った GPU 等

  • RTX 5090 (Windows)
  • RTX 3080 (Windows)
  • RTX 2080 Ti (Windows)
  • RTX 2080 Ti (Ubuntu)
  • Titan V (Ubuntu)
  • Ryzen 9 9950X 内蔵 GPU gfx1036 (Windows)
  • Intel UHD 770 (Windows)

使用コードと対象関数

以下の 4 関数(OpenCL)を各 GPU で評価しました。

  • sin:単精度・標準版
  • cos:単精度・標準版
  • native_sin:単精度・SFU 近似版
  • native_cos:単精度・SFU 近似版

基準値として

  • (double)sin → float へ丸め
  • (double)cos → float へ丸め

も取得しました。キャストは最近接丸め(round-to-nearest-even)をしました。

引数は全探索なので、0x00000000 から 0xffffffff までのすべてを対象にしています。この範囲を浮動小数点数(float)として解釈すると、NaN や Inf、+-0.0、非正規化数も含まれますが、それらも計算対象にしています。関数の出力結果はファイルに分割して保存しました。

ちなみに 32bit 引数全パターンを 4 関数 × GPU 台数でやると未圧縮で481 GBになります!


検証結果と考察

検証1:同一ベンダ(NVIDIA)内で差は出るか? → 出ない

  • RTX 5090 / 3080 / 2080 Ti(Win/Ubuntu)/ Titan V(Ubuntu)の 5 GPU で関数出力内容に違いがあるか比較
  • sin(標準)/cos(標準)/native_sin(SFU)/native_cos(SFU)で検証

すべて bit 単位で完全一致でした(Inf / NaN / 非正規を含む、誤差の出方まで一致)。OS の違いも無関係。
以降、NVIDIA 代表値は RTX 5090で示します。


検証2:SFU近似版の誤差に傾向はあるか? → 引数が0.0, pi/2, pi近傍で精度悪い

いくつかの入力で誤差を観察

まずは点検用に手動でランダム数を引数とし関数誤差を観察します。

x = 1.23

[sin] native: 1 ULP,  [cos] native: 1 ULP   (RTX 5090)

1 ULP 程度に収まっており非常に優秀。

x = 0.8756

[sin] native: 2 ULP,  [cos] native: 2 ULP   (RTX 5090)

— こちらも 2 ULP に収まっており優秀。

x = 0.0123(0 近傍)

[sin] native: 87 ULP, [cos] native: 1 ULP   (RTX 5090)

x = 0.000123(より 0 近傍)

[sin] native: 10,760 ULP, [cos] native: 1 ULP   (RTX 5090)

x = 1.571(≈ pi/2 近傍)

[sin] native: 1 ULP, [cos] native: 6,895 ULP   (RTX 5090)

このように手動入力の結果、0,pi近傍ではnative_sinのULPが上昇、pi/2近傍でnative_cosが上昇している傾向をつかめました。
関数の出力が0.0に近いほうでULPが上昇する ということのようです!

ベンダ横断の比較(抜粋)

引数 NVIDIA GPU AMD (gfx1036) Intel UHD 770
-2.34 4 ULP 1 ULP 0 ULP
0.0123 87 ULP 4 ULP 93 ULP
0.000123 10,760 ULP 4 ULP 6,418 ULP
1.571 6,895 ULP 7,458 ULP 2,800 ULP
3.141598 1,019,989 ULP 196,881 ULP 70,079 ULP

※ ULPはnative_sinとnative_cosで大きい方を採用

  • AMDのnative_sinは例外的に 0 近傍に強い傾向だがあとはNVIDIA同様
  • Intel は NVIDIA, AMDよりやや善戦

検証コード
https://github.com/toropippi/GPU_SinCosULP/blob/main/2_結果表示/peek_trig_result.py


有効桁数が1桁もないことも!?

3.141598 での誤差が非常に大きかったので、良い例(?)としてここに紹介します。

まず入力引数 3.141598 はfloat表現として 3.1415979862213135 として解釈されました。REF.sin がdouble型を使った正解値で TEST.sin が標準版sinの結果です。

  REF.sin:     -5.332631644705543e-06 bits=0xb6b2eef5  (64_sin)
  TEST.sin:    -5.332631644705543e-06 bits=0xb6b2eef5  (64_sin)

どちらもbit表現が一致しており0 ULPであることがわかります。しかし・・

  TEST.native: -4.86879434902221e-06 bits=0xb6a35ea0  (ULP=1019989) 

SFU版sin では正解値 -5.333e-06 に対し -4.869e-06 とかなりずれていることがわかります!って、有効桁数1桁もないじゃないですか!?
言ってしまえばこれは 意地悪ケース なのですが、ひどい例ではこういうこともあると認識しておく必要があります。


検証2(追加):近似版sin(),cos()の境界付近のみ可視化

0.0, pi/2, pi近傍の引数で悪くなるようだったので、これら境界付近の生の値をみてみましょう。native_sin0.0pi付近、native_cospi/2付近の値をグラフにしました。

  • native_sinの0.0~0.000002

プロットは点のみです。NVIDIAとIntelが一定区間内で同値を返していることがわかります!いいのか?と思いますが、近似版関数は精度を犠牲にして速度を最優先にしているので仕方ありません。一方0.0近傍に強いAMDはちゃんと正解に近い値を返しているようです。

  • native_cosのpi/2 ± 0.000002

さっきと比べ点の数が少ないですが、 1.57 ± 0.000002 の極めて狭い範囲を見ているためです。32bit浮動小数点数では隣り合う点と点のあいだに表現可能な値は存在しないので、線で補間しプロットしています。
Intelがいい感じの線を描いている一方、AMDとNVIDIAはガタガタしているのがわかります。

  • native_sinのpi ± 0.000002

同じようにIntelはいい線いってますがNVIDIA、AMDでガタガタしています。

ULP 高値の原因

なぜ関数の出力が 0.0 に近いときに ULP が増大して見えるのかというと、近似版関数では仕様で 絶対誤差のみが保証 されているためです。ちょうどさっきの図が参考になるかと思います。絶対誤差が常に一定とすると、出力される値が小さければ小さいほど 相対誤差 が大きくなりますよね。
ULP は実質的に相対誤差を測る指標であり、浮動小数点数のビット表現に基づいて評価されます。表現値そのものが小さいと、同じ絶対誤差でも仮数部に対して大きなビットの変化として現れ ULP 値が大きく見えるのです。

また、SFUのnative_sinとnative_cosはほぼ同一回路で計算されていると思われます。

SFU 内部では 小区間分割+位相変換 のような実装が想像され[4:1]native_sin(x+pi/2) のULPはそのまま native_cos(x) のULPを見ていることになっていると思われます。これにより(pi/2)(0,pi) の引数でnative_sinとnative_cosのULPがきれいに逆転していることの説明がつきます。

余談:ゲーム開発視点では相対誤差が重視される
Shader で sincos が最も多用されるのは、回転行列やクォータニオンの生成などです。
これらはベクトルの長さや直交性が保たれていることが前提であり、もし cosθsinθ に相対誤差があると、回転を何度も重ねるうちに少しずつ行列の直交性が崩れ、オブジェクトがじわじわスケーリングしたり、姿勢がわずかに傾くといった問題が起こります。

Unity(Shader Model 5.0 想定)では float 型を使えば最大精度の演算になりますが、「より高精度な三角関数に切り替わる修飾子がある」といった公式情報は見つかっていません(AI で調査済み)。

そのため、0°・90°・180° などに非常に近い角度 で回転を繰り返すケース等では実際の描画結果を目視で確認する必要がありそうです。


あなたの環境ではどうですか?

意地悪なコーナーケースも含めた9個の引数でGPU sin()命令の結果を確認できるサイトを作成しました。
WebGL検証サイト(https://toropippi.github.io/webgl-sin-sfu/)

スマホでも確認できます!


検証3:SFU近似版の結果 ULP可視化

[-pi, pi] の平均ULP棒グラフ

続いて引数xを-pi~pi区間に限定して傾向を探ります。区間を400分割し、全探索で得られた全ULPデータを平均しプロットしました。

  • NVIDIA RTX 5090

区間平均のULPは

関数 平均ULP (-pi~pi)
native_sin 96.176 ULP
native_cos 209.38 ULP

でした。
一応区間平均のULPを出してみましたが境界付近が大きく平均を上げているようです。それを除けば近似版三角関数でもほとんどの場合非常に優秀な精度がでている、という認識で良さそうです。
sinのほうは0.0, pi 近傍で、cosのほうはpi/2 近傍でULPが上昇、そして8,9割がたの区間で0~4 ULP程度におさまっていることがわかります。

補足:重み付け平均の意味
浮動小数点数は値の大小によって刻み幅(隣り合う数の間隔)が変わるため、単純平均では小さい値のサンプルが過剰に反映されてしまいます。そこで各サンプルを刻み幅に応じて重み付けしたのが重み付け平均の意味です。x が一様分布しているときの平均ULP を求めています。
本稿では以降、重み付け平均を"平均"と省略しているところがあります。


  • AMD gfx1036

sinのほうはNVIDIAと違い0.0 近傍で優秀な精度となっておりますが pi 近傍ではULPが上昇しています。cosのほうはNVIDIA同様 pi/2 近傍でULPが上昇。
区間平均のULPは

関数 平均ULP (-pi~pi)
native_sin 29.625 ULP
native_cos 138.81 ULP

でした。


  • Intel UHD 770

NVIDIA, AMDと違った結果になりました。sin, cosとも 0.0, pi/2, pi 近傍でULPが上昇。そのかわり半分以上の区間で0~1 ULPと非常に高い精度を叩き出しています。

区間平均のULPは

関数 平均ULP (-pi~pi)
native_sin 43.441 ULP
native_cos 44.250 ULP

でした。

一旦ここまでまとめると

  • 0.0, pi/2, pi 近傍 でsin,cosどちらかor両方のULPが上昇する傾向は3社ともほぼ共通
  • それ以外の区間では0~10 ULPと誤差が非常に小さい → ゲーム用途で十分
  • AMDでは例外的に0.0近傍のsinもULPがほぼ上昇しない(<10 ULP)
  • Intelは実装が他2社と大きく異なるのかsin(pi/2)やcos(0.0)やcos(pi)でもULP上昇
  • そのかわりIntelは全体的に精度が優秀。平均ULPが低め
  • NVIDIAとAMDでは若干AMDが高精度な傾向

[-10pi, 10pi] の平均ULP棒グラフ

周期的にULPが跳ねる傾向が続くか確かめるためにNVIDIAでのみ[-10pi, 10pi]の範囲を調査しました。

  • NVIDIA RTX 5090

周期的に急峻にULPが上昇していることが確認できました。cosも同様だったので割愛


[-100pi, 100pi] の平均ULP棒グラフ

  • NVIDIA RTX 5090

    cosのグラフは同様の傾向だったので割愛してます。

区間平均のULPは

関数 平均ULP (-100pi~100pi)
native_sin 17166.6 ULP
native_cos 10139.3 ULP

でした。

  • AMD gfx1036

区間平均のULPは

関数 平均ULP (-100pi~100pi)
native_sin 10223.6 ULP
native_cos 10123.3 ULP

でした。

  • Intel UHD 770

区間平均のULPは

関数 平均ULP (-100pi~100pi)
native_sin 6924.81 ULP
native_cos 182.745 ULP

でした。

ここまでまとめると

  • NVIDIA, AMDはpi周期で,Intelはpi/2周期でULPが急峻になる傾向が続いている
  • NVIDIA, AMDでは引数が大きくなるとその 1乗に比例 して平均的なULPも上昇してそう?
  • Intelも周期的にULPが急峻に上昇する傾向が続いているが、引数の増大の影響はほぼなし
  • 区間平均から Intelが圧倒的に高精度、次点でAMDが高精度

検証4:SFU近似版の結果 絶対誤差可視化

[-pi, pi] の絶対誤差

ULP高値の説明を補完するため 絶対誤差 についても見ておく必要があります。絶対誤差にもなにか傾向があるのでしょうか。ベンダー差はあるのでしょうか。

  • NVIDIA RTX 5090

  • AMD gfx1036

  • Intel UHD 770

NVIDIAでは公式通り 2^(-20.5) ≒ 6.74e-7 以内の誤差で収まっております(グラフでは最大2.5e-7程度と読み取れる)。AMDも同様です。

しかしIntelはいろいろ想定外です。まず絶対誤差のスケールが1e-5と割とでかめに見えます。ただよく考えると問題ないことがわかります。グラフで誤差が急上昇しているのは sin(±pi/2) cos(0,±pi) のところでありこれは ±1 が出力されるときです。
つまり±1 に対し誤差が0.00001程度なので許容できると考えられ、それ以外の区間は ULPの図 や ベンダ横断の比較(抜粋) で見たとおり優秀な相対誤差を維持しています。Intelは他2社と比べどうも相対誤差を重視しているようです。


[-100pi, 100pi] の絶対誤差

cosは割愛

  • NVIDIA RTX 5090

  • AMD gfx1036

  • Intel UHD 770

3社とも特徴がでています。
まずNVIDIAに注目すると引数が2倍になるごとに絶対誤差も2倍に段差状に上昇しています。これは引数縮約(range reduction)が大きく関係していると考えられます。詳細は後述しますが、少なくともNVIDIA GPUにおいてrange reductionがfloat32精度で行われていることがはっきり分かっています。これでは以下の2点のよう精度に問題がおきます。

  • range reductionに使われる2piが360.00001452314234°とわずかに360°からずれている
  • range reduction後の値もfloat32であり仮数部の精度が引数2倍ごとに階段状に悪化

前者の影響は、引数が大きくなるほどリニア(1乗に比例)に誤差の増大を引き起こします。後者の影響は、2の累乗境界(2, 4, 8, 16, 32...)で段差のように見える誤差を引き起こします。

段差についてもう少し突っ込んで話すと・・・
range reductionで0.15915493667125701904をかける操作が入ります。そしてNVIDIAの図では100や200を少し超えたところに段差があります。これに0.159..をかけると16.0、32.0という数字がでてきます。なのでちょうどfloat32の境界に合致します。

AMDは崖のような段差こそないもののリニアに誤差の増大があるので、やはりrange reductionの精度が良くないのでしょう。

そしてIntel GPU!脅威の平坦でありrange reductionの影響がほぼないことがわかります。ULP比較のときもそうでしたが明らかにIntelだけグラフの形が違います。計算方法(や速度)が根本的に違うのでしょうかね。

Intel GPU の native_sin はソフトウェア寄り?(AI調査)
NVIDIA や AMD の GPU では SFU による 1 サイクル近似計算 (MUFU.SIN など) が使われますが、Intel GPU の native_sin は挙動が異なるようです。広い範囲で誤差が平坦なことから 高精度レンジリダクション(Payne–Hanek 型) を用いるソフトウェア実装(またはマイクロコード)と考えられ、高精度だがレイテンシは長めかもしれません。大量に呼ぶカーネルでは性能影響に注意です。


大きな角度でULPが増える理由と引数縮約(range reduction)について

NVIDIA/AMD では引数が大きいほど ULPが単調増加する一方で、Intelは ほぼ増えない という現象が観測されました。
またNVIDIA(AMDも)なら 標準版sin cos 命令を使えばどんな大きな角度でも2ULP以内が保証されています。これら標準版と近似版の違いはなんでしょうか。

この理解を深めるうえで 引数縮約 (argument reduction) を避けて通ることはできません。三角関数の実装はまず、入力角 \theta を所定の基本区間(例:[-\pi, \pi])へ写す縮約を行います。

例えばNVIDIA GPU で SFU版の sin を呼ぶと、SASS(GPUアセンブリ) は

FMUL.RZ   R9, R5, 0.15915493667125701904  // ≈ 1/(2pi) を掛けて「ターン」に変換
MUFU.SIN  R5, R9                           // 1.0 = 1周を引数とする SIN

のような並びが現れます[4:2]。ここで掛けている 0.15915493667125701904 は float32 に丸められた \tfrac{1}{2\pi} で、16進だと 0x3e22f983 です。
この \tfrac{1}{2\pi}値のズレが、1周が何度かを微妙にずらして結果として 大角度ほど位相誤差(=ULP上昇)が積み上がり ます。

なぜ Intel は平坦だったのか、の答え

  • range reductionの係数精度 は、大きな角度域での位相誤差(ULP)に直結
  • Intel は縮約の実装が他2社と比べ高精度でやっている可能性が高く、位相誤差がほぼ見えない → 平坦に見えていた

どれくらいズレるのか? → 360.00001452314234°

1ターンに相当する角度は

\theta_{\text{full}}(c) = \frac{180}{\pi} \cdot \frac{1}{c}

で計算できます。
float32 の丸めにより、実際には次のように わずかに 360° を外れます

定数 (hex) float32 値 1周の角度 [°] 360°との差 [°]
0x3e22f983 0.15915493667125701904296875 360.000014523142306265438 +0.00001452314
0x3e22f984 0.159154951572418212890625 359.999980817510196189747 -0.00001918249

これが 1 周あたりの誤差です。
|\theta| = 1000000\pi(=500000 周)だと、FMUL.RZ R9, R5, 0.15915493667125701904 の時点で 7.2616° の縮約誤差が発生 していることになります。

トリビア

はじめて 90° 以上の誤差が発生するのは6,197,007周目になります。
これはsin(38936943.33)に相当します。この領域では三角関数の出力に意味はないでしょう


結局float32の精度内で縮約していることが問題なので、たとえプログラマー側がsin(x)をsin(x%pi)と書いても結局同様の縮約誤差が発生するので意味がないことを示唆します。これも今後の一つの検証記事にしていきたいと思います。


検証4.1:NVIDIA保証範囲内(±100pi)でAMD,IntelのSFU版絶対誤差最大値

これまでの検証でNVIDIAの保証範囲内(±100pi)でもNVIDIAに優位性があるようには見えませんでした。しかし他ベンダーが絶対誤差の情報を公式発表してないようなのでここではっきりさせておきましょう。

範囲:±2pi

GPU 関数 最大絶対誤差 近似(2の冪) 発生位置 x(参考)
NVIDIA RTX 5090 native_sin 6.556510925e-07 ≈ 2^-20.54 3.288942575454712
native_cos 7.152557373e-07 ≈ 2^-20.42 4.900111198425293
AMD gfx1036 native_sin 4.768371582e-07 ≈ 2^-21.00 6.0378098487854
native_cos 4.172325134e-07 ≈ 2^-21.19 4.467008590698242
Intel UHD 770 native_sin 3.439188004e-05 ≈ 2^-14.83 1.5625
native_cos 3.439188004e-05 ≈ 2^-14.83 0.008290027268230915

参考
Parallel Thread Execution ISA Version 9.0の記載は最大誤差 2^(-20.5) 程度

範囲:±100pi

GPU 関数 最大絶対誤差 近似(2の冪) 発生位置 x(参考)
NVIDIA RTX 5090 native_sin 3.653950989e-05 ≈ 2^-14.74 311.03460693359375
native_cos 3.659725189e-05 ≈ 2^-14.74 312.6540832519531
AMD gfx1036 native_sin 2.465257421e-05 ≈ 2^-15.31 314.1517639160156
native_cos 2.458691597e-05 ≈ 2^-15.31 312.5647888183594
Intel UHD 770 native_sin 3.439188004e-05 ≈ 2^-14.83 1.5625
native_cos 3.439188004e-05 ≈ 2^-14.83 0.008290027268230915

参考
Parallel Thread Execution ISA Version 9.0の記載は最大誤差 2^(-14.7) 程度


検証5:SFU版の超巨大引数

[-1000000pi, 1000000pi] の平均ULPと絶対誤差

保証外の範囲も見ておきましょう。cosも同様の傾向だったのでsinのほうのみ載せます。(githubリポジトリに全結果があります)

  • NVIDIA RTX 5090 平均ULP

  • NVIDIA RTX 5090 絶対誤差

  • AMD gfx1036 平均ULP

  • AMD gfx1036 絶対誤差

  • Intel UHD 770 平均ULP

  • Intel UHD 770 絶対誤差

NVIDIA、AMDは 引数の1乗に比例 ということでやはり良さそうですね。
Intel さん、、、なんですかこれは・・・


検証6:Intel GPUの段差は? → ±32768piで大きく異なる

Intelの結果がとんでもない崖になっていたので、境目がどこにあるか手動二分探索を試みました。

x> 102942.0
  input  : token='102942.0'
  float32: 102942.0  (bits=0x47c90f00)
  u32    : 0x47c90f00

[sin]
  REF.sin:     -0.990592360496521   bits=0xbf7d9776  (71_sin)
  TEST.sin:    -0.990592360496521   bits=0xbf7d9776  (71_sin)
  TEST.native: -0.990592896938324   bits=0xbf7d977f  (ULP=9)  (71_native_sin)

[cos]
  REF.cos:     -0.13684575259685516 bits=0xbe0c214b  (71_cos)
  TEST.cos:    -0.13684575259685516 bits=0xbe0c214b  (71_cos)
  TEST.native: -0.1368417888879776  bits=0xbe0c2041  (ULP=266)  (71_native_cos)

x> 102943.0
  input  : token='102943.0'
  float32: 102943.0  (bits=0x47c90f80)
  u32    : 0x47c90f80

[sin]
  REF.sin:     -0.6503710746765137  bits=0xbf267eb8  (71_sin)
  TEST.sin:    -0.6503710746765137  bits=0xbf267eb8  (71_sin)
  TEST.native: -0.8414709568023682  bits=0xbf576aa4  (ULP=3206124)  (71_native_sin)

[cos]
  REF.cos:     0.7596166729927063   bits=0x3f42763d  (71_cos)
  TEST.cos:    0.7596166729927063   bits=0x3f42763d  (71_cos)
  TEST.native: 0.7596192955970764   bits=0x3f427669  (ULP=44)  (71_native_cos)

すると引数 102942.0 で 266ULP に対し 102943.0 で 3206124 ULPと一気に悪化していることがわかりました。
この102943.0が一体どんな数字なのか、なんとなく3.14159...で割ってみました。

あ!!

おぼろげながら、浮かんできたんです。32768という数字が

±32768piを境にIntel GPUの native_sin native_cosが大きく破綻していることがわかりました。明らかに意味ありげな数字なのでAIに調査させました。

Intel GPU の native_cos/sin が ±32768pi で誤差急増する理由

実測では、Intel GPU の native_cos/sin は ±32768pi を境に ULP が桁外れに跳ね上がります。これは、固定小数点で x·2/pi を評価する高速レンジリダクションの精度限界と考えられます。
固定小数点の小数部が 16bit の場合、正しく象限を判定できる保証域は ±2^(16-1)·pi/2 = ±32768pi までで、これを超えると象限判定がずれ、近似多項式への入力が 90° 以上ずれてしまい、ULP が爆発します。

直接ズバリ書いたベンダー文書はあるか?

少なくとも公開されている Intel や OpenCL の公式ドキュメントに「±32768pi を境界に誤差が増える」と明記したものは見当たりません。ただし、一般的なレンジリダクション理論(Muller『Elementary Functions』など)では 固定小数点 t ビットでの x·2/pi 評価は ±2^(t-2)·pi が保証域と説明されており、musl や glibc の rem_pio2 実装でも入力の大きさで高速経路/高精度経路を切り替えています。したがって ±32768pi という境界は、実装上のしきい値として理論的にも自然です。

  • Jean-Michel Muller, Elementary Functions: Algorithms and Implementation(レンジリダクションと必要ビット数の解説)
  • Brisebarre et al., “A New Range-Reduction Algorithm”(高精度レンジリダクションの理論)
  • musl __rem_pio2f.c(入力規模で高速経路と大入力経路を切替える実装例)
  • glibc reduce_sincos(しきい値付きレンジリダクション)
  • OpenCL Specification – native_* は範囲・誤差ともに実装依存と明記

ということで 16bit固定小数点 が鍵のようです。時間があるときに考察してみたいと思います。

また意外にも "±32768piを境界に誤差が増える" ことを論じたネットソースはないようです。あらためてChat GPTのDeep reserchで検索しましたが具体的な32768という数字に言及した文献はなさそうでした。(私が最初の発見者か!?)

追加検証でIntel HD Graphics 530でも行い、引数114514の正解値 0.1926048398に対し 0.90929 というデタラメな数字を返すことがわかりました。

ということで、おそらくIntel CPU内蔵GPUでは 引数±32768piを境に精度が壊滅的に悪化するようです。

興味がある人はUnityのShader等をIntel GPU上で動作させてこのような結果が現れるか確認してみてください!


検証7:標準版sin cosの誤差は? → 基本2ULP以内

NVIDIA GPU では標準版sin cosの誤差が 2ULP以内 と明言されていますが AMD/Intel でそのような公開文章は見つけられませんでした。とはいえおそらく2ULP以内に近い高精度を提供してくれていると考えられます。

doubleの結果を正解として何ULPの誤差があるか検証しました。ここでは非正規化数とinfとnanを除いた4,261,412,866個の正規化数を引数として用いました。
https://github.com/toropippi/GPU_SinCosULP/blob/main/2_結果表示/check_highprec_ulp.py

GPU(sin/cos) 最大ULP その時の引数 平均ULP
NVIDIA RTX5090(sin) 1 0.00044363297638483346 0.099970
NVIDIA RTX5090(cos) 2 71476.0703125 0.100751
AMD gfx1036(sin) 2 22.241199493408203 0.099740
AMD gfx1036(cos) 2 1.3181171417236328 0.101446
Intel UHD 770(sin) 4294967296 -0.0 (bits=0x80000000) 1.115203
Intel UHD 770(cos) 2 0.00023394824529532343 0.128965

AMDは想定通り 最大2ULP 以内におさまっており、NVIDIAも公表通り 2ULP 以内におさまっています。平均ULPもとんでもなく優秀ですね!

Intelはsinで+-0.0につまづいてULPが馬鹿高くなっているので忖度版を作成して取り直しました。
出力値が-0.0のものは+0.0に正規化しています。

GPU(sin/cos) 最大ULP その時の引数 平均ULP
Intel UHD 770(sin) 2 1.4148300886154175 0.107329
Intel UHD 770(cos) 2 0.00023394824529532343 0.128965

ということでIntelでも 2 ULP 以内におさまって平均ULPも0.1 と、他2ベンダーと同等の性能であることがわかりました!


検証8:±Inf, ±0.0, NaN, 非正規化数の扱い

±Inf 入力時の出力結果(sin / cos / native_sin / native_cos)

GPU 入力 関数 結果 (print 表現) bit 表現
NVIDIA RTX 5090 +inf sin (標準) nan 0x7fffffff
native_sin nan 0x7fffffff
cos (標準) nan 0x7fffffff
native_cos nan 0x7fffffff
-inf sin (標準) nan 0x7fffffff
native_sin nan 0x7fffffff
cos (標準) nan 0x7fffffff
native_cos nan 0x7fffffff
AMD gfx1036 +inf sin (標準) nan 0x7fc00000
native_sin nan 0xffc00000
cos (標準) nan 0x7fc00000
native_cos nan 0xffc00000
-inf sin (標準) nan 0x7fc00000
native_sin nan 0xffc00000
cos (標準) nan 0x7fc00000
native_cos nan 0xffc00000
Intel UHD 770 +inf sin (標準) nan 0x7fc00000
native_sin nan 0xffc00000
cos (標準) nan 0x7fc00000
native_cos nan 0xffc00000
-inf sin (標準) nan 0x7fc00000
native_sin nan 0xffc00000
cos (標準) nan 0x7fc00000
native_cos nan 0xffc00000
  • NVIDIA は標準版・近似版ともに nan の bit 表現が 0x7fffffff になっている(quiet NaN の全ビット1?)
  • AMD/Intel は標準版が 0x7fc00000(典型的な quiet NaN)、近似版は 0xffc00000(符号ビットが立った NaN
  • 符号付き NaN を返す実装があるのが興味深い(AMD/Intel の native)

±0 入力時の出力結果(sin / cos / native_sin / native_cos)

GPU 入力 関数 結果 (print 表現) bit 表現
NVIDIA RTX 5090 +0.0 sin (標準) 0.0 0x00000000
native_sin 0.0 0x00000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
-0.0 sin (標準) -0.0 0x80000000
native_sin -0.0 0x80000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
AMD gfx1036 +0.0 sin (標準) 0.0 0x00000000
native_sin 0.0 0x00000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
-0.0 sin (標準) -0.0 0x80000000
native_sin -0.0 0x80000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
Intel UHD 770 +0.0 sin (標準) 0.0 0x00000000
native_sin 0.0 0x00000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
-0.0 sin (標準) 0.0(※符号消失) 0x00000000
native_sin -0.0 0x80000000
cos (標準) 1.0 0x3f800000
native_cos 1.0 0x3f800000
  • NVIDIA / AMD は標準版・近似版ともに -0.0 を維持して返却(符号付き 0 を正しく再現
  • Intel標準版 sin-0.0 → +0.0 に変換されており、符号ビットが落ちる実装になっている
  • 一方 Intel native_sin では -0.0 を保持している

NaN 入力時の出力結果(sin / cos / native_sin / native_cos)

入力 NaN の bit 表現: 0x7fc00000

GPU 入力 関数 結果 (print 表現) bit 表現
NVIDIA RTX 5090 nan sin (標準) nan 0x7fffffff
native_sin nan 0x7fffffff
cos (標準) nan 0x7fffffff
native_cos nan 0x7fffffff
AMD gfx1036 nan sin (標準) nan 0x7fc00000
native_sin nan 0x7fc00000
cos (標準) nan 0x7fc00000
native_cos nan 0x7fc00000
Intel UHD 770 nan sin (標準) nan 0x7fc00000
native_sin nan 0x7fc00000
cos (標準) nan 0x7fc00000
native_cos nan 0x7fc00000
  • NVIDIA は標準/近似 ともに 0x7fffffff(実装依存の quiet NaN ペイロード)を返し、AMD/Intel は典型的な quiet NaN 0x7fc00000 を返している

非正規化数 入力時(SFU版)

native_sinの結果と正解値をプロットしています。

  • NVIDIA RTX 5090

  • AMD gfx1036

  • Intel UHD 770

どうもNVIDIAとIntelは0.0を出力しているようです。一方AMDはおそらく重なっている?
手動で確認してみました。

まずNVIDIA RTX 5090のとき

x> 0x00000100
  input  : token='0x00000100'
  float32: 3.587324068671532e-43  (bits=0x00000100)
  u32    : 0x00000100

[sin]
  REF.sin:     3.587324068671532e-43 bits=0x00000100  (0_sin)
  TEST.sin:    3.587324068671532e-43 bits=0x00000100  (0_sin)
  TEST.native: 0.0                  bits=0x00000000  (ULP=256)  (0_native_sin)

[cos]
  REF.cos:     1.0                  bits=0x3f800000  (0_cos)
  TEST.cos:    1.0                  bits=0x3f800000  (0_cos)
  TEST.native: 1.0                  bits=0x3f800000  (ULP=0)  (0_native_cos)

x> 0x007FFFFF
  input  : token='0x007FFFFF'
  float32: 1.1754942106924411e-38  (bits=0x007fffff)
  u32    : 0x007fffff

[sin]
  REF.sin:     1.1754942106924411e-38 bits=0x007fffff  (0_sin)
  TEST.sin:    1.1754942106924411e-38 bits=0x007fffff  (0_sin)
  TEST.native: 0.0                  bits=0x00000000  (ULP=8388607)  (0_native_sin)

[cos]
  REF.cos:     1.0                  bits=0x3f800000  (0_cos)
  TEST.cos:    1.0                  bits=0x3f800000  (0_cos)
  TEST.native: 1.0                  bits=0x3f800000  (ULP=0)  (0_native_cos)

x> 0x807FFFFF
  input  : token='0x807FFFFF'
  float32: -1.1754942106924411e-38  (bits=0x807fffff)
  u32    : 0x807fffff

[sin]
  REF.sin:     -1.1754942106924411e-38 bits=0x807fffff  (128_sin)
  TEST.sin:    -1.1754942106924411e-38 bits=0x807fffff  (128_sin)
  TEST.native: -0.0                 bits=0x80000000  (ULP=8388607)  (128_native_sin)

[cos]
  REF.cos:     1.0                  bits=0x3f800000  (128_cos)
  TEST.cos:    1.0                  bits=0x3f800000  (128_cos)
  TEST.native: 1.0                  bits=0x3f800000  (ULP=0)  (128_native_cos)

このように入力を0x00000100,0x007FFFFF(最大の非正規化数),0x807FFFFFとしてみましたが、native_sinの出力が±0.0に張り付いています!符号情報だけ保持しあとは指数部、仮数部を0にしているだけのようです。
非正規化数は非常に小さい値なのでこれでもいいかもしれませんが。なおIntelも同様の結果でした。

次にAMD gfx1036では

x> 0x00000100
  input  : token='0x00000100'
  float32: 3.587324068671532e-43  (bits=0x00000100)
  u32    : 0x00000100

[sin]
  REF.sin:     3.587324068671532e-43 bits=0x00000100  (0_sin)
  TEST.sin:    3.587324068671532e-43 bits=0x00000100  (0_sin)
  TEST.native: 3.615350037958028e-43 bits=0x00000102  (ULP=2)  (0_native_sin)

[cos]
  REF.cos:     1.0                  bits=0x3f800000  (0_cos)
  TEST.cos:    1.0                  bits=0x3f800000  (0_cos)
  TEST.native: 1.0                  bits=0x3f800000  (ULP=0)  (0_native_cos)

+-0.0ではなくちゃんとした値がでてきました!2 ULPといい感じの精度です。
非正規化数に限定して 誤差ULP をリニアスケールでとってみると

これまたいい感じの精度で近似されていることがわかりました。最大誤差は7 ULPでした。正規化数のとき同様、AMDの0.0近傍のnative_sinには謎のこだわりを感じます。

まとめると

  • NVIDIA / Intel は常に ±0.0 を返すザル設計
  • AMD標準版 sin 相当の精度で 2~3 ULP 程度の平均誤差
  • ちなみにこの領域ではsin(x)=xとしてよく、標準版sinでも単に入力xをそのまま返してくるので、なぜそのような実装がかなわなかったのか・・?

非正規化数 入力時(標準sin版)

上記の検証中、標準版では全例 sin(x)=x として結果が出ているように見えました。
16,777,214個の非正規化数を全て調査したところ、

全ベンダーで sin(x)=x として出力されていることがわかりました。


まとめ集(セクションごと)と総括

長い記事になったのでAIまとめ


序章-方法 — まとめ

  • 主要 3 ベンダ(NVIDIA/AMD/Intel)の複数 GPU を対象
  • sin/cos(標準)と native_sin/native_cos(SFU)を引数 32bit 全探索
  • NaN/Inf/非正規含めすべてを評価、総データ未圧縮 481GB 規模

検証1:同一ベンダ(NVIDIA)内差 — まとめ

  • 複数世代/OS で bit 単位まで完全一致
  • NVIDIA の 近似/標準版の決定論的挙動が確認できた

検証2:SFU 版の誤差傾向 — まとめ

  • 0・pi/2・pi 近傍で ULP が跳ねやすい(出力が 0/±1 近傍=相対誤差が効く)
  • AMD は native_sin0 近傍が強い傾向。Intel は善戦
  • sin(x+pi/2) ≈ cos(x) の関係から、ULP の山/谷が入れ替わる現象を説明

検証3:ULP 可視化(±pi / ±10pi / ±100pi)— まとめ

  • NVIDIA/AMD は pi 周期で鋭いピーク、Intel は pi/2 周期の特徴
  • NVIDIA/AMD は入力が大きくなると 平均 ULP がほぼ線形増加
  • Intel は 入力増大の影響が小さく、多くの区間で 0~1 ULP と高精度

検証4:絶対誤差可視化 — まとめ

  • NVIDIA/AMD は 公式値どおり(±2pi→2^-20.5、±100pi→2^-14.7 近辺)
  • NVIDIA は **段差(2 倍ごとに誤差 2 倍)**が現れ、range reduction を f32 で実施している痕跡
  • Intel は 平坦(±1 近傍で 1e-5 程度の絶対誤差)で、相対誤差最適寄りの設計が示唆

引数縮約(range reduction)の影響・1 周角度のズレ — まとめ

  • NVIDIA の SASS では FMUL x (1/2pi)_f32 → MUFU.SIN係数が f32 丸めで位相ズレを誘発
  • 0x3e22f983/984 のわずかな違いで 1 周角度が ±数 10 ミリ秒角ズレ
  • 大角度ほど**累積位相誤差(=ULP 上昇)**が大きくなる構造。Intel は縮約が高精度っぽく影響が小さい

検証4.1:±100pi 内での絶対誤差最大値 — まとめ

  • NVIDIA:実測は仕様と整合(±2pi/±100pi ともに公称近辺)
  • AMD:±2pi では 2^-21 前後、±100pi でも 2^-15.3 と良好(公式未開示の実測知見)
  • Intel:範囲を広げても ≈2^-14.83 に張り付き(平坦)。高精度縮約の可能性を裏付け

検証5:保証外(±100 万pi)— まとめ

  • NVIDIA/AMD は 平均 ULP がほぼ線形に増加し「予想どおり」
  • Intel は 形状が異常(平坦+ひどい崖)で、±32768piを超えると壊滅的

検証6:Intel の段差(±32768pi)— まとめ

  • Intel の native±32768pi を境に ULP が激増
  • t=16bit 固定小数点x·2/pi を評価する高速縮約の保証域(理論的に自然)と合致
  • ベンダ公式に直書きは見当たらず。本稿は現象の初出級の実測整理か?

検証7:標準版(ソフトウェア)— まとめ

  • NVIDIA/AMD は ≤ 2 ULP を満たす実測。平均 ULP も ≈0.1 と優秀
  • Intel も ±0 の符号正規化後は ≤ 2 ULP に収束(平均 ≈0.1)
  • 標準版は想定通り

検証8:±Inf / ±0.0 / NaN / 非正規 — まとめ

  • NaN:NVIDIA は 0x7fffffff 系、AMD/Intel は 0x7fc00000(native で符号付き NaN を返す実装も)
  • ±0.0:NVIDIA/AMD は -0.0 を保持、Intel の標準版のみ sin(-0.0)符号落ち
  • 非正規(subnormal) のSFU:NVIDIA/Intel の native_sin は全件 ±0 に潰れる。AMD は sin(x)≈x を保ち 0~数 ULP と優秀(?)
  • 非正規(subnormal) の標準版sin:全ベンダで全件 sin(x)=x が成り立つ

総括(全体のまとめ)

  • 近似版を使ってはいけないわけではなく、むしろ特性を理解した上で使うべき!
  • 近似版は高速だが、range reduction の実装差で誤差の出方に明確なベンダ差がある。
    • NVIDIA:良くも悪くも教科書的な実装が見られた。pi 周期でピーク、入力増大で平均 ULP が線形増。f32の痕跡(段差)など
    • AMD:謎の0.0近傍のこだわり。しかし±piではその現象がなくなるので本当に0.0近傍のみ
    • Intel:近似版でもrange reductionが高精度の可能性。pi/2 周期の特徴。しかし引数が±32768pi を超えると事実上使えない(固定小数点縮約の保証域を示唆)
  • 標準版は 3 ベンダとも実測で ≤ 2 ULP に収束(Intel は +0.0 正規化前提)
  • 特殊値では、NaN/±0/非正規の扱いに方針差
脚注
  1. GPU命令レベルで読むマンデルブロー集合 ~Nsight ComputeとSASS解説~ ↩︎

  2. Fermi_Architecture_Whitepaper ↩︎

  3. Hexops: DXC の解説記事 ↩︎

  4. CUDAの__sinf()はSASSでFMUL → MUFU.SINに展開される ↩︎ ↩︎ ↩︎

Discussion