Эффективная реализация log2 (__ m256d) в AVX2

Обычная стратегия основана на идентификаторе log(a*b) = log(a) + log(b), или в данном случае log2( 2^exponent * mantissa) ) = log2( 2^exponent ) + log2(mantissa). Или упрощая, exponent + log2(mantissa). Мантисса имеет очень ограниченный диапазон, от 1,0 до 2,0, поэтому полином для log2(mantissa) должен соответствовать только этому очень ограниченному диапазону. (Или, что то же самое, мантисса = от 0,5 до 1,0 и измените константу коррекции смещения экспоненты на 1).

Расширение ряда Тейлора является хорошей отправной точкой для коэффициентов, но обычно вы хотите минимизировать максимальную абсолютную ошибку (или относительную ошибку) в этом конкретном диапазоне, а коэффициенты ряда Тейлора, вероятно, оставят более низкий или высокий выброс в этом диапазоне. , вместо того, чтобы иметь максимальную положительную ошибку, почти совпадающую с максимальной отрицательной ошибкой. Таким образом, вы можете сделать то, что называется минимаксной подгонкой коэффициентов.

Если важно, чтобы ваша функция оценивала log2(1.0) как точное значение 0.0, вы можете организовать это, фактически используя mantissa-1.0 в качестве полинома, а не постоянный коэффициент. 0.0 ^ n = 0.0. Это также значительно улучшает относительную ошибку для входных данных, близких к 1.0, даже если абсолютная ошибка все еще мала.

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

нацелена на очень высокую точность VCL Агнера _9 для реализации VCL от Агнера Фога. использовать уловки, чтобы избежать ошибки округления, избегая вещей, которые могут привести к добавлению малого и большого числа, когда это возможно. Это несколько затемняет основной дизайн.

Для более быстрого и приблизительного float log() см. Реализацию полинома на http://jrfonseca.blogspot.ca/2008/09/fast-sse2-pow-tables-or-polynomials.html. Он не учитывает МНОГО дополнительных приемов повышения точности, которые использует VCL, поэтому его легче понять. Он использует полиномиальную аппроксимацию мантиссы в диапазоне от 1,0 до 2,0.

(Это настоящий трюк для log() реализаций: вам нужен только многочлен, который работает в небольшом диапазоне.)

Он уже просто делает log2 вместо log, в отличие от VCL, где log-base-e встроен в константы и как он их использует. Его чтение, вероятно, станет хорошей отправной точкой для понимания exponent + polynomial(mantissa) реализаций log().

Даже его версия с наивысшей точностью не является полной float точностью, не говоря уже о double, но вы можете уместить многочлен с большим количеством членов. Или очевидно, что соотношение двух многочленов работает хорошо; это то, что VCL использует для double.

Я получил отличные результаты от переноса функции SSE2 из JRF в AVX2 + FMA (и особенно AVX512 с _mm512_getexp_ps и _mm512_getmant_ps), как только я ее тщательно настроил. (Это было частью коммерческого проекта, поэтому я не думаю, что смогу опубликовать код.) Быстрая приблизительная реализация для float была именно тем, что я хотел.

В моем варианте использования каждый jrf_fastlog() был независимым, поэтому выполнение OOO красиво скрыло задержку FMA, и даже не стоило использовать метод оценки полинома с более высокой ILP и меньшей задержкой, который Функция polynomial_5() VCL использует (" Схема Эстрина ", которая выполняет некоторые не-FMA умножения перед FMA, что приводит к большему количеству инструкций).

VCL Agner Fog теперь лицензируется Apache, поэтому любой проект может просто включать его напрямую. Если вам нужна высокая точность, вы должны просто использовать VCL напрямую. Это только заголовок, просто встроенные функции, поэтому он не приведет к раздуванию вашего двоичного файла.

Функции VCL log float и double находятся в vectormath_exp.h. Алгоритм состоит из двух основных частей:

• извлеките биты экспоненты и преобразуйте это целое число обратно в число с плавающей запятой (после корректировки смещения, которое использует IEEE FP).

• извлеките мантиссу и ИЛИ в некоторых битах экспоненты, чтобы получить вектор из double значений в диапазоне [0.5, 1.0). (Или (0.5, 1.0], я забыл).

  Далее отрегулируйте это с помощью if(mantissa <= SQRT2*0.5) { mantissa += mantissa; exponent++;}, а затем mantissa -= 1.0.

  Используйте полиномиальное приближение к log(x) с точностью около x = 1,0. (Для double log_d() VCL использует соотношение двух полиномов 5-го порядка. @harold говорит, что это часто хорошо для точности. Одно деление, смешанное с большим количеством FMA, обычно не снижает пропускной способности, но оно имеет большую задержку, чем FMA. Использование vrcpps + итерации Ньютона-Рафсона обычно медленнее, чем просто использование vdivps на современном оборудовании. Использование коэффициента также создает больше ILP за счет параллельной оценки двух полиномов более низкого порядка вместо одного полинома высокого порядка и может снизить общую задержку по сравнению с одной длинной цепочкой деплоя для высокого порядка полином (который также будет накапливать значительную ошибку округления в этой длинной цепочке).

Затем добавьте exponent + polynomial_approx_log(mantissa), чтобы получить окончательный результат log (). VCL делает это в несколько этапов, чтобы уменьшить ошибку округления. ln2_lo + ln2_hi = ln(2). Он разделен на малую и большую константу, чтобы уменьшить ошибку округления.

// res is the polynomial(adjusted_mantissa) result
// fe is the float exponent
// x is the adjusted_mantissa.  x2 = x*x;
res  = mul_add(fe, ln2_lo, res);             // res += fe * ln2_lo;
res += nmul_add(x2, 0.5, x);                 // res += x  - 0.5 * x2;
res  = mul_add(fe, ln2_hi, res);             // res += fe * ln2_hi;

Вы можете отказаться от двухэтапного ln2 материала и просто использовать VM_LN2, если вы не стремитесь к точности 0,5 или 1 ulp (или что-то еще, что эта функция фактически обеспечивает; IDK).

Я полагаю, что часть x - 0.5*x2 - это действительно дополнительный полиномиальный член. Это то, что я имел в виду под запеканием логарифмической базы e: вам понадобится коэффициент на этих условиях или избавиться от этой строки и заново подогнать полиномиальные коэффициенты для log2. Вы не можете просто умножить все коэффициенты полинома на константу.

После этого он проверяет отсутствие переполнения, переполнения или денормации и ветвления, если какой-либо элемент в векторе требует специальной обработки для получения правильного NaN или -Inf, а не любого мусора, который мы получили из полинома + экспонента. Если известно, что ваши значения конечны и положительны, вы можете закомментировать эту часть и получить значительное ускорение (даже проверка перед переходом требует нескольких инструкций).

Дальнейшее чтение:

• http://gallium.inria.fr/blog/fast-vectorizable-math-approx/ кое-что о том, как оценивать относительную и абсолютную ошибку в полиномиальном приближении и выполнять минимаксную фиксацию коэффициентов вместо простого использования разложения в ряд Тейлора.

• http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html интересный подход: он набирает float на uint32_t, а преобразует это целое число в float. Поскольку двоичные 32 числа с плавающей запятой IEEE хранят показатель степени в более высоких битах, чем мантисса, результирующий float в основном представляет собой значение показателя, масштабированное на 1 << 23, но также содержащее информацию из мантиссы.

  Затем он использует выражение с парой коэффициентов, чтобы исправить ситуацию и получить приближение log(). Он включает деление на (constant + mantissa), чтобы исправить загрязнение мантиссы при преобразовании битового шаблона с плавающей запятой в float. Я обнаружил, что векторизованная версия этого была медленнее и менее точна с AVX2 на HSW и SKL, чем JRF fastlog с полиномами 4-го порядка. (Особенно при использовании его как части быстрого arcsinh, в котором также используется блок деления для vsqrtps.)

Наконец, вот мой лучший результат, который на Ryzen 1800X @ 3,6 ГГц дает около 0,8 миллиарда логарифмов в секунду (200 миллионов векторов по 4 логарифма в каждом) в одном потоке и является точным до нескольких последних бит в мантиссе. Спойлер: посмотрим напоследок, как увеличить производительность до 0,87 миллиарда логарифмов в секунду.

Особые случаи: отрицательные числа, отрицательная бесконечность и NaN с битом отрицательного знака обрабатываются так, как если бы они были очень близки к 0 (что приводит к появлению большого количества отрицательных логарифмических значений для мусора). Положительная бесконечность и NaN с положительным битом знака приводят к логарифму около 1024. Если вам не нравится, как обрабатываются особые случаи, можно добавить код, который проверяет их и делает то, что вам больше подходит. Это замедлит вычисления.

namespace {
  // The limit is 19 because we process only high 32 bits of doubles, and out of
  //   20 bits of mantissa there, 1 bit is used for rounding.
  constexpr uint8_t cnLog2TblBits = 10; // 1024 numbers times 8 bytes = 8KB.
  constexpr uint16_t cZeroExp = 1023;
  const __m256i gDoubleNotExp = _mm256_set1_epi64x(~(0x7ffULL << 52));
  const __m256d gDoubleExp0 = _mm256_castsi256_pd(_mm256_set1_epi64x(1023ULL << 52));
  const __m256i cAvxExp2YMask = _mm256_set1_epi64x(
    ~((1ULL << (52-cnLog2TblBits)) - 1) );
  const __m256d cPlusBit = _mm256_castsi256_pd(_mm256_set1_epi64x(
    1ULL << (52 - cnLog2TblBits - 1)));
  const __m256d gCommMul1 = _mm256_set1_pd(2.0 / 0.693147180559945309417); // 2.0/ln(2)
  const __m256i gHigh32Permute = _mm256_set_epi32(0, 0, 0, 0, 7, 5, 3, 1);
  const __m128i cSseMantTblMask = _mm_set1_epi32((1 << cnLog2TblBits) - 1);
  const __m128i gExpNorm0 = _mm_set1_epi32(1023);
  // plus |cnLog2TblBits|th highest mantissa bit
  double gPlusLog2Table[1 << cnLog2TblBits];
} // anonymous namespace

void InitLog2Table() {
  for(uint32_t i=0; i<(1<<cnLog2TblBits); i++) {
    const uint64_t iZp = (uint64_t(cZeroExp) << 52)
      | (uint64_t(i) << (52 - cnLog2TblBits)) | (1ULL << (52 - cnLog2TblBits - 1));
    const double zp = *reinterpret_cast<const double*>(&iZp);
    const double l2zp = std::log2(zp);
    gPlusLog2Table[i] = l2zp;
  }
}

__m256d __vectorcall Log2TblPlus(__m256d x) {
  const __m256d zClearExp = _mm256_and_pd(_mm256_castsi256_pd(gDoubleNotExp), x);
  const __m256d z = _mm256_or_pd(zClearExp, gDoubleExp0);

  const __m128i high32 = _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(
    _mm256_castpd_si256(x), gHigh32Permute));
  // This requires that x is non-negative, because the sign bit is not cleared before
  //   computing the exponent.
  const __m128i exps32 = _mm_srai_epi32(high32, 20);
  const __m128i normExps = _mm_sub_epi32(exps32, gExpNorm0);

  // Compute y as approximately equal to log2(z)
  const __m128i indexes = _mm_and_si128(cSseMantTblMask,
    _mm_srai_epi32(high32, 20 - cnLog2TblBits));
  const __m256d y = _mm256_i32gather_pd(gPlusLog2Table, indexes,
    /*number of bytes per item*/ 8);
  // Compute A as z/exp2(y)
  const __m256d exp2_Y = _mm256_or_pd(
    cPlusBit, _mm256_and_pd(z, _mm256_castsi256_pd(cAvxExp2YMask)));

  // Calculate t=(A-1)/(A+1). Both numerator and denominator would be divided by exp2_Y
  const __m256d tNum = _mm256_sub_pd(z, exp2_Y);
  const __m256d tDen = _mm256_add_pd(z, exp2_Y);

  // Compute the first polynomial term from "More efficient series" of https://en.wikipedia.org/wiki/Logarithm#Power_series
  const __m256d t = _mm256_div_pd(tNum, tDen);

  const __m256d log2_z = _mm256_fmadd_pd(t, gCommMul1, y);

  // Leading integer part for the logarithm
  const __m256d leading = _mm256_cvtepi32_pd(normExps);

  const __m256d log2_x = _mm256_add_pd(log2_z, leading);
  return log2_x;
}

Он использует комбинацию подхода таблицы поиска и полинома 1-й степени, в основном описанного в Википедии (ссылка находится в комментариях к коду). Я могу позволить себе выделить здесь 8 КБ кэша L1 (что составляет половину от 16 КБ кеша L1, доступного на логическое ядро), потому что вычисление логарифмов на самом деле является узким местом для меня, и нет ничего, что требует кеш L1.

Однако, если вам нужно больше кеша L1 для других нужд, вы можете уменьшить объем кеша, используемый алгоритмом логарифмирования, уменьшив cnLog2TblBits до, например, 5 за счет уменьшения точности вычисления логарифма.

Или, чтобы сохранить высокую точность, вы можете увеличить количество полиномиальных членов, добавив:

namespace {
  // ...
  const __m256d gCoeff1 = _mm256_set1_pd(1.0 / 3);
  const __m256d gCoeff2 = _mm256_set1_pd(1.0 / 5);
  const __m256d gCoeff3 = _mm256_set1_pd(1.0 / 7);
  const __m256d gCoeff4 = _mm256_set1_pd(1.0 / 9);
  const __m256d gCoeff5 = _mm256_set1_pd(1.0 / 11);
}

А затем измените хвост Log2TblPlus() после строки const __m256d t = _mm256_div_pd(tNum, tDen);:

  const __m256d t2 = _mm256_mul_pd(t, t); // t**2

  const __m256d t3 = _mm256_mul_pd(t, t2); // t**3
  const __m256d terms01 = _mm256_fmadd_pd(gCoeff1, t3, t);
  const __m256d t5 = _mm256_mul_pd(t3, t2); // t**5
  const __m256d terms012 = _mm256_fmadd_pd(gCoeff2, t5, terms01);
  const __m256d t7 = _mm256_mul_pd(t5, t2); // t**7
  const __m256d terms0123 = _mm256_fmadd_pd(gCoeff3, t7, terms012);
  const __m256d t9 = _mm256_mul_pd(t7, t2); // t**9
  const __m256d terms01234 = _mm256_fmadd_pd(gCoeff4, t9, terms0123);
  const __m256d t11 = _mm256_mul_pd(t9, t2); // t**11
  const __m256d terms012345 = _mm256_fmadd_pd(gCoeff5, t11, terms01234);

  const __m256d log2_z = _mm256_fmadd_pd(terms012345, gCommMul1, y);

Затем следует комментарий // Leading integer part for the logarithm и остальные без изменений.

Обычно вам не нужно так много терминов, даже для нескольких битов таблицы, я просто предоставил коэффициенты и вычисления для справки. Скорее всего, если cnLog2TblBits==5, вам ничего не понадобится, кроме terms012. Но я таких замеров не проводил, нужно поэкспериментировать, что вам подходит.

Очевидно, что чем меньше полиномиальных членов вы вычисляете, тем быстрее будут вычисления.

ИЗМЕНИТЬ: этот вопрос В какой ситуации инструкции по сбору AVX2 будут быстрее, чем индивидуальная загрузка данных? предполагает, что вы можете улучшить производительность, если

const __m256d y = _mm256_i32gather_pd(gPlusLog2Table, indexes,
  /*number of bytes per item*/ 8);

заменяется на

const __m256d y = _mm256_set_pd(gPlusLog2Table[indexes.m128i_u32[3]],
  gPlusLog2Table[indexes.m128i_u32[2]],
  gPlusLog2Table[indexes.m128i_u32[1]],
  gPlusLog2Table[indexes.m128i_u32[0]]);

Для моей реализации он экономит около 1,5 цикла, уменьшая общее количество циклов для вычисления 4 логарифмов с 18 до 16,5, таким образом, производительность возрастает до 0,87 миллиарда логарифмов в секунду. Я оставляю текущую реализацию как есть, потому что она более идиоматична и должна быть быстрее, когда процессоры начнут правильно выполнять gather операции (с объединением, как это делают графические процессоры).

РЕДАКТИРОВАТЬ2: на процессорах Ryzen (но не на Intel) вы можете получить немного большее ускорение (около 0,5 цикла), заменив

const __m128i high32 = _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(
  _mm256_castpd_si256(x), gHigh32Permute));

с участием

  const __m128 hiLane = _mm_castpd_ps(_mm256_extractf128_pd(x, 1));
  const __m128 loLane = _mm_castpd_ps(_mm256_castpd256_pd128(x));
  const __m128i high32 = _mm_castps_si128(_mm_shuffle_ps(loLane, hiLane,
    _MM_SHUFFLE(3, 1, 3, 1)));