正弦余弦モジュラ拡張精度演算

Question

私は、正弦/余弦の多くの実装では、いわゆる拡張モジュール精度演算を見てきました。しかし、それは何のためですか？ たとえば、cephes implemetationでは、[0、pi/4]の範囲に縮小した後で、精度を向上させるためにこのモジュラ精度演算を行っています。コード】以下

：DP1、DP2及びDP3は、いくつかのハードコーディング係数である

z = ((x - y * DP1) - y * DP2) - y * DP3; 

。 これらの係数を数学的に見つけるにはどうすればよいですか？私はbig numのための "モジュラーエクステンション算術"の目的を理解しましたが、ここでその正確な目的は何ですか？

Answer 1

三角関数の引数削減の文脈では、本書で紹介されているテクニックであるCody-Waite議論の削減があります：William J. CodyとWilliam Waite、ソフトウェアマニュアル、 Prentice-Hall、1980。中間計算ではsubtractive cancellationであるにもかかわらず、一定の大きさまでの議論に対して、正確な減少議論を達成することが目標です。この目的のために、関連定数は、複数の減少する大きさの和（ここでは、DP1,DP2,DP3）を使用することによって、ネイティブ精度を超えるで表され、最下位のものを除くすべての中間生成物が丸め誤差なしで計算することができます。

例として、IEEE-754 binary32（単精度）のsin（113）の計算を考えてみましょう。典型的な引数の削減は、概念的にはi=rintf(x/(π/2)); reduced_x = x-i*(π/2)を計算するでしょう。 π/ 2に最も近い数字のbinary32は0x1.921fb6p+0です。 i=72を計算すると、0x1.c463acp+6に丸められ、これは引数x=0x1.c40000p+6に近いです。減算の間、いくつかの先頭のビットはキャンセルされ、reduced_x = -0x1.8eb000p-4で終了します。再正規化によって導入された末尾のゼロに注意してください。これらのゼロビットは有用な情報を持たない。縮小された引数に正確な近似を適用するとsin(x) = -0x1.8e0eeap-4、真の結果は-0x1.8e0e9d39...p-4です。大きな相対誤差と大きなulpエラーで巻き上げます。

2段階のCody-Waite引数削減を使用してこれを修正できます。たとえば、pio2_hi = 0x1.921f00p+0とpio2_lo = 0x1.6a8886p-17を使用できます。末尾の8つのゼロビットは、単精度表現でpio2_hiであることに注意してください。これにより、8ビット整数iを乗算することができます。i * pio2_hiの単精度数が得られます。 ((x - i * pio2_hi) - i * pio2_lo)を計算すると、reduced_x = -0x1.8eafb4p-4となるので、sin(x) = -0x1.8e0e9ep-4という非常に正確な結果が得られます。

定数に合計を分割する最も良い方法は、与えられた引数範囲に対して減法的な相殺の対象となるビットの最大数（対数πの整数倍に基づいて処理する必要があるiの大きさに依存します）/2は整数に達することができます）、およびパフォーマンスの考慮事項。典型的な実際の使用事例には、2〜4段のCody-Waite削減スキームが含まれます。融合多重加算（FMA）が利用可能であるため、後続のゼロビットが少ない構成定数を使用することができます。この論文を参照してください：Sylvie Boldo、Marc Daumas、Ren-Cang Li、「正式に議論された、積和加算による引数の削減」 IEEE Transactions onコンピュータ,58：1139-1145,2009を参照してください。fmaf()を使用した作業例については、one of my previous answersのコードを参照してください。