Быстрая SSE экспоненциальная низкая точность с использованием операций двойной точности

Я ищу экспоненциальную функцию быстрой SSE с низкой точностью (~ 1e-3).

Я сталкивался с этим великим ответ:

/* 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);
}

По материалам работы Никола Н. Шраудольфа: Н. Н. Шраудольфа. «Быстрое, компактное приближение экспоненциальной функции». Нейронные вычисления, 11 (4), май 1999, с. 853-862.

Теперь мне нужна версия с двойной точностью: __m128d FastExpSSE (__m128d x),
Это потому, что я не контролирую точность ввода и вывода, которая оказывается двойной точностью, и два преобразования double -> float, затем float -> double потребляют 50% ресурсов ЦП.

Какие изменения будут необходимы?

Я наивно пробовал это:

__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]);
}
}

6

Решение

Нечто подобное должно делать работу. Вам нужно настроить 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 значение для некоторых значений.

3

Другие решения

Других решений пока нет …

По вопросам рекламы ammmcru@yandex.ru
Adblock
detector