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

Question

我做了一个函数来计算atan2(y,x)的定点近似.问题在于运行整个函数所需的~83个循环,70个循环(在AMD FX-6100上使用gcc 4.9.1 mingw-w64 -O3进行编译)完全由一个简单的64位整数除法完成!遗憾的是,该分裂的任何条款都不变.我能加速分裂吗？有什么方法可以删除它吗？

我想我需要这个除法,因为我用1D查找表近似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在其浮点执行单元中仅具有分频器(它们在硅区域方面非常昂贵),因此需要将整数转换为浮点并再返回.因此,只需用浮动除法代替整数1就不太可能有所帮助.您需要重构代码以使用浮点代替以利用更快的浮点除法.

如果您能够重构代码RCPSS,那么如果您不需要确切的答案,您也可以从大致的浮点互惠指令中受益.它的延迟大约为5个周期.