c - 从固定点 atan2() 近似中删除缓慢的 int64 除法

Question

我创建了一个函数来计算 atan2(y, x) 的定点近似值。问题在于，在运行整个函数所需的约 83 个周期中，70 个周期(在 AMD FX-6100 上使用 gcc 4.9.1 mingw-w64 -O3 编译)完全由一个简单的 64 位整数除法占用!遗憾的是，该部门的所有条款都不是不变的。我可以加快 split 本身吗？有什么方法可以删除它吗？

我想我需要这个除法因为我用一维查找表近似 atan2(y, x) 我需要将 x,y 表示的点的距离归一化为单位圆或单位正方形(我选择一个单位“菱形”，它是一个旋转 45° 的单位正方形，它在正象限上提供了相当均匀的精度)。所以除法找到 (|y|-|x|)/(|y|+|x|)。请注意，除数是 32 位的，而分子是一个 32 位的数字，右移 29 位，因此除法的结果有 29 位小数位。也不能使用浮点除法，因为此函数要求不使用浮点运算。

有什么想法吗？我想不出有什么可以改进的(而且我不明白为什么仅仅一个部门就需要 70 个周期)。完整功能供引用:

int32_t fpatan2(int32_t y, int32_t x)       // does the equivalent of atan2(y, x)/2pi, y and x are integers, not fixed point
{
    #include "fpatan.h" // includes the atan LUT as generated by tablegen.exe, the entry bit precision (prec), LUT size power (lutsp) and how many max bits |b-a| takes (abdp)
    const uint32_t outfmt = 32; // final output format in s0.outfmt
    const uint32_t ofs=30-outfmt, ds=29, ish=ds-lutsp, ip=30-prec, tp=30+abdp-prec, tmask = (1<<ish)-1, tbd=(ish-tp);   // ds is the division shift, the shift for the index, bit precision of the interpolation, the mask, the precision for t and how to shift from p to t
    const uint32_t halfof = 1UL<<(outfmt-1);    // represents 0.5 in the output format, which since it is in turns means half a circle
    const uint32_t pds=ds-lutsp;    // division shift and post-division shift
    uint32_t lutind, p, t, d;
    int32_t a, b, xa, ya, xs, ys, div, r;

    xs = x >> 31;       // equivalent of fabs()
    xa = (x^xs) - xs;
    ys = y >> 31;
    ya = (y^ys) - ys;

    d = ya+xa;
    if (d==0)       // if both y and x are 0 then they add up to 0 and we must return 0
        return 0;

    // the following does 0.5 * (1. - (y-x) / (y+x))
    // (y+x) is u1.31, (y-x) is s0.31, div is in s1.29

    div = ((int64_t) (ya-xa)<<ds) / d;  // '/d' normalises distance to the unit diamond, immediate result of division is always <= +/-1^ds
    p = ((1UL<<ds) - div) >> 1;     // before shift the format is s2.29. position in u1.29

    lutind = p >> ish;      // index for the LUT
    t = (p & tmask) >> tbd;     // interpolator between two LUT entries

    a = fpatan_lut[lutind];
    b = fpatan_lut[lutind+1];
    r = (((b-a) * (int32_t) t) >> abdp) + (a<<ip);  // linear interpolation of a and b by t in s0.32 format

    // Quadrants
    if (xs)             // if x was negative
        r = halfof - r;     // r = 0.5 - r

    r = (r^ys) - ys;        // if y was negative then r is negated

    return r;
}

Answer 1

不幸的是，对于 x86 CPU 上的 64 位整数除法，70 个周期的延迟是典型的。浮点除法通常具有大约一半或更少的延迟。成本增加的原因是现代 CPU 在其浮点执行单元中只有除法器(它们在硅面积方面非常昂贵)，因此需要将整数转换为浮点并再次转换回来。因此，仅用浮点除法代替整数除法不太可能有帮助。您需要重构您的代码以使用浮点而不是利用更快的浮点除法。

如果您能够重构您的代码，如果您不需要确切的答案，您也可以从近似浮点倒数指令 RCPSS 中获益。它有大约 5 个周期的延迟。