0

次のリンクhttp://rosettacode.org/wiki/Nth_root#Cから double 値の C で n 番目のルートを検索する関数を変換して、 AVX を使用して一度に 8 つの float の n 番目のルートを検索しようとしています。

そのコードの一部は DBL_EPSILON * 10 を使用しています。ただし、これを変換して AVX で float を使用する場合、FLT_EPSILON*1000 を使用する必要があります。そうしないと、コードがハングして収束しません。FLT_EPSILON を印刷すると、オーダー 1E-7 であることがわかります。しかし、このリンクhttp://www.cplusplus.com/reference/cfloat/ には、1E-5 である必要があると書かれています。DBL_EPSILON を印刷すると 1E-16 ですが、リンクには 1E-9 のみである必要があると記載されています。どうしたの？

これまでのコードは次のとおりです (完全には最適化されていません)。

#include <stdio.h>
#include <float.h>
#include <immintrin.h>                 // AVX

inline double abs_(double x) { 
    return x >= 0 ? x : -x; 
}

double pow_(double x, int e)
{
    double ret = 1;
    for (ret = 1; e; x *= x, e >>= 1) {
        if ((e & 1)) ret *= x;
    }
    return ret;
}

double root(double a, int n)
{
    double d, x = 1;
    x = a/n;
    if (!a) return 0;
    //if (n < 1 || (a < 0 && !(n&1))) return 0./0.; /* NaN */

    int cnt = 0;
    do {
        cnt++;  
        d = (a / pow_(x, n - 1) - x) / n;
        x+= d;
    } while (abs_(d) >= abs_(x) * (DBL_EPSILON * 10));
    printf("%d\n", cnt);    

    return x;
}


__m256 pow_avx(__m256 x, int e) {
    __m256 ret = _mm256_set1_ps(1.0f);
    for (; e; x = _mm256_mul_ps(x,x), e >>= 1) {
        if ((e & 1)) ret = _mm256_mul_ps(x,ret);
    }
    return ret;
}

inline __m256 abs_avx (__m256 x) { 
    return _mm256_max_ps(_mm256_sub_ps(_mm256_setzero_ps(), x), x);
    //return x >= 0 ? x : -x; 
}

int get_mask(const __m256 d, const __m256 x) {
    __m256 ad = abs_avx(d);
    __m256 ax = abs_avx(x);
    __m256i mask = _mm256_castps_si256(_mm256_cmp_ps(ad, ax, _CMP_GT_OQ));
    return _mm_movemask_epi8(_mm256_castsi256_si128(mask)) + _mm_movemask_epi8(_mm256_extractf128_si256(mask,1));
}

__m256 root_avx(__m256 a, int n) {
    printf("%e\n", FLT_EPSILON);
    printf("%e\n", DBL_EPSILON);
    printf("%e\n", FLT_EPSILON*1000.0f);
    __m256 d;
    __m256 x = _mm256_set1_ps(1.0f);
    //if (!a) return 0;
    //if (n < 1 || (a < 0 && !(n&1))) return 0./0.; /* NaN */
    __m256 in = _mm256_set1_ps(1.0f/n);
    __m256 xtmp;
    do {
        d = _mm256_rcp_ps(pow_avx(x, n - 1));
        d = _mm256_sub_ps(_mm256_mul_ps(a,d),x);
        d = _mm256_mul_ps(d,in);
        //d = (a / pow_avx(x, n - 1) - x) / n;
        x = _mm256_add_ps(x, d); //x+= d;
        xtmp =_mm256_mul_ps(x, _mm256_set1_ps(FLT_EPSILON*100.0f));
    //} while (abs_(d) >= abs_(x) * (DBL_EPSILON * 10));
    } while (get_mask(d, xtmp)); 
    return x;
}

int main()
{
    __m256 d = _mm256_set1_ps(16.0f);
    __m256 out = root_avx(d, 4);
    float result[8];
    int i;
    _mm256_storeu_ps(result, out);

    for(i=0; i<8; i++) {
        printf("%f\n", result[i]);
    } printf("\n");

    //double x = 16;
    //printf("root(%g, 15) = %g\n", x, root(x, 4));

    //double x = pow_(-3.14159, 15);
    //printf("root(%g, 15) = %g\n", x, root(x, 15));
    return 0;
}
4

2 に答える 2

3

_mm256_rcp_ps命令にマップされる はrcpps、おおよその逆数のみを実行します。Intel 64 and IA-32 Architectures Software Developer's Manualによると、その相対誤差は最大 1.5•2 -12である可能性があります。これは、ルート ファインダを正確に収束させるには不十分100*FLT_EPSILONです。

次のような正確な除算を使用できます。

d = pow_avx(x, n-1);
d = _mm256_sub_ps(_mm256_div_ps(a, d), x);

または、逆数の推定にいくつかの調整ステップを追加します。

ちなみに、コンパイラが SIMD オブジェクトでの通常の C 演算子の使用をサポートしている場合は、代わりに通常の C 演算子の使用を検討してください。

d = pow_avx(x, n-1);
d = a/d - x;
于 2013-06-14T16:38:16.030 に答える
1

1e-5は、単に C 標準が実装で使用できる最大値ですFLT_EPSILON。実際には、2 -23のイプシロンを持つ IEEE-754 単精度を使用します。これは約1e-7です。

于 2013-06-14T14:22:31.407 に答える