Skip to main content
 首页 » 编程设计

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

2024年02月20日21cloudgamer

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倍。

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

请您参考如下方法:

通常的策略基于身份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使用了除法单元。)