代码之家 › 专栏 › 技术社区 › ThreeStarProgrammer57

双精度运算的快速sse低精度指数

exponential simd sse precision c++

ThreeStarProgrammer57 · 技术社区 · 7 年前

我正在寻找一个快速的sse低精度(~1e-3)指数函数。

我遇到了这个伟大的 answer :

/* max. rel. error = 3.55959567e-2 on [-87.33654, 88.72283] */
__m128 FastExpSse (__m128 x)
{
    __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
    __m128i b = _mm_set1_epi32 (127 * (1 << 23) - 298765);
    __m128i t = _mm_add_epi32 (_mm_cvtps_epi32 (_mm_mul_ps (a, x)), b);
    return _mm_castsi128_ps (t);
}

根据Nicol N.Schraudolph的工作:N.N.Schraudolph。”指数函数的快速紧致逼近〉,《神经计算》,11(4),1999年5月,第853-862页。

现在我需要一个“双精度”版本: __m128d FastExpSSE (__m128d x) . 这是因为我不控制输入和输出精度,这恰好是双精度的,两个转换是双浮点,然后是双浮点,这消耗了50%的CPU资源。

需要做什么改变?

我天真地试过:

__m128i double_to_uint64(__m128d x) {
    x = _mm_add_pd(x, _mm_set1_pd(0x0010000000000000));
    return _mm_xor_si128(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0010000000000000))
    );
}

__m128d FastExpSseDouble(__m128d x) {

    #define S 52
    #define C (1llu << S) / log(2)

    __m128d a = _mm_set1_pd(C); /* (1 << 52) / log(2) */
    __m128i b = _mm_set1_epi64x(127 * (1llu << S) - 298765llu << 29);

    auto y = double_to_uint64(_mm_mul_pd(a, x));

    __m128i t = _mm_add_epi64(y, b);
    return _mm_castsi128_pd(t);
}

当然这是垃圾,因为我不知道我在做什么…

编辑:

关于50%个因素,这是一个非常粗略的估计,比较加速(相对于STD::EXP),将一个单精度数字(大)向量转换为双精度数字列表的加速(不是很好)。

这是我使用的代码:

// gives the result in place
void FastExpSseVector(std::vector<double> & v) { //vector with several millions elements

    const auto I = v.size();

    const auto N = (I / 4) * 4;

    for (int n = 0; n < N; n += 4) {

        float a[4] = { float(v[n]), float(v[n + 1]), float(v[n + 2]), float(v[n + 3]) };

        __m128 x;
        x = _mm_load_ps(a);

        auto r = FastExpSse(x);

        _mm_store_ps(a, r);

        v[n]     = a[0];
        v[n + 1] = a[1];
        v[n + 2] = a[2];
        v[n + 3] = a[3];
    }

    for (int n = N; n < I; ++n) {
        v[n] = FastExp(v[n]);
    }

}

如果我有这个“双精度”版本我会做些什么:

void FastExpSseVectorDouble(std::vector<double> & v) {

    const auto I = v.size();

    const auto N = (I / 2) * 2;

    for (int n = 0; n < N; n += 2) {
        __m128d x;
        x = _mm_load_pd(&v[n]);
        auto r = FastExpSseDouble(x);

        _mm_store_pd(&v[n], r);
    }

    for (int n = N; n < I; ++n) {
        v[n] = FastExp(v[n]);
    }
}

1 回复 | 直到 7 年前

chtz 7 年前

像这样的事情应该能解决问题。你需要调整 1.05 常数以获得较低的最大误差——我懒得这么做:

__m128d fastexp(const __m128d &x)
{
    __m128d scaled = _mm_add_pd(_mm_mul_pd(x, _mm_set1_pd(1.0/std::log(2.0)) ), _mm_set1_pd(3*1024.0-1.05));

    return _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(scaled), 11));
}

这只得到了2.5%的相对精度——为了获得更好的精度,可能需要添加第二个项。

此外,对于溢出或下溢的值,这将导致未指定的值,可以通过夹持 scaled 值到某些值。