Быстрая реализация экспоненциальной функции с использованием AVX

Question

Быстрая реализация экспоненциальной функции с использованием AVX

Я ищу эффективное (быстрое) приближение экспоненциальной функции, работающей с элементами AVX (плавающая точка одинарной точности). А именно - __m256 _mm256_exp_ps( __m256 x ) без SVML.

Относительная точность должна быть примерно такой: ~1e-6 или ~20 битов мантиссы (1 часть в 2^20).

Я был бы счастлив, если бы он был написан в стиле C с использованием встроенных функций Intel.
Код должен быть переносимым (Windows, macOS, Linux, MSVC, ICC, GCC и т. Д.).

Это похоже на быструю реализацию экспоненциальной функции с использованием SSE, но этот вопрос ищет очень быстро с низкой точностью (текущий ответ дает точность около 1e-3).

Также этот вопрос ищет AVX / AVX2 (и FMA). Но учтите, что ответы на оба вопроса легко переносятся между SSE4 __m128 или AVX2 __m256, поэтому будущие читатели должны выбирать, основываясь на требуемой компромиссе между точностью и производительностью.

14

x86 simd avx exponential avx2

Источник

user195787 19 фев '18 в 10:08

5 ответов

Решение

Так как быстрое вычисление exp() требует манипулирования полем экспоненты операндов с плавающей точкой IEEE-754, AVX не очень подходит для этого вычисления, так как в нем отсутствуют целочисленные операции. Поэтому я сосредоточусь на AVX2, Поддержка слияния-умножения - технически особенность, отличная от AVX2 поэтому я предоставляю два кодовых пути, с использованием и без использования FMA, контролируемых макросом USE_FMA,

Код ниже вычисляет exp() почти до желаемой точности 10 ^-6. Использование FMA здесь не дает каких-либо существенных улучшений, но оно должно обеспечивать преимущество в производительности на платформах, которые его поддерживают.

Алгоритм, использованный в предыдущем ответе для реализации SSE с более низкой точностью, не является полностью расширяемым для довольно точной реализации, поскольку он содержит некоторые вычисления с плохими числовыми свойствами, которые, однако, не имеют значения в этом контексте. Вместо вычисления e ^x = 2 ⁱ * 2 ^f, с f в [0,1] или f в [-½, ½] выгодно вычислять e ^x = 2 ⁱ * e ^f с f в более узком интервале [-½log 2, ½log 2], где log обозначает натуральный логарифм.

Для этого сначала вычислим i = rint(x * log2(e)), затем f = x - log(2) * i, Важно, что в последнем вычислении должна использоваться точность, превышающая нативную, для предоставления точного сокращенного аргумента, который должен быть передан в основное приближение. Для этого мы используем схему Коди-Уэйта, впервые опубликованную в WJ Cody & W. Waite, "Руководство по программному обеспечению для элементарных функций", Prentice Hall 1980. Константа log(2) разбивается на "большую" часть большего величина и "низкая" часть гораздо меньшей величины, которая содержит разницу между "высокой" частью и математической константой.

Верхняя часть выбирается с достаточным количеством завершающих нулевых битов в мантиссе, так что произведение i с "высокой" частью точно представимо в родной точности. Здесь я выбрал "старшую" часть с восемью завершающими нулевыми битами, как i наверняка уместится в восемь бит.

По сути, мы вычисляем f = x - i * log(2) _high - i * log(2) _low. Этот приведенный аргумент передается в базовое приближение, которое является полиномиальным минимаксным приближением, а результат масштабируется на 2 ^i, как и в предыдущем ответе.

#include <immintrin.h>

#define USE_FMA 0

/* compute exp(x) for x in [-87.33654f, 88.72283] 
   maximum relative error: 3.1575e-6 (USE_FMA = 0); 3.1533e-6 (USE_FMA = 1)
*/
__m256 faster_more_accurate_exp_avx2 (__m256 x)
{
    __m256 t, f, p, r;
    __m256i i, j;

    const __m256 l2e = _mm256_set1_ps (1.442695041f); /* log2(e) */
    const __m256 l2h = _mm256_set1_ps (-6.93145752e-1f); /* -log(2)_hi */
    const __m256 l2l = _mm256_set1_ps (-1.42860677e-6f); /* -log(2)_lo */
    /* coefficients for core approximation to exp() in [-log(2)/2, log(2)/2] */
    const __m256 c0 =  _mm256_set1_ps (0.041944388f);
    const __m256 c1 =  _mm256_set1_ps (0.168006673f);
    const __m256 c2 =  _mm256_set1_ps (0.499999940f);
    const __m256 c3 =  _mm256_set1_ps (0.999956906f);
    const __m256 c4 =  _mm256_set1_ps (0.999999642f);

    /* exp(x) = 2^i * e^f; i = rint (log2(e) * x), f = x - log(2) * i */
    t = _mm256_mul_ps (x, l2e);      /* t = log2(e) * x */
    r = _mm256_round_ps (t, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC); /* r = rint (t) */

#if USE_FMA
    f = _mm256_fmadd_ps (r, l2h, x); /* x - log(2)_hi * r */
    f = _mm256_fmadd_ps (r, l2l, f); /* f = x - log(2)_hi * r - log(2)_lo * r */
#else // USE_FMA
    p = _mm256_mul_ps (r, l2h);      /* log(2)_hi * r */
    f = _mm256_add_ps (x, p);        /* x - log(2)_hi * r */
    p = _mm256_mul_ps (r, l2l);      /* log(2)_lo * r */
    f = _mm256_add_ps (f, p);        /* f = x - log(2)_hi * r - log(2)_lo * r */
#endif // USE_FMA

    i = _mm256_cvtps_epi32(t);       /* i = (int)rint(t) */

    /* p ~= exp (f), -log(2)/2 <= f <= log(2)/2 */
    p = c0;                          /* c0 */
#if USE_FMA
    p = _mm256_fmadd_ps (p, f, c1);  /* c0*f+c1 */
    p = _mm256_fmadd_ps (p, f, c2);  /* (c0*f+c1)*f+c2 */
    p = _mm256_fmadd_ps (p, f, c3);  /* ((c0*f+c1)*f+c2)*f+c3 */
    p = _mm256_fmadd_ps (p, f, c4);  /* (((c0*f+c1)*f+c2)*f+c3)*f+c4 ~= exp(f) */
#else // USE_FMA
    p = _mm256_mul_ps (p, f);        /* c0*f */
    p = _mm256_add_ps (p, c1);       /* c0*f+c1 */
    p = _mm256_mul_ps (p, f);        /* (c0*f+c1)*f */
    p = _mm256_add_ps (p, c2);       /* (c0*f+c1)*f+c2 */
    p = _mm256_mul_ps (p, f);        /* ((c0*f+c1)*f+c2)*f */
    p = _mm256_add_ps (p, c3);       /* ((c0*f+c1)*f+c2)*f+c3 */
    p = _mm256_mul_ps (p, f);        /* (((c0*f+c1)*f+c2)*f+c3)*f */
    p = _mm256_add_ps (p, c4);       /* (((c0*f+c1)*f+c2)*f+c3)*f+c4 ~= exp(f) */
#endif // USE_FMA

    /* exp(x) = 2^i * p */
    j = _mm256_slli_epi32 (i, 23); /* i << 23 */
    r = _mm256_castsi256_ps (_mm256_add_epi32 (j, _mm256_castps_si256 (p))); /* r = p * 2^i */

    return r;
}

Если требуется более высокая точность, степень полиномиальной аппроксимации может быть увеличена на единицу с использованием следующего набора коэффициентов:

/* maximum relative error: 1.7428e-7 (USE_FMA = 0); 1.6586e-7 (USE_FMA = 1) */
const __m256 c0 =  _mm256_set1_ps (0.008301110f);
const __m256 c1 =  _mm256_set1_ps (0.041906696f);
const __m256 c2 =  _mm256_set1_ps (0.166674897f);
const __m256 c3 =  _mm256_set1_ps (0.499990642f);
const __m256 c4 =  _mm256_set1_ps (0.999999762f);
const __m256 c5 =  _mm256_set1_ps (1.000000000f);

10

Источник

user780717 03 мар '18 в 23:56

Я много играл с этим и обнаружил этот, который имеет относительную точность около ~1-07e и простой для преобразования в векторные инструкции. Имея всего 4 константы, 5 умножений и 1 деление, это в два раза быстрее, чем встроенная функция exp().

float fast_exp(float x)
{
    const float c1 = 0.007972914726F;
    const float c2 = 0.1385283768F;
    const float c3 = 2.885390043F;
    const float c4 = 1.442695022F;      
    x *= c4; //convert to 2^(x)
    int intPart = (int)x;
    x -= intPart;
    float xx = x * x;
    float a = x + c1 * xx * x;
    float b = c3 + c2 * xx;
    float res = (b + a) / (b - a);
    reinterpret_cast<int &>(res) += intPart << 23; // res *= 2^(intPart)
    return res;
}

Конвертация в AVX (обновлено)

__m256 _mm256_exp_ps(__m256 _x)
{
    __m256 c1 = _mm256_set1_ps(0.007972914726F);
    __m256 c2 = _mm256_set1_ps(0.1385283768F);
    __m256 c3 = _mm256_set1_ps(2.885390043F);
    __m256 c4 = _mm256_set1_ps(1.442695022F);
    __m256 x = _mm256_mul_ps(_x, c4); //convert to 2^(x)
    __m256 intPartf = _mm256_round_ps(x, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
    x = _mm256_sub_ps(x, intPartf);
    __m256 xx = _mm256_mul_ps(x, x);
    __m256 a = _mm256_add_ps(x, _mm256_mul_ps(c1, _mm256_mul_ps(xx, x))); //can be improved with FMA
    __m256 b = _mm256_add_ps(c3, _mm256_mul_ps(c2, xx));
    __m256 res = _mm256_div_ps(_mm256_add_ps(b, a), _mm256_sub_ps(b, a));
    __m256i intPart = _mm256_cvtps_epi32(intPartf); //res = 2^intPart. Can be improved with AVX2!
    __m128i ii0 = _mm_slli_epi32(_mm256_castsi256_si128(intPart), 23);
    __m128i ii1 = _mm_slli_epi32(_mm256_extractf128_si256(intPart, 1), 23);     
    __m128i res_0 = _mm_add_epi32(ii0, _mm256_castsi256_si128(_mm256_castps_si256(res)));
    __m128i res_1 = _mm_add_epi32(ii1, _mm256_extractf128_si256(_mm256_castps_si256(res), 1));
    return _mm256_insertf128_ps(_mm256_castsi256_ps(_mm256_castsi128_si256(res_0)), _mm_castsi128_ps(res_1), 1);
}

6

Источник

user2604941 08 фев '20 в 00:09

Вы можете приблизить экспоненту с помощью ряда Тейлора:

exp(z) = 1 + z + pow(z,2)/2 + pow(z,3)/6 + pow(z,4)/24 + ...

Для этого вам понадобятся только операции сложения и умножения из AVX. Коэффициенты, такие как 1/2, 1/6, 1/24 и т. Д. Быстрее, если жестко закодированы, а затем умножены на, а не разделены.

Возьмите столько членов последовательности, сколько требует ваша точность. Обратите внимание, что вы получите относительную ошибку: для малых z это может быть 1e-6 в абсолюте, но для больших z это будет больше чем 1e-6 в абсолюте еще abs(E-E1)/abs(E) - 1 меньше чем 1e-6 (где E точный показатель и E1 это то, что вы получаете с приближением).

ОБНОВЛЕНИЕ: Как @Peter Cordes упомянул в комментарии, точность может быть улучшена путем разделения возведения в степень целочисленной и дробной частей, обработки целочисленной части путем манипулирования полем экспоненты двоичного float представление (которое основано на 2^ х, а не е ^ х). Тогда ваша серия Тейлора должна только минимизировать ошибку на небольшом диапазоне.

0

Источник

user1915854 19 фев '18 в 11:34

Для нормализованных входных данных ([-1,1]) вы можете использовать полиномиальную аппроксимацию:

          // compute Simd exp() at a time (only optimized for Type=float)
    template<typename Type, int Simd>
    inline
    void expFast(float * const __restrict__ data, float * const __restrict__ result) noexcept
    {

        alignas(64)
        Type resultData[Simd];

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    Type(0.0001972591916103993980868836)*data[i] + Type(0.001433947376170863208244555);
        }

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    resultData[i]*data[i] + Type(0.008338950118885968265658448);
        }

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    resultData[i]*data[i] + Type(0.04164162895364054151059463);
        }

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    resultData[i]*data[i] + Type(0.1666645212581130408580066);
        }

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    resultData[i]*data[i] + Type(0.5000045184212300597437206);
        }

        
        for(int i=0;i<Simd;i++)
        {
            resultData[i] =    resultData[i]*data[i] + Type(0.9999999756072401879691824);
        }

        
        for(int i=0;i<Simd;i++)
        {
            result[i] =    resultData[i]*data[i] + Type(0.999999818912344906607359);
        }

    }

Он имеет среднюю ошибку 0,5 ULPS и максимальную ошибку 10 ULPS для 64 миллионов точек, выбранных между -1 и 1. Он имеет 10-кратное ускорение по сравнению со std::exp на AVX1 (бульдозер).

Я думаю, вы можете комбинировать эту функцию с целыми умножениями, чтобы поддерживать все степени. Но простая часть умножения требует O(logN) вместо O(N), чтобы быть достаточно быстрой для больших мощностей. Например, если вы вычислили x^10, то потребуется всего 1 дополнительная операция с самим собой, чтобы получить x^20 вместо 10 дополнительных операций путем умножения на x.

При использовании в цикле компилятор генерирует это:

      .L2:
    vmovaps zmm1, ZMMWORD PTR [rax]
    add     rax, 64
    vmovaps zmm0, zmm1
    vfmadd132ps     zmm0, zmm8, zmm9
    vfmadd132ps     zmm0, zmm7, zmm1
    vfmadd132ps     zmm0, zmm6, zmm1
    vfmadd132ps     zmm0, zmm5, zmm1
    vfmadd132ps     zmm0, zmm4, zmm1
    vfmadd132ps     zmm0, zmm3, zmm1
    vfmadd132ps     zmm0, zmm2, zmm1
    vmovaps ZMMWORD PTR [rax-64], zmm0
    cmp     rax, rdx
    jne     .L2

Я думаю, что это достаточно быстро, чтобы сэкономить несколько циклов для обработки целочисленных степеней ввода, возможно, до предела числа с плавающей запятой (10^38).

0

Источник

user1470092 03 май '22 в 16:55

Другие вопросы по тегам x86 simd avx exponential avx2

user2439725 19 фев '18 в 15:27 2018-02-19 15:27 · Accepted Answer · 2018-02-19 15:27

exp функция из avx_mathfun использует уменьшение диапазона в сочетании с полиномом, подобным аппроксимации Чебышева, для вычисления 8 exp -s параллельно с инструкциями AVX. Используйте правильные настройки компилятора, чтобы убедиться, что addps а также mulps где возможно, слиты с инструкциями FMA.

Это довольно просто адаптировать оригинал exp код из avx_mathfun в переносимый (через разные компиляторы) встроенный код C / AVX2. Оригинальный код использует атрибуты выравнивания стиля gcc и оригинальные макросы. Модифицированный код, который использует стандарт _mm256_set1_ps() вместо этого ниже небольшой тестовый код и таблица. Модифицированный код требует AVX2.

Следующий код используется для простого теста:

int main(){
    int i;
    float xv[8];
    float yv[8];
    __m256 x = _mm256_setr_ps(1.0f, 2.0f, 3.0f ,4.0f ,5.0f, 6.0f, 7.0f, 8.0f);
    __m256 y = exp256_ps(x);
    _mm256_store_ps(xv,x);
    _mm256_store_ps(yv,y);

    for (i=0;i<8;i++){
        printf("i = %i, x = %e, y = %e \n",i,xv[i],yv[i]);
    }
    return 0;
}

Вывод вроде бы в порядке:

i = 0, x = 1.000000e+00, y = 2.718282e+00 
i = 1, x = 2.000000e+00, y = 7.389056e+00 
i = 2, x = 3.000000e+00, y = 2.008554e+01 
i = 3, x = 4.000000e+00, y = 5.459815e+01 
i = 4, x = 5.000000e+00, y = 1.484132e+02 
i = 5, x = 6.000000e+00, y = 4.034288e+02 
i = 6, x = 7.000000e+00, y = 1.096633e+03 
i = 7, x = 8.000000e+00, y = 2.980958e+03

Модифицированный код (AVX2):

#include <stdio.h>
#include <immintrin.h>
/*     gcc -O3 -m64 -Wall -mavx2 -march=broadwell  expc.c    */

__m256 exp256_ps(__m256 x) {
/* Modified code. The original code is here: https://github.com/reyoung/avx_mathfun

   AVX implementation of exp
   Based on "sse_mathfun.h", by Julien Pommier
   http://gruntthepeon.free.fr/ssemath/
   Copyright (C) 2012 Giovanni Garberoglio
   Interdisciplinary Laboratory for Computational Science (LISC)
   Fondazione Bruno Kessler and University of Trento
   via Sommarive, 18
   I-38123 Trento (Italy)
  This software is provided 'as-is', without any express or implied
  warranty.  In no event will the authors be held liable for any damages
  arising from the use of this software.
  Permission is granted to anyone to use this software for any purpose,
  including commercial applications, and to alter it and redistribute it
  freely, subject to the following restrictions:
  1. The origin of this software must not be misrepresented; you must not
     claim that you wrote the original software. If you use this software
     in a product, an acknowledgment in the product documentation would be
     appreciated but is not required.
  2. Altered source versions must be plainly marked as such, and must not be
     misrepresented as being the original software.
  3. This notice may not be removed or altered from any source distribution.
  (this is the zlib license)
*/
/* 
  To increase the compatibility across different compilers the original code is
  converted to plain AVX2 intrinsics code without ingenious macro's,
  gcc style alignment attributes etc. The modified code requires AVX2
*/
__m256   exp_hi        = _mm256_set1_ps(88.3762626647949f);
__m256   exp_lo        = _mm256_set1_ps(-88.3762626647949f);

__m256   cephes_LOG2EF = _mm256_set1_ps(1.44269504088896341);
__m256   cephes_exp_C1 = _mm256_set1_ps(0.693359375);
__m256   cephes_exp_C2 = _mm256_set1_ps(-2.12194440e-4);

__m256   cephes_exp_p0 = _mm256_set1_ps(1.9875691500E-4);
__m256   cephes_exp_p1 = _mm256_set1_ps(1.3981999507E-3);
__m256   cephes_exp_p2 = _mm256_set1_ps(8.3334519073E-3);
__m256   cephes_exp_p3 = _mm256_set1_ps(4.1665795894E-2);
__m256   cephes_exp_p4 = _mm256_set1_ps(1.6666665459E-1);
__m256   cephes_exp_p5 = _mm256_set1_ps(5.0000001201E-1);
__m256   tmp           = _mm256_setzero_ps(), fx;
__m256i  imm0;
__m256   one           = _mm256_set1_ps(1.0f);

        x     = _mm256_min_ps(x, exp_hi);
        x     = _mm256_max_ps(x, exp_lo);

  /* express exp(x) as exp(g + n*log(2)) */
        fx    = _mm256_mul_ps(x, cephes_LOG2EF);
        fx    = _mm256_add_ps(fx, _mm256_set1_ps(0.5f));
        tmp   = _mm256_floor_ps(fx);
__m256  mask  = _mm256_cmp_ps(tmp, fx, _CMP_GT_OS);    
        mask  = _mm256_and_ps(mask, one);
        fx    = _mm256_sub_ps(tmp, mask);
        tmp   = _mm256_mul_ps(fx, cephes_exp_C1);
__m256  z     = _mm256_mul_ps(fx, cephes_exp_C2);
        x     = _mm256_sub_ps(x, tmp);
        x     = _mm256_sub_ps(x, z);
        z     = _mm256_mul_ps(x,x);

__m256  y     = cephes_exp_p0;
        y     = _mm256_mul_ps(y, x);
        y     = _mm256_add_ps(y, cephes_exp_p1);
        y     = _mm256_mul_ps(y, x);
        y     = _mm256_add_ps(y, cephes_exp_p2);
        y     = _mm256_mul_ps(y, x);
        y     = _mm256_add_ps(y, cephes_exp_p3);
        y     = _mm256_mul_ps(y, x);
        y     = _mm256_add_ps(y, cephes_exp_p4);
        y     = _mm256_mul_ps(y, x);
        y     = _mm256_add_ps(y, cephes_exp_p5);
        y     = _mm256_mul_ps(y, z);
        y     = _mm256_add_ps(y, x);
        y     = _mm256_add_ps(y, one);

  /* build 2^n */
        imm0  = _mm256_cvttps_epi32(fx);
        imm0  = _mm256_add_epi32(imm0, _mm256_set1_epi32(0x7f));
        imm0  = _mm256_slli_epi32(imm0, 23);
__m256  pow2n = _mm256_castsi256_ps(imm0);
        y     = _mm256_mul_ps(y, pow2n);
        return y;
}

int main(){
    int i;
    float xv[8];
    float yv[8];
    __m256 x = _mm256_setr_ps(1.0f, 2.0f, 3.0f ,4.0f ,5.0f, 6.0f, 7.0f, 8.0f);
    __m256 y = exp256_ps(x);
    _mm256_store_ps(xv,x);
    _mm256_store_ps(yv,y);

    for (i=0;i<8;i++){
        printf("i = %i, x = %e, y = %e \n",i,xv[i],yv[i]);
    }
    return 0;
}

Как отмечает @Peter Cordes, должна быть возможность заменить _mm256_floor_ps(fx + 0.5f) от _mm256_round_ps(fx), Кроме того, mask = _mm256_cmp_ps(tmp, fx, _CMP_GT_OS); и следующие две строки кажутся излишними. Дальнейшая оптимизация возможна путем объединения cephes_exp_C1 а также cephes_exp_C2 в inv_LOG2EF, Это приводит к следующему коду, который не был тщательно протестирован!

#include <stdio.h>
#include <immintrin.h>
#include <math.h>
/*    gcc -O3 -m64 -Wall -mavx2 -march=broadwell  expc.c -lm     */

__m256 exp256_ps(__m256 x) {
/* Modified code from this source: https://github.com/reyoung/avx_mathfun

   AVX implementation of exp
   Based on "sse_mathfun.h", by Julien Pommier
   http://gruntthepeon.free.fr/ssemath/
   Copyright (C) 2012 Giovanni Garberoglio
   Interdisciplinary Laboratory for Computational Science (LISC)
   Fondazione Bruno Kessler and University of Trento
   via Sommarive, 18
   I-38123 Trento (Italy)
  This software is provided 'as-is', without any express or implied
  warranty.  In no event will the authors be held liable for any damages
  arising from the use of this software.
  Permission is granted to anyone to use this software for any purpose,
  including commercial applications, and to alter it and redistribute it
  freely, subject to the following restrictions:
  1. The origin of this software must not be misrepresented; you must not
     claim that you wrote the original software. If you use this software
     in a product, an acknowledgment in the product documentation would be
     appreciated but is not required.
  2. Altered source versions must be plainly marked as such, and must not be
     misrepresented as being the original software.
  3. This notice may not be removed or altered from any source distribution.
  (this is the zlib license)

*/
/* 
  To increase the compatibility across different compilers the original code is
  converted to plain AVX2 intrinsics code without ingenious macro's,
  gcc style alignment attributes etc.
  Moreover, the part "express exp(x) as exp(g + n*log(2))" has been significantly simplified.
  This modified code is not thoroughly tested!
*/


__m256   exp_hi        = _mm256_set1_ps(88.3762626647949f);
__m256   exp_lo        = _mm256_set1_ps(-88.3762626647949f);

__m256   cephes_LOG2EF = _mm256_set1_ps(1.44269504088896341f);
__m256   inv_LOG2EF    = _mm256_set1_ps(0.693147180559945f);

__m256   cephes_exp_p0 = _mm256_set1_ps(1.9875691500E-4);
__m256   cephes_exp_p1 = _mm256_set1_ps(1.3981999507E-3);
__m256   cephes_exp_p2 = _mm256_set1_ps(8.3334519073E-3);
__m256   cephes_exp_p3 = _mm256_set1_ps(4.1665795894E-2);
__m256   cephes_exp_p4 = _mm256_set1_ps(1.6666665459E-1);
__m256   cephes_exp_p5 = _mm256_set1_ps(5.0000001201E-1);
__m256   fx;
__m256i  imm0;
__m256   one           = _mm256_set1_ps(1.0f);

        x     = _mm256_min_ps(x, exp_hi);
        x     = _mm256_max_ps(x, exp_lo);

  /* express exp(x) as exp(g + n*log(2)) */
        fx     = _mm256_mul_ps(x, cephes_LOG2EF);
        fx     = _mm256_round_ps(fx, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC);
__m256  z      = _mm256_mul_ps(fx, inv_LOG2EF);
        x      = _mm256_sub_ps(x, z);
        z      = _mm256_mul_ps(x,x);

__m256  y      = cephes_exp_p0;
        y      = _mm256_mul_ps(y, x);
        y      = _mm256_add_ps(y, cephes_exp_p1);
        y      = _mm256_mul_ps(y, x);
        y      = _mm256_add_ps(y, cephes_exp_p2);
        y      = _mm256_mul_ps(y, x);
        y      = _mm256_add_ps(y, cephes_exp_p3);
        y      = _mm256_mul_ps(y, x);
        y      = _mm256_add_ps(y, cephes_exp_p4);
        y      = _mm256_mul_ps(y, x);
        y      = _mm256_add_ps(y, cephes_exp_p5);
        y      = _mm256_mul_ps(y, z);
        y      = _mm256_add_ps(y, x);
        y      = _mm256_add_ps(y, one);

  /* build 2^n */
        imm0   = _mm256_cvttps_epi32(fx);
        imm0   = _mm256_add_epi32(imm0, _mm256_set1_epi32(0x7f));
        imm0   = _mm256_slli_epi32(imm0, 23);
__m256  pow2n  = _mm256_castsi256_ps(imm0);
        y      = _mm256_mul_ps(y, pow2n);
        return y;
}

int main(){
    int i;
    float xv[8];
    float yv[8];
    __m256 x = _mm256_setr_ps(11.0f, -12.0f, 13.0f ,-14.0f ,15.0f, -16.0f, 17.0f, -18.0f);
    __m256 y = exp256_ps(x);
    _mm256_store_ps(xv,x);
    _mm256_store_ps(yv,y);

 /* compare exp256_ps with the double precision exp from math.h, 
    print the relative error             */
    printf("i      x                     y = exp256_ps(x)      double precision exp        relative error\n\n");
    for (i=0;i<8;i++){ 
        printf("i = %i  x =%16.9e   y =%16.9e   exp_dbl =%16.9e   rel_err =%16.9e\n",
           i,xv[i],yv[i],exp((double)(xv[i])),
           ((double)(yv[i])-exp((double)(xv[i])))/exp((double)(xv[i])) );
    }
    return 0;
}

Следующая таблица дает представление о точности в определенных точках, сравнивая exp256_ps с двойной точностью exp от math.h, Относительная ошибка в последнем столбце.

i      x                     y = exp256_ps(x)      double precision exp        relative error

i = 0  x = 1.000000000e+00   y = 2.718281746e+00   exp_dbl = 2.718281828e+00   rel_err =-3.036785947e-08
i = 1  x =-2.000000000e+00   y = 1.353352815e-01   exp_dbl = 1.353352832e-01   rel_err =-1.289636419e-08
i = 2  x = 3.000000000e+00   y = 2.008553696e+01   exp_dbl = 2.008553692e+01   rel_err = 1.672817689e-09
i = 3  x =-4.000000000e+00   y = 1.831563935e-02   exp_dbl = 1.831563889e-02   rel_err = 2.501162103e-08
i = 4  x = 5.000000000e+00   y = 1.484131622e+02   exp_dbl = 1.484131591e+02   rel_err = 2.108215155e-08
i = 5  x =-6.000000000e+00   y = 2.478752285e-03   exp_dbl = 2.478752177e-03   rel_err = 4.380257261e-08
i = 6  x = 7.000000000e+00   y = 1.096633179e+03   exp_dbl = 1.096633158e+03   rel_err = 1.849522682e-08
i = 7  x =-8.000000000e+00   y = 3.354626242e-04   exp_dbl = 3.354626279e-04   rel_err =-1.101575118e-08