3

私はまだSSEの使用にかなり慣れていない2*Piため、次数の倍精度入力のモジュロを実装しようとして1e8います(その結果は、いくつかのベクトル化された三角計算に供給されます)。

mod(x, 2*Pi) = x - floor(x/(2*Pi))*2*Piコードでの私の現在の試みは、次のようなアイデアに基づいています。

#define _PD_CONST(Name, Val)                                            \
static const double _pd_##Name[2] __attribute__((aligned(16))) = { Val, Val }  

_PD_CONST(2Pi, 6.283185307179586);  /* = 2*pi  */  
_PD_CONST(recip_2Pi, 0.159154943091895); /* = 1/(2*pi)  */

void vec_mod_2pi(const double * vec, int Size, double * modAns)
{
    __m128d sse_a, sse_b, sse_c;
    int i;
    int k = 0;
    double t = 0;

    unsigned int initial_mode;
    initial_mode = _MM_GET_ROUNDING_MODE();

    _MM_SET_ROUNDING_MODE(_MM_ROUND_DOWN);

    for (i = 0; i < Size; i += 2)
    {
        sse_a = _mm_loadu_pd(vec+i);
        sse_b = _mm_mul_pd( _mm_cvtepi32_pd( _mm_cvtpd_epi32( _mm_mul_pd(sse_a, *(__m128d*)_pd_recip_2Pi) ) ), *(__m128d*)_pd_2Pi);
        sse_c = _mm_sub_pd(sse_a, sse_b);
        _mm_storeu_pd(modAns+i,sse_c);
    }

    k = i-2;
    for (i = 0; i < Size%2; i++)
    {
        t = (double)((int)(vec[k+i] * 0.159154943091895)) * 6.283185307179586;
        modAns[k+i] = vec[k+i] - t;
    }

    _MM_SET_ROUNDING_MODE(initial_mode);
}

残念ながら、これは現在、いくつかNaNの回答とともに多くの回答を返しています (私が目指していた->1.128e119の範囲外のものもあります!)。私が間違っているのは、.02*Pifloor

どこが間違っていて、それを改善する方法を誰かが提案できますか?

PSそのコードの形式について申し訳ありません。ここに質問を投稿したのは初めてで、コードブロック内に空の行を表示して読みやすくすることができないようです。

4

2 に答える 2

7

なんらかの精度が必要な場合、単純なアルゴリズムはひどく悪いです。正確な範囲縮小アルゴリズムについては、Ng et al., ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit (Wayback Machine から入手可能: 2012-12-24 )などを参照してください。

于 2012-07-20T09:08:05.067 に答える
1

大きな引数の場合、 Hayne-Panek アルゴリズムが通常使用されます。ただし、Hayne-Panek の論文は非常に読みにくいため、より理解しやすい説明については、Handbook of Floating-Point Arithmetic の第 11 章を参照することをお勧めします。

于 2012-07-20T11:00:43.797 に答える