优化的低精度近似`rootn(x,n)`

Question

rootn (float_t x, int_t n)是一个计算第n个根x ^1/n的函数,并且受一些编程语言(如 OpenCL)的支持.当使用IEEE-754浮点数时n,假设只x需要处理归一化的操作数,则可以基于对基础位模式的简单操作生成任何有效的低精度起始近似值.

二进制指数root (x, n)将是二进制指数的1/n x.IEEE-754浮点数的指数字段是有偏差的.我们可以简单地将偏置指数除以n,然后应用偏移量来补偿先前忽略的偏差,而不是对指数进行不偏置,将其除去,并重新偏置结果.此外,我们可以简单地划分整个操作数x,重新解释为整数,而不是提取,然后除去指数字段.找到所需的偏移是微不足道的,因为1的参数将返回1的结果n.

如果我们有两个辅助函数可供我们使用,__int_as_float()它将IEEE-754重新解释binary32为int32,并且__float_as_int()重新解释int32操作数为binary32,则我们以rootn (x, n)直接的方式得出以下低精度近似:

rootn (x, n) ~= __int_as_float((int)(__float_as_int(1.0f)*(1.0-1.0/n)) + __float_as_int(x)/n)

通过恒定除数的整数除法的公知优化,__float_as_int (x) / n可以将整数除法简化为移位或乘法加移位.一些有用的例子是:

rootn (x,  2) ~= __int_as_float (0x1fc00000 + __float_as_int (x) / 2)  // sqrt (x)
rootn (x,  3) ~= __int_as_float (0x2a555556 + __float_as_int (x) / 3)  // cbrt (x)
rootn (x, -1) ~= __int_as_float (0x7f000000 - __float_as_int (x) / 1)  // rcp (x)
rootn (x, -2) ~= __int_as_float (0x5f400000 - __float_as_int (x) / 2)  // rsqrt (x)
rootn (x, -3) ~= __int_as_float (0x54aaaaaa - __float_as_int (x) / 3)  // rcbrt (x)

对于所有这些近似,只有当x= 2 ^n*m为整数时,结果才是精确的m.否则,与真实的数学结果相比,近似将提供高估.我们可以通过略微减小偏移来将最大相对误差减半,从而导致低估和高估的平衡混合.这可以通过二进制搜索最佳偏移量来轻松实现,该最佳偏移量使用区间[1,2 ⁿ ] 中的所有浮点数作为测试用例.这样做,我们发现:

rootn (x, 2) ~= __int_as_float (0x1fbb4f2e + __float_as_int(x)/2) // max rel err = 3.47474e-2
rootn (x, 3) ~= __int_as_float (0x2a51067f + __float_as_int(x)/3) // max rel err = 3.15547e-2
rootn (x,-1) ~= __int_as_float (0x7ef311c2 - __float_as_int(x)/1) // max rel err = 5.05103e-2
rootn (x,-2) ~= __int_as_float (0x5f37642f - __float_as_int(x)/2) // max rel err = 3.42128e-2
rootn (x,-3) ~= __int_as_float (0x54a232a3 - __float_as_int(x)/3) // max rel err = 3.42405e-2

有些人可能会注意到,计算rootn (x,-2)基本上是Quake快速反平方根的初始部分.

基于观察原始原始偏移与最优偏移之间的差异以最小化最大相对误差,我可以为二次校正制定启发式,从而制定最终的优化偏移值.

但是,我想知道是否有可能通过一些闭式公式确定最佳偏移量,使得相对误差的最大绝对值max(|((x,n) - x ^1/n)/ x ^1/n |),x在[1,2 ⁿ ]中最小化.为了便于说明,我们可以限制为binary32(IEEE-754单精度)数字.

我知道一般来说,对于极小极大近似没有封闭形式的解决方案,但是我的印象是,对于像第n个根这样的代数函数的多项式近似,确实存在闭式解.在这种情况下,我们有(分段)线性近似.

Answer 1

这里有一些Octave(MATLAB)代码,它根据下面的猜想计算正n的偏移量.似乎适用于2和3,但我怀疑当n太大时,其中一个假设会崩溃.现在没时间调查了.

% finds the best offset for positive n
% assuming that the conjectures below are true
function oint = offset(n)
% find the worst upward error with no offset
efrac = (1 / log(2) - 1) * n;
evec = [floor(efrac) ceil(efrac)];
rootnvec = 2 .^ (evec / n);
[abserr i] = max((1 + evec / n) ./ rootnvec);
relerr = abserr - 1;
rootnx = rootnvec(i);
% conjecture: z such that approx(z, n) = 1
% should have the worst downward error
fun = @(o) relerr - o / rootnx + (1 / (1 + o * n) ^ (1 / n) - 1);
oreal = fzero(fun, [0 1]);
oint = round((127 * (1 - 1 / n) - oreal) * 2 ^ 23);

只有正n的部分答案 - 假设我们没有向下调整偏移量,我只是通过推测最坏的估计来略微蚕食.

让我们定义一个理想化的近似版本x in [1, 2^n).

rootn-A(x, n) = 1 + floor(lg(x))/n + ((x/2^floor(lg(x))) - 1) / n
                    ^^^^^^^^^^^^^^   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                    contribution of       contribution of the
                     the exponent             significand

我们想要最大化rootn-A(x, n) / x^(1/n).

实验似乎是当2 x的幂时最大值出现.在这种情况下,有效数字项为零floor(lg(x)) = lg(x),因此我们可以最大化

(1 + lg(x)/n) / x^(1/n).

替代y = lg(x)/n的,我们可以最大限度地(1 + y) / 2^y为y in [0, 1)这样n*y是一个整数.降低完整性条件,这是一个微积分练习,表明这个凹函数的最大值是y = 1/log(2) - 1,大约0.4426950408889634.因此x,2的幂的最大值是x = 2^floor((1/log(2) - 1) * n)或者x = 2^ceil((1/log(2) - 1) * n).我猜想其中一个实际上是全球最大值.

在低估结束,看来我们想要x这样的输出rootn(x, n)是1,至少小n.希望以后会有更多的事情发生.