在AVX2中高效实现log2(__ m256d)

Question

__m256d _mm256_log2_pd (__m256d a)除了英特尔之外,SVML 还没有其他编译器可用,他们表示其性能在AMD处理器上是有缺陷的.在g ++ - 4.8中缺少AVX日志内在函数(_mm256_log_ps)中的一些互联网实现？和SSE和AVX的SIMD数学库,但它们似乎比AVX2更多的SSE.还有Agner Fog的矢量库,但是它是一个包含更多东西的大型库,它只是向量log2,所以从它的实现中很难找出向量log2操作的基本部分.

那么有人可以解释如何有效地实现log2()4个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);
}

与C++和汇编中的Logarithm结果相比,当前矢量实现速度快4倍,速度快std::log2()2.5倍std::log().

具体而言,使用以下近似公式:

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的VCL实现旨在实现非常高的准确性,使用技巧避免舍入错误,避免可能导致在可能的情况下添加小数量的事物.这在一定程度上掩盖了基本设计.

(Agner没有VCL的正式回购,但是那个定期从上游tarball版本更新).

有关更快的更近似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和的AVX512 _mm512_getmant_ps)后得到了很好的结果,一旦我仔细调整它.(它是商业项目的一部分,所以我认为我不能发布代码.)快速的近似实现float正是我想要的.

在我的用例中,每个jrf_fastlog()都是独立的,因此OOO执行很好地隐藏了FMA延迟,甚至不值得使用VCL polynomial_5()函数使用的更高ILP更短时延多项式评估方法("Estrin方案",它做了一些非FMA在FMA之前相乘,导致更多的总指令).

如果您的项目可以使用GPL许可代码并且您需要高精度,则应该直接使用VCL.它只是标题,因此只有您实际使用的部分才会成为二进制文件的一部分.

VCL的log浮动和双重功能都在vectormath_exp.h.该算法有两个主要部分:

• 提取指数位并将该整数转换回浮点数(在调整IEEE FP使用的偏差之后).

• 在一些指数位中提取尾数和OR,以得到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,VCL log_d()使用了两个五阶多项式的比率. @ harold说这通常对精度有好处.一个与大量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;

你可以放弃两步的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一个有趣的方法:它键入puns a floatto uint32_t,并将该整数转换为float.由于IEEE binary32浮点数将指数存储在比尾数更高的位中,因此得到的结果float主要表示指数的值,按比例缩放1 << 23,但也包含来自尾数的信息.

  然后它使用具有几个系数的表达式来修复并获得log()近似值.它包括(constant + mantissa)在将浮点位模式转换为时校正尾数污染的除法float.我发现使用AVX2在HSW和SKL上的矢量化版本比使用4阶多项式的JRF fastlog更慢且更不准确.(特别是当它用作快速的一部分时,arcsinh它也使用除法单元vsqrtps.)

Answer 2

最后，这是我的最佳结果，在 Ryzen 1800X @3.6GHz 上，在单个线程中每秒给出约 8 亿个对数（每个对数有 2 亿个向量），并且精确到尾数中的最后几位。剧透：最后看看如何将性能提高到每秒 8.7 亿对数。

特殊情况：负数、负无穷大和NaN带负号位的 s 被视为非常接近 0（导致一些垃圾大负“对数”值）。正无穷大和NaN带正号位的 s 会产生大约 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;
}

它使用查找表方法和一阶多项式的组合，主要在维基百科上描述（链接位于代码注释中）。我可以在这里分配 8KB 的 L1 缓存（这是每个逻辑核心可用的 16KB 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，因此性能提高到每秒 8.7 亿对数。我将保持当前的实现不变，因为一旦 CPU 开始gather正确执行操作（像 GPU 那样进行合并），它会更惯用并且速度会更快。

EDIT2：在 Ryzen CPU（但不是 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)));