加速牛頓找到第n個根的方法

Question

讓我用一個陳述來預測這個問題;此代碼按預期工作，但它的速度非常緩慢。有沒有辦法讓它的牛頓方法更快地收斂，或者設置一個__m256 var等於單個float的方法，而不會混淆float數組等等？

__m256 nthRoot(__m256 a, int root){ 

#define aligned __declspec(align(16)) float 

// uses the calculation 
// n_x+1 = (1/root)*(root * x + a/pow(x,root)) 


//initial numbers 
aligned r[8]; 
aligned iN[8]; 
aligned mN[8]; 

//Function I made to fill arrays 
/* 
template<class T> 
void FillArray(T a[],T b) 
{ 
int n = sizeof(a)/sizeof(T); 
for(int i = 0; i < n; a[i++] = b); 
}*/ 

//fills the arrays 
FillArray(iN,(1.0f/(float)root)); 
FillArray(mN,(float)(root-1)); 
FillArray(r,(float)root); 

//loads the arrays into the sse componenets 
__m256 R = _mm256_load_ps(r); 
__m256 Ni = _mm256_load_ps(iN); 
__m256 Nm = _mm256_load_ps(mN); 

    //sets initaial guess to 1/(a * root) 
__m256 x = _mm256_rcp_ps(_mm256_mul_ps(R,a)); 

for(int i = 0; i < 20 ; i ++){ 
    __m256 tmpx = x; 
    for(int k = 0 ; k < root -2 ; k++){ 
     tmpx = _mm256_mul_ps(x,tmpx); 
    } 
    //f over f' 
    __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
    //fmac with Ni*X+tar 
    tar = _mm256_fmadd_ps(Nm,x,tar); 
    //Multipled by Ni 
    x = _mm256_mul_ps(Ni,tar); 
} 
return x; 
} 

編輯＃1

__m256 SSEnthRoot(__m256 a, int root){ 

__m256 R = _mm256_set1_ps((float)root); 
__m256 Ni = _mm256_set1_ps((1.0f)/((float)root)); 
__m256 Nm = _mm256_set1_ps((float)(root -1)); 

__m256 x = _mm256_mul_ps(a,_mm256_rcp_ps(R)); 

for(int i = 0; i < 10 ; i ++){ 
    __m256 tmpx = x; 
    for(int k = 0 ; k < root -2 ; k++){ 
     tmpx = _mm256_mul_ps(x,tmpx); 
    } 
    //f over f' 
    __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
    //mult nm x then add tar because my compiler stoped thinking that fmadd is a valid instruction 
    tar = _mm256_add_ps(_mm256_mul_ps(Nm,x),tar); 
    //Multiplied by the inverse of power 
    x = _mm256_mul_ps(Ni,tar); 
} 

return x; 
} 

任何提示或指針（未存儲器的那種），以使其成爲牛頓方法更快收斂，將不勝感激。

編輯＃2 _mm256_set1_ps刪除（）函數調用_mm256_rcp_ps（），因爲我已經裝了什麼我需要爲R的倒數

__m256 SSEnthRoot(__m256 a, int root){ 
__m256 R = _mm256_set1_ps((float)root); 
__m256 Ni = _mm256_rcp_ps(R); 
__m256 Nm = _mm256_set1_ps((float)(root -1)); 

__m256 x = _mm256_mul_ps(a,Ni); 
for(int i = 0; i < 20 ; i ++){ 
    __m256 tmpx = x; 
    for(int k = 0 ; k < root -2 ; k++) 
     tmpx = _mm256_mul_ps(x,tmpx); 
    //f over f' 
    __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
    //fmac with Ni*X+tar 
      //my compiler believes in fmac again 
    tar = _mm256_fmadd_ps(Nm,x,tar); 
    //Multiplied by the inverse of power 
    x = _mm256_mul_ps(Ni,tar); 
} 
return x; 
} 

編輯＃3

__m256 SSEnthRoot(__m256 a, int root){ 
__m256 Ni = _mm256_set1_ps(1.0f/(float)root); 
__m256 Nm = _mm256_set1_ps((float)(root -1)); 
__m256 x = _mm256_mul_ps(a,Ni); 
for(int i = 0; i < 20 ; i ++){ 
    __m256 tmpx = x; 
    for(int k = 0 ; k < root -2 ; k++) 
     tmpx = _mm256_mul_ps(x,tmpx); 
    __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
    tar = _mm256_fmadd_ps(Nm,x,tar); 
    x = _mm256_mul_ps(Ni,tar); 
} 
return x; 
} 

Answer 1

您的pow功能效率低下。

for(int k = 0 ; k < root -2 ; k++) 
     tmpx = _mm256_mul_ps(x,tmpx); 

在你的例子中，你正在採取第29根。您需要pow(x, 29-1) = x^28。目前你使用27次乘法，但只用6次乘法就可以做到這一點。

x^28 = (x^4)*(x^8)*(x^16) 
x^4 = y -> 2 multiplications 
x^8 = y*y = z -> 1 multiplication 
x^16 = z^2 = w-> 1 multiplications 
y*z*w -> 2 multiplications 
6 multiplications in total 

這裏是你代碼的改進版本，這是大約快一倍我的系統上。它使用了我創建的新功能pow_avx_fast，它使用AVX一次爲8個浮動塊生成x^n。它確實例如x^28乘以6次而不是27. 請參閱我的答案。我發現一個版本可以在一定的容差範圍內找到結果xacc。如果收斂發生得很快，這可能會更快。

inline __m256 pow_avx_fast(__m256 x, const int n) { 
    //n must be greater than zero 
    if(n%2 == 0) { 
     return pow_avx_fast(_mm256_mul_ps(x, x), n/2); 
    } 
    else { 
     if(n>1) return _mm256_mul_ps(x,pow_avx_fast(_mm256_mul_ps(x, x), (n-1)/2)); 
     return x; 
    } 
} 

inline __m256 SSEnthRoot_fast(__m256 a, int root) { 
    // n_x+1 = (1/root)*((root-1) * x + a/pow(x,root-1)) 
    __m256 R = _mm256_set1_ps((float)root); 
    __m256 Ni = _mm256_rcp_ps(R); 
    __m256 Nm = _mm256_set1_ps((float)(root -1)); 

    __m256 x = _mm256_mul_ps(a,Ni); 
    for(int i = 0; i < 20 ; i ++) { 
     __m256 tmpx = pow_avx_fast(x, root-1); 
     //f over f' 
     __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
     //fmac with Ni*X+tar 
     //tar = _mm256_fmadd_ps(Nm,x,tar); 
     tar = _mm256_add_ps(_mm256_mul_ps(Nm,x),tar); 
     //Multiplied by the inverse of power 
     x = _mm256_mul_ps(Ni,tar); 
    } 
    return x; 
} 

欲瞭解更多信息，如何編寫一個高效pow功能看到這些鏈接http://en.wikipedia.org/wiki/Addition-chain_exponentiation和 http://en.wikipedia.org/wiki/Exponentiation_by_squaring

而且，你的初始猜測可能不是那麼好。這裏是標量代碼來找到基於你的方法的第n根（但使用數學可能比你更快的函數pow）。它需要大約50次迭代來解決16的第4根（這是2）。對於你使用它的20次迭代，它會返回超過4000，而在接近2.0的地方不會。所以你需要調整你的方法來做足夠的迭代，以確保在一定的容差範圍內有合理的答案。

float fx(float a, int n, float x) { 
    return 1.0f/n * ((n-1)*x + a/pow(x, n-1)); 
} 
float scalar_nthRoot_v2(float a, int root) { 
    //sets initaial guess to 1/(a * root) 
    float x = 1.0f/(a*root); 
    printf("x0 %f\n", x); 
    for(int i = 0; i<50; i++) { 
     x = fx(a, root, x); 
     printf("x %f\n", x); 
    } 
    return x; 
} 

我從這裏得到了牛頓法的公式。 http://en.wikipedia.org/wiki/Nth_root_algorithm

這裏是一個版本的功能，這給了一定的公差範圍內xacc的結果或退出後，如果nmax迭代不收斂。如果收斂發生在少於20次的迭代中，此函數可能比您的方法快得多。它要求所有八個浮體一次會聚。換句話說，如果七個會聚，一個不會，那麼其他七個必須等待不會收斂的那個。這就是SIMD的問題（在GPU上），但總的來說，它比沒有SIMD的情況下還要快。

int get_mask(const __m256 dx, const float xacc) { 
    __m256i mask = _mm256_castps_si256(_mm256_cmp_ps(dx, _mm256_set1_ps(xacc), _CMP_GT_OQ)); 
    return _mm_movemask_epi8(_mm256_castsi256_si128(mask)) + _mm_movemask_epi8(_mm256_extractf128_si256(mask,1)); 
} 

inline __m256 SSEnthRoot_fast_xacc(const __m256 a, const int root, const int nmax, float xacc) { 
    // n_x+1 = (1/root)*(root * x + a/pow(x,root)) 
    __m256 R = _mm256_set1_ps((float)root); 
    __m256 Ni = _mm256_rcp_ps(R); 
    //__m256 Ni = _mm256_set1_ps(1.0f/root); 
    __m256 Nm = _mm256_set1_ps((float)(root -1)); 

    __m256 x = _mm256_mul_ps(a,Ni); 

    for(int i = 0; i <nmax ; i ++) { 
     __m256 tmpx = pow_avx_fast(x, root-1); 
     __m256 tar = _mm256_mul_ps(a,_mm256_rcp_ps(tmpx)); 
     //tar = _mm256_fmadd_ps(Nm,x,tar); 
     tar = _mm256_add_ps(_mm256_mul_ps(Nm,x),tar); 
     tmpx = _mm256_mul_ps(Ni,tar); 
     __m256 dx = _mm256_sub_ps(tmpx,x); 
     dx = _mm256_max_ps(_mm256_sub_ps(_mm256_setzero_ps(), dx), dx); //fabs(dx) 
     int cnt = get_mask(dx, xacc); 
     if(cnt == 0) return x; 
     x = tmpx; 
    } 
    return x; //at least one value out of eight did not converge by nmax. 
} 

下面是avx的pow函數的更一般的版本，它也適用於n < = 0。

__m256 pow_avx(__m256 x, const int n) { 
    if(n<0) { 
     return pow_avx(_mm256_rcp_ps(x), -n); 
    } 
    else if(n == 0) { 
     return _mm256_set1_ps(1.0f); 
    } 
    else if(n == 1) { 
     return x; 
    } 
    else if(n%2 ==0) { 
     return pow_avx(_mm256_mul_ps(x, x), n/2); 
    } 
    else { 
     return _mm256_mul_ps(x,pow_avx(_mm256_mul_ps(x, x), (n-1)/2)); 
    } 
} 

一些其他建議

您可以使用它找到的第n根SIMD數學庫。 SIMD math libraries for SSE and AVX

對於英特爾，您可以使用昂貴且封閉源代碼的SVML（英特爾的OpenCL驅動程序使用SVML，因此您可以免費獲得）。對於AMD，您可以使用免費但封閉源代碼的LIBM。有幾個開源SIMD數學庫，如http://software-lisc.fbk.eu/avx_mathfun/和https://bitbucket.org/eschnett/vecmathlib/wiki/Home

Answer 2

要將__m256載體的所有元素設置爲單個值：

__m256 v = _mm256_set1_ps(1.0f); 

或在您的具體情況：

__m256 R = _mm256_set1_ps((float)root); 
__m256 Ni = _mm256_set1_ps((1.0f/(float)root)); 
__m256 Nm = _mm256_set1_ps((float)(root-1)); 

一旦你做出這個改變很明顯，你可以擺脫的FillArray東西。

Answer 3

也許你應該在日誌域中這樣做。

pow(a,1/root) == exp(log(x) /root) 

朱利安POMMIER有sse_mathfun.h有SSE，SSE2日誌和記錄功能，但我不能說我用那些特別。這些技術可以擴展到avx。