c++ - 在AVX2中有效实现log2(__ m256d)

Question

SVML的__m256d _mm256_log2_pd (__m256d a)在除Intel以外的其他编译器上不可用，他们说它的性能在AMD处理器上受到限制。互联网上有一些AVX log intrinsics (_mm256_log_ps) missing in g++-4.8?和SIMD math libraries for SSE and AVX引用的实现，但是它们似乎比AVX2更具SSE。还有Agner Fog's vector library，但是它是一个很大的库，其中包含的内容远远超过了vector log2，因此从实现中很难仅找出vector log2操作的基本部分。

因此，有人可以解释一下如何有效地对4个log2()数字的 vector 实现double操作吗？即类似于__m256d _mm256_log2_pd (__m256d a)的功能，但可用于其他编译器，并且对于AMD和Intel处理器都相当有效。

编辑:在我当前的特定情况下，数字是0到1之间的概率，并且对数用于熵计算:对i的所有P[i]*log(P[i])求和的求和。 P[i]的浮点指数范围很大，因此数字可以接近0。我不确定精度，因此会考虑以30个尾数开头的任何解决方案，尤其是可调整的解决方案是首选。

EDIT2:这是到目前为止，根据https://en.wikipedia.org/wiki/Logarithm#Power_series的“更有效的系列”，我的实现。如何改善？ (性能和准确性都需要提高)

namespace {
  const __m256i gDoubleExpMask = _mm256_set1_epi64x(0x7ffULL << 52);
  const __m256i gDoubleExp0 = _mm256_set1_epi64x(1023ULL << 52);
  const __m256i gTo32bitExp = _mm256_set_epi32(0, 0, 0, 0, 6, 4, 2, 0);
  const __m128i gExpNormalizer = _mm_set1_epi32(1023);
  //TODO: some 128-bit variable or two 64-bit variables here?
  const __m256d gCommMul = _mm256_set1_pd(2.0 / 0.693147180559945309417); // 2.0/ln(2)
  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 gVect1 = _mm256_set1_pd(1.0);
}

__m256d __vectorcall Log2(__m256d x) {
  const __m256i exps64 = _mm256_srli_epi64(_mm256_and_si256(gDoubleExpMask, _mm256_castpd_si256(x)), 52);
  const __m256i exps32_avx = _mm256_permutevar8x32_epi32(exps64, gTo32bitExp);
  const __m128i exps32_sse = _mm256_castsi256_si128(exps32_avx);
  const __m128i normExps = _mm_sub_epi32(exps32_sse, gExpNormalizer);
  const __m256d expsPD = _mm256_cvtepi32_pd(normExps);
  const __m256d y = _mm256_or_pd(_mm256_castsi256_pd(gDoubleExp0),
    _mm256_andnot_pd(_mm256_castsi256_pd(gDoubleExpMask), x));

  // Calculate t=(y-1)/(y+1) and t**2
  const __m256d tNum = _mm256_sub_pd(y, gVect1);
  const __m256d tDen = _mm256_add_pd(y, gVect1);
  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 log2_y = _mm256_mul_pd(terms01234, gCommMul);
  const __m256d log2_x = _mm256_add_pd(log2_y, expsPD);

  return log2_x;
}

到目前为止，我的实现每秒可以执行405 268 490次操作，而且直到第8位数字为止看起来还是很精确的。使用以下功能衡量性能:

#include <chrono>
#include <cmath>
#include <cstdio>
#include <immintrin.h>

// ... Log2() implementation here

const int64_t cnLogs = 100 * 1000 * 1000;

void BenchmarkLog2Vect() {
  __m256d sums = _mm256_setzero_pd();
  auto start = std::chrono::high_resolution_clock::now();
  for (int64_t i = 1; i <= cnLogs; i += 4) {
    const __m256d x = _mm256_set_pd(double(i+3), double(i+2), double(i+1), double(i));
    const __m256d logs = Log2(x);
    sums = _mm256_add_pd(sums, logs);
  }
  auto elapsed = std::chrono::high_resolution_clock::now() - start;
  double nSec = 1e-6 * std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count();
  double sum = sums.m256d_f64[0] + sums.m256d_f64[1] + sums.m256d_f64[2] + sums.m256d_f64[3];
  printf("Vect Log2: %.3lf Ops/sec calculated %.3lf\n", cnLogs / nSec, sum);
}

与 Logarithm in C++ and assembly的结果相比，当前的 vector 实现比 std::log2()快4倍，比 std::log()快2.5倍。

具体来说，使用以下近似公式:

Answer 1

通常的策略基于身份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的输入的相对误差。

您需要它有多精确，以及在什么输入范围内？像往常一样，在精度和速度之间需要权衡取舍，但幸运的是，通过例如再增加一个多项式项(并重新拟合系数)，或通过舍弃一些舍入误差来避免。

Agner Fog's VCL implementation of log_d() 的目标是达到极高的准确性，它使用技巧来避免舍入错误，从而避免了可能会导致加或减的数量增加的事情。这使基本设计有些模糊。

有关更快更近似的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的方式。

将JRF的SSE2函数移植到AVX2 + FMA(尤其是带有_mm512_getexp_ps和_mm512_getmant_ps的AVX512)后，我进行了仔细的调整，就获得了出色的结果。 (这是一个商业项目的一部分，所以我认为我不能发布代码。)float的快速近似实现正是我想要的。

在我的用例中，每个jrf_fastlog()是独立的，因此OOO执行很好地隐藏了FMA延迟，甚至不值得使用VCL's polynomial_5() function使用的更高ILP较短延迟多项式评估方法("Estrin's scheme"，它执行一些非FMA在FMA之前相乘，从而产生更多的总指令)。

Agner Fog的VCL现在已获得Apache许可，因此任何项目都可以直接包含它。如果需要高精度，则应直接使用VCL。它仅是标题，只是内联函数，因此不会膨胀您的二进制文件。

VCL的log float和double函数位于 vectormath_exp.h 中。该算法有两个主要部分:

提取指数位并将该整数转换回浮点数(针对IEEE FP使用的偏差进行调整之后)。

提取尾数和某些指数位中的OR，以获取double范围内[0.5, 1.0)值的 vector 。 (或者(0.5, 1.0]，我忘记了)。

进一步使用if(mantissa <= SQRT2*0.5) { mantissa += mantissa; exponent++;}和mantissa -= 1.0进行调整。

使用log(x)的多项式近似值，其精确度约为x = 1.0。 (对于double，VCL的log_d()使用两个5阶多项式的比率。@harold says this is often good for precision。一个除法器与许多FMA混合通常不会影响吞吐量，但确实比FMA具有更高的延迟。使用vrcpps + Newton- Raphson迭代通常比仅在现代硬件上使用vdivps慢，使用比率还可以通过并行评估两个低阶多项式(而不是一个高阶多项式)来创建更多的ILP，并且与一个较长的dep链相比，可以降低总延迟。一个高阶多项式(也会沿该长链累积明显的舍入误差)。

然后添加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;

如果您不希望精度达到0.5或1 ulp(或此功能实际提供的任何功能； IDK)，则可以删除两步ln2的东西，仅使用VM_LN2。

我猜x - 0.5*x2部分确实是一个额外的多项式术语。这就是我引入对数基数e的意思:您需要在这些条件上使用一个系数，或者要摆脱那条线并重新拟合log2的多项式系数。您不能只将所有多项式系数乘以一个常数。

之后，它将检查下溢，溢出或异常，并分支该 vector 中的任何元素是否需要特殊处理以产生适当的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。由于IEEE binary32 float将指数存储在比尾数更高的位中，因此生成的float主要表示指数的值，由1 << 23缩放，但也包含尾数信息。

然后，它使用具有几个系数的表达式来修正问题并获得log()近似值。它包括对(constant + mantissa)的除法，以在将float位模式转换为float时校正尾数污染。我发现，与带有4阶多项式的JRF faSTLog相比，在HSW和SKL上使用AVX2进行矢量化处理的速度较慢且准确性较低。 (尤其是将其用作快速arcsinh的一部分时，它也对vsqrtps使用了除法单元。)