为了账号安全,请及时绑定邮箱和手机立即绑定

请教如何在AVX2中高效实现log2(__ m256d)

/ 猿问

请教如何在AVX2中高效实现log2(__ m256d)

C++
largeQ 2019-12-31 01:00:59

SVML __m256d _mm256_log2_pd (__m256d a)在除Intel以外的其他编译器上不可用,他们说它的性能在AMD处理器上受到限制。g ++-4.8中缺少 Internet上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的 “更有效的系列” 。如何改善?(同时需要提高性能和准确性)

到目前为止,我的实现每秒可以执行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 ++和汇编语言中对数的结果相比,当前向量实现的速度是C + +的 4倍std::log2()和2.5倍std::log()。


查看完整描述

3 回答

?
慕桂英4014372

通常的策略是基于身份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代替了logVCL,而log-base-e被引入常量及其使用方式。阅读它可能是理解.NET 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的logfloat和double函数在中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()使用两个5阶多项式的比率。 @ 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一个有趣的方法:这类型双关语一float来uint32_t,和转换是整数float。由于IEEE binary32浮点数将指数存储在比尾数更高的位中,因此结果float大部分表示指数的值,以进行缩放1 << 23,但还包含尾数信息。


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



查看完整回答
反对 2020-01-06
?
交互式爱情

最后,这是我的最佳结果,它在Ryzen 1800X @ 3.6GHz上在单个线程中每秒提供约8亿个对数(每个向量有4个对数的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;

}

它结合了查找表方法和一阶多项式,主要在Wikipedia上进行了描述(链接在代码注释中)。我可以在这里分配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()after line 的尾部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)));



查看完整回答
反对 2020-01-06
?
哈士奇WWW

对于Intel CPU,使用i64gather_pd和256b向量将更有效,而不是压缩为__m128ifor i32gather_pd。例如使用exps32 = _mm_srli_epi64(x, 20 + 32);。(为什么要使用算术移位而不是逻辑移位?是否需要符号位广播)?不过,提取high32可能对Ryzen很有好处,因为256b矢量指令是2 oups,而不是1 oups。因此,您花费2 oups提取以节省4 oups的__m128i指令,而不是__m256i

查看完整回答
反对 2020-01-06

添加回答

回复

举报

0/150
提交
取消
意见反馈 帮助中心 APP下载
官方微信