私はまだ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
の範囲外のものもあります!)。私が間違っているのは、.0
2*Pi
floor
どこが間違っていて、それを改善する方法を誰かが提案できますか?
PSそのコードの形式について申し訳ありません。ここに質問を投稿したのは初めてで、コードブロック内に空の行を表示して読みやすくすることができないようです。