在 C++ 中以优化方式计算 0<=x<=1 和 1<y<2 的 pow(x,y)答案

【问题标题】：Calculating in an optimized way pow(x,y) for 0<=x<=1 and 1<y<2 in C++在 C++ 中以优化方式计算 0<=x<=1 和 1<y<2 的 pow(x,y)
【发布时间】：2016-12-05 11:40:53
【问题描述】：

对于我的一个项目，我需要对pow(x,y) 进行非常快速的计算。希望它被限制在一个精确的领域，但它也需要足够的内存效率，如果它不够快的话。

就像我说的，它在 x 介于 0 和 1 之间，y 介于 1 和 2 之间的短范围内。因此，它必须在整个范围上足够精确，以便在递归调用时（和不要在一个数字上拖延）

如果你们遇到这样的事情或有建议......

【问题讨论】：

@Mehrdad 已编辑
exp(y*log(x)) 速度太慢怎么办？ ——你认为什么精度是足够的？ -- 试试x^y = x*(1+(x-1))^(y-1) = x*(1+(y-1)*(x-1)+(y-1)*(y-2)/2*(x-1)^2 + O(x-1)^3 作为近似值。
如果有人想出比 C++ 库中 pow() 的本机实现更好的东西，我会感到惊讶，这可能会使用本机 FPU 指令进行计算。
@SamVarshavchik ：如果您查看像 cise.ufl.edu/~cop4600/cgi-bin/lxr/http/source.cgi/lib/math/… 或 github.com/evanphx/ulysses-libc/blob/master/src/math/pow.c 这样的 pow 的 libc 实现的实际来源，您会感到惊讶。
@SamVarshavchik ：我的意思是，这取决于您使用“更好”解决的质量问题。通常的实现力求获得最精确的结果，因此可以通过牺牲一些准确性来获得一些速度。

标签： c++ c++11 math optimization

【解决方案1】：

将pow(x,y) 替换为

exp(y*log(x))

这有时也是 C 库实现，但对于许多整数输入会给出轻微失真的结果。因此，pow(x,y) 实现通常有开销来捕捉指数为 0 和 1 的微弱幂，整数幂，并进行一些其他转换以获得最精确的结果。减少这种开销可能已经足够加速了。

【讨论】：

你认为使用查找表优化日志会带来更多好处，因为 x 是一个定义的范围？顺便说一句，这就是我想要的答案
这完全取决于处理器并且必须进行测试。 exp(y*log(x)) 需要在 FPU 堆栈上加载两个参数，执行两个 FPU 内部操作并存储结果。查找表需要协调 FPU、CPU 和处理器缓存，这似乎更复杂。但这是针对标准 PC 架构的。在微处理器上，其他限制可能会导致其他偏好。
根据我的经验，这里建议的替换通常会将现代数学库的指数运算性能提高两倍（对于针对 x87 FPU 的旧库可能会降低）。
顺便说一句，STD log 和 exp 是否也有开销？
一个人必须编译一个小程序到 asm 来检查这个，但我认为如果有一个 FPU 命令，它将被使用，唯一的开销是堆栈操作和等待周期。

【解决方案2】：

对于x 和p 在这样的有界区间内，最佳方法是运行exp2() 和log2() 的两个优化近似值。由于来自exp2() 极小极大多项式的误差将呈指数复合，但对于0 < x ≤ 1，输入始终为p*log2(x) ≤ 0，因此该误差将表现良好。注意：x = 0, log2(0) = -∞。下面是powf() 例程的示例代码，以及1 ≤ p ≤ 2 和log2() 和exp2 的多重幂次方的误差图和多阶极小极大近似值，希望能满足您的绝对误差容限。

float aPowf(float x, float p){
    union { float f; uint32_t u; } L2x, e2e;
    float pL2, pL2r, pL2i, E2;

    if(x == 0)  return(0.0f);

    /* Calculate log2(x) Approximation */
    L2x.f = x;
    pL2 = (uint8_t)(L2x.u >> 23) - 127;                 // log2(m*2^e) = log2(m) + e
    L2x.u = (L2x.u & ~(0xFF << 23)) | (0x7F << 23);

    // Approximate log2(x) over 1 <= x < 2, use fma() fused multiply accumulate function for efficient evaluation
    // pL2 += -0xf.4fb9dp-4 + L2x.f;    // 4.303568e-2
    // pL2 += -0x1.acc4cap0 + L2x.f *  (0x2.065084p0 + L2x.f *  (-0x5.847fe8p-4));  // 4.9397e-3
    // pL2 += -0x2.2753ccp0 + L2x.f *  (0x3.0c426p0 + L2x.f *  (-0x1.0d47dap0 + L2x.f *  0x2.88306cp-4));   // 6.3717e-4
    pL2 += -0x2.834a9p0 + L2x.f *  (0x4.11f1d8p0 + L2x.f *  (-0x2.1ee4fcp0 + L2x.f *  (0xa.52841p-4 + L2x.f *  (-0x1.4e4cf6p-4)))); // 8.761e-5
    // pL2 += -0x2.cce408p0 + L2x.f *  (0x5.177808p0 + L2x.f *  (-0x3.8cfd5cp0 + L2x.f *  (0x1.a19084p0 + L2x.f *  (-0x6.aa30dp-4 + L2x.f *  0xb.7cb75p-8))));  // 1.2542058e-5
    // pL2 += -0x3.0a3514p0 + L2x.f *  (0x6.1cbb88p0 + L2x.f *  (-0x5.5737a8p0 + L2x.f *  (0x3.490a04p0 + L2x.f *  (-0x1.442ae8p0 + L2x.f *  (0x4.66497p-4 + L2x.f *  (-0x6.925fe8p-8))))));    // 1.8533e-6
    // pL2 += -0x3.3eb71cp0 + L2x.f *  (0x7.2194ep0 + L2x.f *  (-0x7.7cf968p0 + L2x.f *  (0x5.c642f8p0 + L2x.f *  (-0x2.faeb44p0 + L2x.f *  (0xf.9e012p-4 + L2x.f *  (-0x2.ef86f8p-4 + L2x.f *  0x3.dc524p-8)))))); // 2.831e-7
    // pL2 += -0x3.6c382p0 + L2x.f *  (0x8.23b47p0 + L2x.f *  (-0x9.f803dp0 + L2x.f *  (0x9.3b4f3p0 + L2x.f *  (-0x5.f739ep0 + L2x.f *  (0x2.9cb704p0 + L2x.f *  (-0xb.d395dp-4 + L2x.f *  (0x1.f3e2p-4 + L2x.f *  (-0x2.49964p-8))))))));  // 4.7028674e-8

    pL2 *= p;
//  if(pL2 <= -128)     return(0.0f);

    /* Calculate exp2(p*log2(x)) */
    pL2i = floorf(pL2);
    pL2r = pL2 - pL2i;
    e2e.u = ((int)pL2i + 127) << 23;

    // Approximate exp2(x) over 0 <= x < 1, use fma() fused multiply accumulate function for efficient evaluation.
    // E2 = 0xf.4fb9dp-4 + pL2r;    // 4.303568e-2
    // E2 = 0x1.00a246p0 + pL2r * (0xa.6aafdp-4 + pL2r * 0x5.81078p-4); // 2.4761e-3
    // E2 = 0xf.ff8fcp-4 + pL2r * (0xb.24b0ap-4 + pL2r * (0x3.96e39cp-4 + pL2r * 0x1.446bc8p-4));   // 1.0705e-4
    E2 = 0x1.00003ep0 + pL2r * (0xb.1663cp-4 + pL2r * (0x3.ddbffp-4 + pL2r * (0xd.3b9afp-8 + pL2r * 0x3.81ade8p-8)));   // 3.706393e-6
    // E2 = 0xf.ffffep-4 + pL2r * (0xb.1729bp-4 + pL2r * (0x3.d79b5cp-4 + pL2r * (0xe.4d721p-8 + pL2r * (0x2.49e21p-8 + pL2r * 0x7.c5b598p-12))));  // 1.192e-7
    // E2 = 0x1.p0 + pL2r * (0xb.17215p-4 + pL2r * (0x3.d7fb5p-4 + pL2r * (0xe.34192p-8 + pL2r * (0x2.7a7828p-8 + pL2r * (0x5.15bd08p-12 + pL2r * 0xe.48db2p-16)))));   // 2.9105833e-9
    // E2 = 0x1.p0 + pL2r * (0xb.17218p-4 + pL2r * (0x3.d7f7acp-4 + pL2r * (0xe.35916p-8 + pL2r * (0x2.761acp-8 + pL2r * (0x5.7e9f9p-12 + pL2r * (0x9.70c6ap-16 + pL2r * 0x1.666008p-16))))));  // 8.10693e-11
    // E2 = 0x1.p0 + pL2r * (0xb.17218p-4 + pL2r * (0x3.d7f7b8p-4 + pL2r * (0xe.35874p-8 + pL2r * (0x2.764dccp-8 + pL2r * (0x5.76b95p-12 + pL2r * (0xa.15ca6p-16 + pL2r * (0xf.94e0dp-20 + pL2r * 0x1.cc690cp-20)))))));    // 3.9714478e-11

    return(E2 * e2e.f);
}

一旦选择了适当的极小极大近似值，请确保使用融合乘法累加运算 fma() [这是单周期指令] 来实现 horner 多项式求值。

可以通过引入 LUT 来进一步提高精度，以减少范围以提高对数函数的精度。

【讨论】：