Самый простой способ - использовать экспоненциальное приближение.Один возможный случай, основанный на этом пределе
Для n = 256 = 2^8:
__m128d fastExp1(__m128d x)
{
__m128d ret = _mm_mul_pd(_mm_set1_pd(1.0 / 256), x);
ret = _mm_add_pd(_mm_set1_pd(1.0), ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
ret = _mm_mul_pd(ret, ret);
return ret;
}
Другая идея - разложение по полиному.В частности, разложение в ряд Тейлора:
__m128d fastExp2(__m128d x)
{
const __m128d a0 = _mm_set1_pd(1.0);
const __m128d a1 = _mm_set1_pd(1.0);
const __m128d a2 = _mm_set1_pd(1.0 / 2);
const __m128d a3 = _mm_set1_pd(1.0 / 2 / 3);
const __m128d a4 = _mm_set1_pd(1.0 / 2 / 3 / 4);
const __m128d a5 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5);
const __m128d a6 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5 / 6);
const __m128d a7 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5 / 6 / 7);
__m128d ret = _mm_fmadd_pd(a7, x, a6);
ret = _mm_fmadd_pd(ret, x, a5);
// If fma extention is not present use
// ret = _mm_add_pd(_mm_mul_pd(ret, x), a5);
ret = _mm_fmadd_pd(ret, x, a4);
ret = _mm_fmadd_pd(ret, x, a3);
ret = _mm_fmadd_pd(ret, x, a2);
ret = _mm_fmadd_pd(ret, x, a1);
ret = _mm_fmadd_pd(ret, x, a0);
return ret;
}
Обратите внимание, что при одинаковом количестве членов разложения вы можете получить лучшее приближение, если приблизитефункция для определенного диапазона x, используя, например, метод наименьших квадратов.
Все эти методы работают в очень ограниченном диапазоне x , но с непрерывными производными, которые могут быть важны в некоторых случаях,
Существует хитрость для аппроксимации показателя в очень широком диапазоне , но с заметными кусочно-линейными областями.Он основан на реинтерпретации целых чисел как чисел с плавающей точкой.Для более точного описания я рекомендую следующие ссылки:
Кусочно-линейное приближение к экспоненте и логарифму
Быстрое, компактное приближение экспоненциальной функции
Возможная реализация этого подхода:
__m128d fastExp3(__m128d x)
{
const __m128d a = _mm_set1_pd(1.0 / M_LN2);
const __m128d b = _mm_set1_pd(3 * 1024.0 - 1.05);
__m128d t = _mm_fmadd_pd(x, a, b);
return _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(t), 11));
}
Несмотря на простоту и широкий диапазон x для этого метода, будьте осторожны при использовании в математике.В небольших областях это дает кусочную аппроксимацию, которая может нарушить чувствительные алгоритмы, особенно те, которые используют дифференцирование.
Чтобы сравнить точность различных методов, посмотрите на графику.Первый график сделан для диапазона x = [0..1).Как видите, наилучшее приближение в этом случае дается методом fastExp2(x), чуть хуже, но приемлемым является fastExp1(x).Наихудшее приближение дает fastExp3(x) - кусочная структура заметна, разрывы первой производной есть.
В диапазоне x = [0..10) *Метод 1054 * обеспечивает наилучшую аппроксимацию, чуть хуже - аппроксимацию, заданную fastExp1(x) - при том же количестве вычислений он обеспечивает больший порядок, чем fastExp2(x).
Следующим шагом является повышение точности алгоритма fastExp3(x).Самый простой способ значительно повысить точность - использовать равенство exp(x) = exp(x/2)/exp(-x/2). Хотя это увеличивает объем вычислений, оно значительно уменьшает погрешность из-за взаимной компенсации ошибок при делении.
__m128d fastExp5(__m128d x)
{
const __m128d ap = _mm_set1_pd(0.5 / M_LN2);
const __m128d an = _mm_set1_pd(-0.5 / M_LN2);
const __m128d b = _mm_set1_pd(3 * 1024.0 - 1.05);
__m128d tp = _mm_fmadd_pd(x, ap, b);
__m128d tn = _mm_fmadd_pd(x, an, b);
tp = _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(tp), 11));
tn = _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(tn), 11));
return _mm_div_pd(tp, tn);
}
Еще большей точности можно достичь, комбинируя методы из алгоритмов fastExp1(x) или fastExp2(x) и fastExp3(x) с использованием равенства exp(x+dx) = exp(x) *exp(dx).Как показано выше, первый множитель может быть вычислен аналогично подходу fastExp3(x), для второго умножителя может использоваться метод fastExp1(x) или fastExp2(x).Поиск оптимального решения в этом случае является довольно сложной задачей, и я бы рекомендовал взглянуть на реализацию в библиотеках, предложенных в ответах.