x86_64平台下快速计算浮点数2的幂

Question

x86_64平台下快速计算浮点数2的幂

c++floating-pointvectorizationx86-64pow

8

是否有一种快速的方法将 2.0 提升到某个浮点数幂 x？我的意思是比使用 pow(2.0, x) 更快，最好能够与 AVX2 向量化。

对于整数，相应的方法是 1<<n，但它只适用于整数 n。

- Serge Rogatch

1

只是一个快速的想法，怎么样做 1<<n 并转换为双精度浮点数。如果这听起来很愚蠢，请原谅我 :-) - Malice

@Malice，问题在于“x”不是整数，而是浮点数。 - Serge Rogatch

1

对于整数，这只是写指数项。小数部分需要进行对数映射（或反向映射）到非指数部分。但实际上，如果可能实现大幅加速，则x^y = 2^(lg x * y)；因此，您可以希望的最大加速是消除一个基于2的对数和一个乘法。你希望做多少个2的幂次方？数百万个吗？ - Yakk - Adam Nevraumont

1

了解IEEE-754标准，您可以自己进行位操作... https://www.h-schmidt.net/FloatConverter/IEEE754.html - Tommy Andersen

这里可以找到使用“float”的天真方法：https://ideone.com/Bv0Y6V。虽然我不知道它是否更快，但显然这是针对整数幂的。 - Tommy Andersen

显示剩余8条评论

3个回答

3

对于整数次幂，您可以使用std::ldexp函数：

 double x = std::ldexp(1.0, k); // 2 to the power of k

这种方法比使用1<<k和强制类型转换更好，因为它不会出现中间溢出，并且还支持负数幂。

- Simon Byrne

挑战在于将2提高到浮点数幂x，其中x不是整数k。 - Serge Rogatch

如果k是一个整数且无需检查是否超出范围，您可以通过位移和加法手动将偏置指数填入IEEE 754 double中（通过偏置指数设置尾数），这样可以很好地进行自动向量化。 - Peter Cordes

2

TL;DR：将整数和小数分开，使用泰勒展开计算2^fraction部分，将其加回整数部分，然后使用许多人提到的转换为整数的魔术技巧。由于长时间的依赖关系计算了连续的幂，因此泰勒展开速度较慢，但该过程仅需要3个寄存器，这意味着您可以并行处理4个向量以允许流水线处理。

我知道这是一个古老的问题，但它从未得到很好的答案（甚至没有真正满足原问题者要求的任何答案）- 最佳答案是在评论中的一篇论文链接，该论文解释了泰勒级数和转换为整数技巧，并详细介绍了一种使用多项式调整您已经具有的小数部分以获取所需的2^fract的方法，但它很复杂，没有包括使其工作所需的实际数字。

原问题者想要每秒进行400,000,000次2^x计算，最少需要40位精度。我的要求有点不同-主要是我不需要那么高的精度，因此我尚未尝试使用double，但我在笔记本电脑上的10年旧的第三代i7移动芯片上以16位最小精度获得了1,730,000,000个2^x每秒的性能，在单个线程中，我相信在同一硬件上，具有双倍精度和系列中的两倍术语，其中1/4是可行的。

泰勒级数每步计算出e^x的更高精度，并且如下所示：

1 + x + x²/2! + x³/3! + (x^4)/4! + ...

要使其计算任何其他基数的幂，请将所需的幂乘以ln（power）。换句话说：

N ^ x = e^(x * ln(N))

或对于我们：

2 ^ x = e^（x * ln(2)）

如果使用大量的x值，它需要很长时间才能收敛，但由于我们对整数部分有一个技巧，因此我们永远不会处理x => 1.0。我们在6个术语（达到x ^ 6）中获得16位精度，在12个术语中获得40位，精度越高，离0越近，基本过程如下：

vector register x, y, z;   /// x is our input, y will become our output
y = floor(x);    /// 3 round toward -infinity
x -= y;          /// 3 split integer and fraction
x *= ln(2);      /// 5 ln(2) = 0.693147180559945...
y += x;          /// summing the terms direct into the integer part
z = x * x;       /// 5
y += z * (1/2);  /// this is a single fma instruction but ...
z *= x;          /// 5
y += z * (1/6);  /// if your processor lacks it, the working register ...
z *= x;          /// 5
y += z * (1/24); /// can be shared between the 4 parallel calculation chains
z *= x;          /// 5
y += z * (1/120); /// the broadcasted constants can be shared too
z *= x;           /// 5
y += z * (1/720); /// 8 I stopped here but you could easily carry on.
y *= 2^(mantissa bits);  /// 5        8388608.0 for float, 
                         /// 4503599627370496.0 for double   
y = convert-float-to-int(y); /// 3 creates an integer where the original
                             /// integer part appears where the exponent bits
                             /// of a float value would be, and our added
                             /// fraction part becomes the mantissa, gaining
                             /// the 1.0 it was missing in the process
y = integer-addition(y, interpret-bits-as-integer(1.0)); /// 3
                             /// add the exponent bias to complete our result
                             /// this can be done as a float addition before
                             /// the convert-to-int, but bits of the precision
                             /// of the result get lost, because all of the
                             /// mantissa and exponent bits have to fit into the
                             /// mantissa before conversion

每行可以转换为一个向量指令（在没有fma的旧处理器上为两个），但许多人必须等待近期步骤的结果完成（指令延迟-某些行后面的数字表示它们是链的一部分以及必须等待结果的周期数），这个链限制了这个6步版本的最小处理时间为55个周期。但是浮点硬件实际上可以比这快得多，用于进行24位乘法的5个加法阶段可以同时处理不同的值。在AVX寄存器文件中有16个寄存器，可以编码4个这样的计算链以单线程运行，允许在同一时间计算4个向量的结果...大约55个时钟周期。

使用双倍会减半每个链的结果数，但是翻倍的项数并不会使所需时间翻倍，因为起始步骤和完成步骤相同。它导致大约85个周期的延迟链，因此理论上可以在第三代3Ghz英特尔处理器的单个线程上处理560,000,000个40位精确的2^double结果（包括读取输入和写入输出）。新的处理器已经稍微降低了乘法的延迟，但是舍入和浮点到整数转换的延迟已经增加，因此总体上仍然大致相同。

- Jeremy

是的，多项式逼近在尾数范围内以及指数的整数部分中的处理通常是一个不错的选择。但是，在多项式求值方面，通常要使用霍纳规则，这只是一连串的FMA。请参见Fastest Implementation of Exponential Function Using AVX 了解实际C语言内置函数的工作实现。在您的伪代码中，1/6作为C表达式等于0，而不像1.f/6。 - Peter Cordes

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- DAle · Accepted Answer

有一个标准的std::exp2(double n)

计算 n 次方的 2 。

在特定的环境中，exp2(x) 可能不比 pow(2.0, x) 更快，但它比一般的 pow 更具体。