Modulo 2 * Pi使用SSE/SSE2

我對使用SSE還很新，並且試圖爲1e8（其結果將被饋送到一些向量化的trig）的雙精度輸入執行模的2*Pi計算）。Modulo 2 * Pi使用SSE/SSE2

我的代碼當前的嘗試是根據各地的想法，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的雙倍詮釋到雙倍的轉換。

任何人都可以提出我哪裏出了問題，以及如何改善它？

P.S.對於這段代碼的格式感到抱歉，這是我第一次在這裏發佈一個問題，似乎無法讓代碼塊中的空行讓它變得可讀。

來源

2012-07-20 BigA

也許我應該讓它更清楚一點 - 我不是超級關心結果的準確性（事實上，如果它以單精度返回，那很好）。我最關心的是a）讓它在SSE中運行，這樣我可以擺脫目前在兩個其他SSE塊之間存在的不必要的存儲/負載，並且b）儘可能快地運行它。 – BigA 2012-07-22 13:06:57

如果你想要任何一種精度，簡單的算法是非常糟糕的。對於精度範圍縮減算法，參見例如Ng et al., ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit 。

來源

2012-07-20 09:08:05 janneb

我需要再次閱讀這篇文章，以確保我已經理解了它，但該方法不需要比XMM寄存器中更多的位嗎？這種方法在SSE中可行嗎？ – BigA 2012-07-20 11:39:11

@BigA：的確，您無法一次將所有位整合到單個XMM註冊表中。你必須做更復雜的事情。 – janneb 2012-07-20 11:44:12

簡單的「fmod 2pi」算法的一個優點是不容忽視的是，由此引入的舍入誤差通常會抵消角度計算中較早的舍入誤差。例如，如果'Sin'執行自變量減少mod'Math.Pi'，則表達式'Math.Sin（x *（2 * Math.Pi））'的結果通常會更準確，而不是執行減少模π[對於x的某些值，減少模式π會略微好一些，但對於某些情況來說，它會變得更糟]。 – supercat 2014-06-04 15:50:00

對於較大的參數Hayne-Panek algorithm通常使用。然而，Hayne-Panek論文很難閱讀，我建議看看Chapter 11 in the Handbook of Floating-Point Arithmetic以獲得更易於理解的解釋。

來源

2012-07-20 11:00:43

Modulo 2 * Pi使用SSE/SSE2

回答

相關問題