【发布时间】：2018-01-27 22:24:22
【问题描述】：
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 ，但是它是一个大型库，其中包含更多的东西，只是矢量 log2，所以从它的实现中很难找出矢量 log2 操作的基本部分。
那么有人能解释一下如何有效地为 4 个double 数字的向量实现log2() 操作吗？ IE。就像 __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 的结果相比，当前向量实现比std::log2() 快4 倍，比std::log() 快2.5 倍。
【问题讨论】：
-
你不能只使用avx_mathfun中的log函数，然后将结果乘以所需的常量吗？
-
@PaulR，这是给
float，而不是double。至少我不知道如何获得像cephes_log_p0这样的常量来表示双数：github.com/reyoung/avx_mathfun/blob/master/avx_mathfun.h -
啊 - 没有注意到 - 在这种情况下，我建议寻找 SSE 解决方案（例如，参见 this question and its answers） - 将 SSE 实现扩展到 AVX 应该很容易。
-
@SergeRogatch 您可以使用您最喜欢的多项式拟合方法来生成这些常数。
-
英特尔最近将 AVX2 和 FMA 优化的数学函数添加到 glibc phoronix.com/…
标签： c++ algorithm floating-point logarithm avx2