前に完全精度expf関数を作った(高速な完全精度 expf 関数の作り方 - よーる)ので、今度は完全精度sincosf関数を作っていきます。
作戦
まず、引数の絶対値が小さい時のsin関数とcos関数を(taylor展開などを利用して)作ります。 引数の絶対値が大きい場合、引数をと分解します。その後の値に応じて先に作ったsin関数やcos関数を使い分ければ答えが求まります。
ただし、この分解は簡単にできることではありません。 今回は、ある程度簡単に実行可能な、引数の絶対値がそれほど大きくない(4×108程度以下)場合のみを作ることにします。
ベースとなるTaylor展開ベースの実装
まず、sin関数とcos関数をtaylor展開で求めるコードを書きます。double
精度で求めて最後にfloat
に丸める作戦なので、そこまで精度は必要ではありません。
double cos_taylor(double x) { double x2 = x * x; double a = -1./87178291200; a = std::fma(a, x2, 1./479001600); a = std::fma(a, x2, -1./3628800); a = std::fma(a, x2, 1./40320); a = std::fma(a, x2, -1./720); a = std::fma(a, x2, 1./24); a = std::fma(a, x2, -1./2); return std::fma(a, x2, 1); } double sin_taylor(double x) { double x2 = x * x; double a = -1./1307674368000; a = std::fma(a, x2, 1./6227020800); a = std::fma(a, x2, -1./39916800); a = std::fma(a, x2, 1./362880); a = std::fma(a, x2, -1./5040); a = std::fma(a, x2, 1./120); a = std::fma(a, x2, -1./6); return std::fma(a, x2 * x, x); }
で割ったあまりを求める
引数をと分解するためには、以下のような計算を行います。 これは、完全精度expf関数を作る時に行った、引数をと分解するのと同様の手順です。
constexpr double harf_pi_hi = 0x1.921fb5p+0; constexpr double harf_pi_lo = 0x1.110b4611a6263p-26; constexpr double inv_harf_pi = 0x1.45f306dc9c883p-1; double k = std::fma(x, inv_harf_pi, 0x3.p+51) - 0x3.p+51; assert( -0x1.p+28 <= k && k <= 0x1.p+28 ); double t = std::fma(-harf_pi_lo, k, std::fma(-harf_pi_hi, k, x));
ここで、harf_pi_hi
、harf_pi_lo
、inv_harf_pi
は、以下の条件を満たすdouble
の定数です。
harf_pi_hi
の下位Nビットは0harf_pi_hi
の値とharf_pi_lo
の値を足し合わせると、に近いinv_harf_pi
の値はに近い
harf_pi_hi
の下位Nビットを0にするのは、k
倍したときに誤差が生じないようにするためです。
ここでN=28を選べば、先のコードになります。
N=28の時、x
の絶対値が4.2e+8
程度以下であればassert
が満たされます。
なぜN=28を選んだか
なぜN=28を選ぶ必要があるのでしょうか。例えば、Nをもっと大きくすれば、対応できるx
の範囲が広がりよさそうに思われます。
しかし、実際には、Nを大きくするとの有効桁数が小さくなり、それに起因する誤差が増加し、完全精度を達成できなくなります。
このことを説明してみます。
double
の精度は53ビットですが、二つのdouble
の和で表された数の精度は107ビットです。これは、下位を担当するdouble
の符号ビットも有効活用されるためです。
ここで、Nビットを強制的に0にしなければならないので、の実効的な精度は107-Nビットです。
を計算するとき、桁落ちが発生します。ととで打ち消しあう部分の長さはおおよそN+1ビットなので、最終結果には106-2Nビットしか残りません。
N=29の場合、小数部には48ビットしか使われないということになります。これはdouble
の精度を下回っています。
実際、0x1.8db252p+25
をで割ったあまりを求めてみると、正確には0x1.5292b563fb139p-5
程度になるべきところ、0x1.5292b563fb120p-5
と求まってしまいます。なんと25ULPも誤差が発生しています。
運の悪いことにcos(0x1.8db252p+25) = 0x1.527a08ffffff0p-5
とfloat
に丸めるときの丸め境界から16ULP程度しか離れていないため、ここの計算を間違えてしまいます。
よって、このような問題を発生させないためにも、N=28よりも大きくすることはできません(N=28であっても、運よくこのような問題が発生していないだけであり、やや精度不足であることには変わりありません)。
絶対値がそれほど大きくない場合の実装
void sincos_reduction(double x, double* s, double* c) { constexpr double harf_pi_hi = 0x1.921fb5p+0; constexpr double harf_pi_lo = 0x1.110b4611a6263p-26; constexpr double inv_harf_pi = 0x1.45f306dc9c883p-1; double k = std::fma(x, inv_harf_pi, 0x3.p+51) - 0x3.p+51; assert( -0x1.p+28 <= k && k <= 0x1.p+28 ); double t = std::fma(-harf_pi_lo, k, std::fma(-harf_pi_hi, k, x)); std::int32_t k_ = static_cast<std::int32_t>(k); switch(k_ & 3) { case 0: *s = sin_taylor(t), *c = cos_taylor(t); break; case 1: *s = cos_taylor(t), *c = -sin_taylor(t); break; case 2: *s = -sin_taylor(t), *c = -cos_taylor(t); break; case 3: *s =-cos_taylor(t), *c = sin_taylor(t); break; } }
8/28 追記
入力が±0x1.33333p+13
の時、sin(x) ~ ∓0x1.63f4bafffffffa5cp-2
であり、完全精度なら∓0x1.63f4bap-2
を返すべきですが、上記コードは∓0x1.63f4bcp-2
を返すので完全精度になっていませんでした。
これはテイラー展開の誤差由来やリダクション部分の誤差由来ではなく、double
を経由していること、つまり二回丸めていることが原因でした。sin(x) ~ ∓0x1.63f4bafffffffa5cp-2
はdouble
に丸めると∓0x1.63f4bb0000000p-2
に丸められ、float
の丸め境界上の値に丸まってしまいます。これをfloat
に丸めると偶数丸めになるため最終桁が0になる∓0x1.63f4bcp-2
に丸められ、正しくなくなります。