快速近似浮点数除法

Question

在现代处理器上，浮点数除法在倒数吞吐量方面比浮点乘法慢一个数量级。

我想知道是否有任何算法可以计算快速近似值x/y，在给定一定的假设和容差水平的情况下。例如，如果您假设0 < x < y，并且愿意接受任何真实值的输出在真实值的10%以内，则是否有比内置的FDIV操作更快的算法？

- dshin

---

你考虑过改变算法的结构以避免使用除法，而不是尝试找到更快的除法技巧吗？ - Apples

每个分割的10%误差在重复使用时可能会导致指数级的错误，这值得吗？可能不值得。 - WorldSEnder

这是此链接的副本吗？我不知道它是否更快，但至少有了一些东西。 - harold

1. 掩盖除数的小数位，2) 取其8个最高有效位，3) 查找其倒数，4) 用倒数乘以被除数，5) 将除数的指数加到乘积的指数部分。

- user3528438

你可以使用RCPSS指令快速1/X除法（倒数）。 - phuclv

Answer 1

希望这能帮到你，因为这可能是你所寻找的最接近的内容。

__inline__ double __attribute__((const)) divide( double y, double x ) {
                                    // calculates y/x
    union {
        double dbl;
        unsigned long long ull;
    } u;
    u.dbl = x;                      // x = x
    u.ull = ( 0xbfcdd6a18f6a6f52ULL - u.ull ) >> (unsigned char)1;
                                    // pow( x, -0.5 )
    u.dbl *= u.dbl;                 // pow( pow(x,-0.5), 2 ) = pow( x, -1 ) = 1.0/x
    return u.dbl * y;               // (1.0/x) * y = y/x
}

参见：
关于倒数逼近的另一篇文章。
维基百科页面。

- Jack G

Answer 2

FDIV通常比FMUL慢得多，因为它不能像乘法一样被分段处理，需要多个时钟周期进行迭代收敛的硬件查找过程。

最简单的方法是将除数y与除数的倒数x的乘积识别为除法。不那么直观的部分是记住浮点值x = m * 2 ^ e及其倒数x^-1 = (1/m)*2^(-e) = (2/m)*2^(-e-1) = p * 2^q的近似值p = 2/m = 3-x，对于1<=m<2。这给出了逆函数的粗略分段线性逼近，但我们可以通过使用迭代的牛顿根查找方法来改善该逼近。

让w = f(x) = 1/x，这个函数的反函数f(x)通过解出关于w的方程x = f^(-1)(w) = 1/w得到。为了用根查找方法改进输出，我们必须首先创建一个函数，其零点反映所需的输出，即g(w) = 1/w - x, d/dw(g(w)) = -1/w^2。 w[n+1]= w[n] - g(w[n])/g'(w[n]) = w[n] + w[n]^2 * (1/w[n] - x) = w[n] * (2 - x*w[n]) w[n+1] = w[n] * (2 - x*w[n]), 当w[n]=1/x时，w[n+1]=1/x*(2-x*1/x)=1/x 然后将这些组件相加，以获得最终的代码片段：

float inv_fast(float x) {
    union { float f; int i; } v;
    float w, sx;
    int m;

    sx = (x < 0) ? -1:1;
    x = sx * x;

    v.i = (int)(0x7EF127EA - *(uint32_t *)&x);
    w = x * v.f;

    // Efficient Iterative Approximation Improvement in horner polynomial form.
    v.f = v.f * (2 - w);     // Single iteration, Err = -3.36e-3 * 2^(-flr(log2(x)))
    // v.f = v.f * ( 4 + w * (-6 + w * (4 - w)));  // Second iteration, Err = -1.13e-5 * 2^(-flr(log2(x)))
    // v.f = v.f * (8 + w * (-28 + w * (56 + w * (-70 + w *(56 + w * (-28 + w * (8 - w)))))));  // Third Iteration, Err = +-6.8e-8 *  2^(-flr(log2(x)))

    return v.f * sx;
}

- nimig18