正弦関数の例を使用した浮動小数点数の正確で高速な計算。パート3：固定点

講義のサイクルを続けます（パート1とパート2）。パート2では、libmライブラリ内にあるものを確認しました。この作業では、do_sin関数を少し変更して、精度と速度を向上させます。この関数をもう一度引用しますdo_sin）：







前の記事、パート132-145に示されているように。 [0.126、0.855469]の範囲内のxに対して実行されます。上手。与えられた制限内で、より正確に、そしておそらくより速くなる関数を書いてみましょう。



私たちがそれを使用する方法はかなり明白です。計算の精度を拡張して、小数点以下の桁数を増やす必要があります。明らかな解決策は、long doubleタイプを選択し、それを数えてから、元に戻すことです。精度の面では解決策は良いはずですが、パフォーマンスの面では問題があるかもしれません。結局のところ、long doubleはかなりエキゾチックなタイプのデータであり、最新のプロセッサでのサポートは優先事項ではありません。 x86_64では、このデータタイプのSSE / AVX命令は機能しません。数学的コプロセッサーは「吹き飛ばされ」ます。



では、何を選ぶべきですか？引数と関数の限界をもう一度見てみましょう。



それらは1.0領域にあります。それら。実際、浮動小数点は必要ありません。関数を計算するときに64ビットの整数を使用しましょう。これにより、元の精度に10〜11ビットが追加されます。これらの番号を操作する方法を理解しましょう。この形式の数値はa / dとして表されます。ここで、aは整数であり、dはすべての変数に対して定数を選択し、コンピューターのメモリではなく「メモリに」格納する除数です。以下は、そのような番号のいくつかの操作です。

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

ご覧のとおり、複雑なことは何もありません。最後の式は、任意の整数による乗算を示しています。サイズNの2つの符号なし整数変数を乗算した結果は、多くの場合、2 * Nまでのサイズの数になることにも注意してください。追加すると、最大1ビットのオーバーフローが発生する可能性があります。



除数dを選択してみましょう。明らかに、バイナリの世界では、除算せずに、たとえばレジスタを移動するために、2の累乗として選択するのが最善です。 2の累乗を選択する必要がありますか？乗算機の命令でヒントを見つけてください。たとえば、x86システムの標準MUL命令は、2つのレジスタを乗算し、結果を2つのレジスタにも書き込みます。ここで、レジスタの1つは結果の「上部」で、2番目は下部です。



たとえば、64ビットの数値が2つある場合、結果は128ビットの数値が2つの64ビットレジスタに書き込まれます。RH-「大文字」の場合、RL-「小文字」の場合¹と呼びましょう。次に、数学的に、結果は次のように書くことができます。$$$R=R_H \cdot 2^{64}+R_L$..。ここで、上記の式を使用して、次の乗算を記述します。$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

そして、これらの2つの固定小数点数を乗算した結果がレジスタであることがわかります $$$R=R_H$..。



Aarch64システムの場合はさらに簡単です。「UMULH」命令は2つのレジスタを乗算し、乗算の「上位」部分を3番目のレジスタに書き込みます。



じゃあ。固定ポイント番号を指定しましたが、まだ問題が1つあります。負の数。テイラーシリーズでは、展開は可変符号で行われます。この問題に対処するために、Goner法に従って多項式を計算するための式を次の形式に変換します。

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

数学的に元の式とまったく同じであることを確認してください。しかし、各括弧にはいくつかの形式があります$$$1/(2n + 1)! - x^2\cdot(\cdots)$常にポジティブ。それら。この変換により、式を符号なし整数として評価できます。

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

どこ $$$tsx[i]=1/i!$固定ポイント形式。今回は、便宜上、すべてのコードをfast_sine githubに投稿し、clangおよびarmとの互換性のためにquadmathを削除しました。そして、誤差の計算方法を少し変更しました。



標準正弦関数と固定点関数の比較を以下の2つの表に示します。最初の表は計算精度を示しています（x86_64とARMで完全に同じです）。2番目の表はパフォーマンスの比較です。







テスト中に、MPFRライブラリを使用して「真の」正弦値が計算されました..。最大ULP値は、「真の」値からの最大偏差と見なされました。エラーの割合-私たちまたはlibmバイナリによる正弦関数の計算値が、2倍の正弦値に切り上げられた値と一致しなかったケースの数。偏差の平均値は、計算エラーの「方向」を示します。つまり、値の過大評価または過小評価です。表からわかるように、私たちの関数は正弦値を過大評価する傾向があります。これは修正できます！ tsx値はテイラーシリーズの係数と正確に等しくなければならないと誰が言いましたか？かなり明白なアイデアは、近似の精度を向上させ、エラーの定数成分を取り除くために、係数の値を変更することを示唆しています。このようなバリエーションを正しく作成することは非常に困難です。しかし、私たちは試すことができます。例を見てみましょうtsx係数の配列（tsx [3]）の4番目の値で、最後の数値aをfに変更します。プログラムを再起動して、精度（sin_e7a）を確認しましょう。ほら、少しですが増えました！このメソッドをピギーバンクに追加します。



それでは、パフォーマンスを見てみましょう。テストのために、手元にあるi5 mobileと、わずかにオーバークロックされた4番目のラズベリー（Raspberry PI 4 8GB）、両方のシステムのUbuntu 20.04x64ディストリビューションのGCC10を使用しました。







私はこれらの測定でより正確なふりをしません。プロセッサの負荷に応じて、数十パーセントの変動が可能です。主な結論はこのようにすることができます。整数演算に切り替えても、最新のプロセッサ²ではパフォーマンスが向上しません。最新のプロセッサの驚異的な数のトランジスタにより、複雑な計算をすばやく実行できます。ただし、Intel Atomなどのプロセッサや弱いコントローラでは、このアプローチによってパフォーマンスが大幅に向上すると思います。どう思いますか？



このアプローチにより精度が向上しましたが、この精度の向上は有用というよりも興味深いようです。パフォーマンスの観点から、このアプローチは、たとえばIoTに見られます。しかし、高性能コンピューティングの場合、それはもはや主流ではありません。今日の世界では、SSE / AVX / CUDAは並列関数計算を使用することを好みます。そして浮動小数点演算で。 MUL機能の並列アナログはありません。機能自体はむしろ伝統への賛辞です。



次の章では、計算にAVXを効果的に使用する方法について説明します。繰り返しになりますが、libmコードにアクセスして、改善してみましょう。



¹私が知っているプロセッサには、そのような名前のレジスタはありません。たとえば、名前が選択されています。

²ここで、私のARMには最新バージョンの数学コプロセッサが装備されていることに注意してください。プロセッサが浮動小数点計算をエミュレートした場合、結果は劇的に異なる可能性があります。