Очень быстрая функция приближенного логарифма (натуральный логарифм) в с ++?

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

На алгоритмическом уровне проверьте, действительно ли все вызовы трансцендентных функций необходимы. Возможно, существует математическое преобразование, которое требует меньшего количества вызовов функций или преобразует трансцендентные функции в алгебраические операции. Возможно ли, что какие-либо вызовы трансцендентных функций избыточны, например потому что вычисления без надобности включают и выключают логарифмическое пространство? Если требования к точности невысоки, можно ли все вычисления выполнить с одинарной точностью, используя float вместо double? На большинстве аппаратных платформ отказ от double вычислений может привести к значительному увеличению производительности.

Компиляторы, как правило, предлагают множество переключателей, влияющих на производительность кода с большим количеством чисел. Помимо увеличения общего уровня оптимизации до -O3, часто есть способ отключить ненормальную поддержку, то есть включить режим сброса до нуля или FTZ. Это дает преимущества в производительности на различных аппаратных платформах. Кроме того, часто существует быстрый математический флаг, использование которого приводит к небольшому снижению точности и устраняет накладные расходы на обработку особых случаев, таких как NaN и бесконечности, а также обработку errno. Некоторые компиляторы также поддерживают автоматическую векторизацию кода и поставляются с математической библиотекой SIMD, например компилятор Intel.

Специальная реализация функции логарифмирования обычно включает разделение двоичного аргумента с плавающей запятой x на экспоненту e и мантиссу m, так что x = m * 2 e, следовательно log(x) = log(2) * e + log(m). m выбирается так, чтобы он был близок к единице, поскольку это обеспечивает эффективные аппроксимации, например log(m) = log(1+f) = log1p(f) с помощью минимаксного многочлена приближение.

C ++ предоставляет функцию frexp() для разделения операнда с плавающей запятой на мантиссу и экспоненту, но на практике обычно используются более быстрые машинно-зависимые методы, которые манипулируют данными с плавающей запятой на битовом уровне, повторно интерпретируя их как целые числа того же размера. Приведенный ниже код для логарифма с одинарной точностью logf() демонстрирует оба варианта. Функции __int_as_float() и __float_as_int() обеспечивают переинтерпретацию int32_t в число с плавающей запятой IEEE-754 binary32 и наоборот. Этот код в значительной степени полагается на объединенную операцию умножения-сложения FMA, поддерживаемую непосредственно аппаратным обеспечением на большинстве современных процессоров, CPU или GPU. На платформах, где fmaf() отображается на программную эмуляцию, этот код будет неприемлемо медленным.

#include <cmath>
#include <cstdint>
#include <cstring>

float __int_as_float (int32_t a) { float r; memcpy (&r, &a, sizeof r); return r;}
int32_t __float_as_int (float a) { int32_t r; memcpy (&r, &a, sizeof r); return r;}

/* compute natural logarithm, maximum error 0.85089 ulps */
float my_logf (float a)
{
    float i, m, r, s, t;
    int e;

#if PORTABLE
    m = frexpf (a, &e);
    if (m < 0.666666667f) {
        m = m + m;
        e = e - 1;
    }
    i = (float)e;
#else // PORTABLE
    i = 0.0f;
    if (a < 1.175494351e-38f){ // 0x1.0p-126
        a = a * 8388608.0f; // 0x1.0p+23
        i = -23.0f;
    }
    e = (__float_as_int (a) - __float_as_int (0.666666667f)) & 0xff800000;
    m = __int_as_float (__float_as_int (a) - e);
    i = fmaf ((float)e, 1.19209290e-7f, i); // 0x1.0p-23
#endif // PORTABLE
    /* m in [2/3, 4/3] */
    m = m - 1.0f;
    s = m * m;
    /* Compute log1p(m) for m in [-1/3, 1/3] */
    r =             -0.130310059f;  // -0x1.0ae000p-3
    t =              0.140869141f;  //  0x1.208000p-3
    r = fmaf (r, s, -0.121483512f); // -0x1.f198b2p-4
    t = fmaf (t, s,  0.139814854f); //  0x1.1e5740p-3
    r = fmaf (r, s, -0.166846126f); // -0x1.55b36cp-3
    t = fmaf (t, s,  0.200120345f); //  0x1.99d8b2p-3
    r = fmaf (r, s, -0.249996200f); // -0x1.fffe02p-3
    r = fmaf (t, m, r);
    r = fmaf (r, m,  0.333331972f); //  0x1.5554fap-2
    r = fmaf (r, m, -0.500000000f); // -0x1.000000p-1  
    r = fmaf (r, s, m);
    r = fmaf (i,  0.693147182f, r); //  0x1.62e430p-1 // log(2)
    if (!((a > 0.0f) && (a < INFINITY))) {
        r = a + a;  // silence NaNs if necessary
        if (a  < 0.0f) r = INFINITY - INFINITY; //  NaN
        if (a == 0.0f) r = -INFINITY;
    }
    return r;
}

Как отмечено в комментарии к коду, приведенная выше реализация обеспечивает точно округленные результаты с одинарной точностью и имеет дело с исключительными случаями, соответствующими стандарту с плавающей запятой IEEE-754. Производительность можно еще больше повысить, исключив поддержку особых случаев, исключив поддержку денормальных аргументов и снизив точность. Это приводит к следующему примерному варианту:

/* natural log on [0x1.f7a5ecp-127, 0x1.fffffep127]. Maximum relative error 9.4529e-5 */
float my_faster_logf (float a)
{
    float m, r, s, t, i, f;
    int32_t e;

    e = (__float_as_int (a) - 0x3f2aaaab) & 0xff800000;
    m = __int_as_float (__float_as_int (a) - e);
    i = (float)e * 1.19209290e-7f; // 0x1.0p-23
    /* m in [2/3, 4/3] */
    f = m - 1.0f;
    s = f * f;
    /* Compute log1p(f) for f in [-1/3, 1/3] */
    r = fmaf (0.230836749f, f, -0.279208571f); // 0x1.d8c0f0p-3, -0x1.1de8dap-2
    t = fmaf (0.331826031f, f, -0.498910338f); // 0x1.53ca34p-2, -0x1.fee25ap-2
    r = fmaf (r, s, t);
    r = fmaf (r, s, f);
    r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2) 
    return r;
}

Взгляните на это обсуждение, принятый ответ относится к реализация функции для вычисления логарифмов на основе разложения Цекендорфа.

В комментариях к файлу реализации обсуждается сложность и некоторые уловки для достижения O (1).

Надеюсь это поможет!

Я векторизовал ответ @njuffa. натуральный журнал, работает с AVX2:

inline __m256 mm256_fmaf(__m256 a, __m256 b, __m256 c){
    return _mm256_add_ps(_mm256_mul_ps(a, b), c);
}

//https://stackoverflow.com/a/39822314/9007125
//https://stackoverflow.com/a/65537754/9007125
// vectorized version of the answer by njuffa
/* natural log on [0x1.f7a5ecp-127, 0x1.fffffep127]. Maximum relative error 9.4529e-5 */
inline __m256 fast_log_sse(__m256 a){

    __m256i aInt = *(__m256i*)(&a);
    __m256i e =    _mm256_sub_epi32( aInt,  _mm256_set1_epi32(0x3f2aaaab));
            e =    _mm256_and_si256( e,  _mm256_set1_epi32(0xff800000) );
        
    __m256i subtr =  _mm256_sub_epi32(aInt, e);
    __m256 m =  *(__m256*)&subtr;

    __m256 i =  _mm256_mul_ps( _mm256_cvtepi32_ps(e), _mm256_set1_ps(1.19209290e-7f));// 0x1.0p-23
    /* m in [2/3, 4/3] */
    __m256 f =  _mm256_sub_ps( m,  _mm256_set1_ps(1.0f) );
    __m256 s =  _mm256_mul_ps(f, f); 
    /* Compute log1p(f) for f in [-1/3, 1/3] */
    __m256 r =  mm256_fmaf( _mm256_set1_ps(0.230836749f),  f,  _mm256_set1_ps(-0.279208571f) );// 0x1.d8c0f0p-3, -0x1.1de8dap-2
    __m256 t =  mm256_fmaf( _mm256_set1_ps(0.331826031f),  f,  _mm256_set1_ps(-0.498910338f) );// 0x1.53ca34p-2, -0x1.fee25ap-2

           r =  mm256_fmaf(r, s, t);
           r =  mm256_fmaf(r, s, f);
           r =  mm256_fmaf(i, _mm256_set1_ps(0.693147182f),  r);  // 0x1.62e430p-1 // log(2)
    return r;
}

Это зависит от того, насколько точным вы должны быть. Часто журнал вызывается, чтобы получить представление о величине числа, что вы можете сделать по существу бесплатно, исследуя поле экспоненты числа с плавающей запятой. Это тоже ваше первое приближение. Я добавлю плагин для моей книги «Основные алгоритмы», в которой объясняется, как реализовать математические функции стандартной библиотеки из первых принципов.