前に完全精度expf関数を作った（高速な完全精度 expf 関数の作り方 - よーる）ので、今度は完全精度sincosf関数を作っていきます。

作戦

まず、引数の絶対値が小さい時のsin関数とcos関数を（taylor展開などを利用して）作ります。引数の絶対値が大きい場合、引数 $x$ を $x = t +\frac\pi 2 k$ と分解します。その後 $k$ の値に応じて先に作ったsin関数やcos関数を使い分ければ答えが求まります。

ただし、この分解は簡単にできることではありません。今回は、ある程度簡単に実行可能な、引数の絶対値がそれほど大きくない（4×10⁸程度以下）場合のみを作ることにします。

ベースとなる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);
}

$\frac\pi 2$ で割ったあまりを求める

引数 $x$ を $x = t + \frac\pi 2 k$ と分解するためには、以下のような計算を行います。これは、完全精度expf関数を作る時に行った、引数 $x$ を $x = t + \frac{\log 2}{2^N}k$ と分解するのと同様の手順です。

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ビットは0
harf_pi_hiの値とharf_pi_loの値を足し合わせると、 $\frac\pi 2$ に近い
inv_harf_piの値は $\frac 2\pi$ に近い

harf_pi_hiの下位Nビットを0にするのは、k倍したときに誤差が生じないようにするためです。ここでN=28を選べば、先のコードになります。 N=28の時、xの絶対値が4.2e+8程度以下であればassertが満たされます。

なぜN=28を選んだか

なぜN=28を選ぶ必要があるのでしょうか。例えば、Nをもっと大きくすれば、対応できるxの範囲が広がりよさそうに思われます。しかし、実際には、Nを大きくすると $\frac\pi 2$ の有効桁数が小さくなり、それに起因する誤差が増加し、完全精度を達成できなくなります。このことを説明してみます。

doubleの精度は53ビットですが、二つのdoubleの和で表された数の精度は107ビットです。これは、下位を担当するdoubleの符号ビットも有効活用されるためです。ここで、Nビットを強制的に0にしなければならないので、 $\frac\pi 2$ の実効的な精度は107-Nビットです。

$t = x - \frac\pi 2 k$ を計算するとき、桁落ちが発生します。 $\frac\pi 2 k$ と $x$ とで打ち消しあう部分の長さはおおよそN+1ビットなので、最終結果には106-2Nビットしか残りません。

N=29の場合、小数部には48ビットしか使われないということになります。これはdoubleの精度を下回っています。実際、0x1.8db252p+25を $\frac\pi 2$ で割ったあまりを求めてみると、正確には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に丸められ、正しくなくなります。