rootn (float_t x, int_t n)是一个计算第n个根x1 / n的函数，并受某些编程语言(例如OpenCL)支持。当使用IEEE-754浮点数时，假设只需要对标准化的操作数n进行处理，则可以基于对基础位模式的简单操作来生成任何x的有效的低精度起始近似值。root (x, n)的二进制指数将是x的二进制指数的1 / n。 IEEE-754浮点数的指数字段有偏。我们可以简单地将偏置的指数除以n来代替对指数进行无偏，除法和重新偏置结果，然后应用偏移量以补偿先前被忽略的偏斜。此外，除了提取然后划分指数字段外，我们还可以简单地划分整个操作数x，将其重新解释为整数。所需的偏移量很容易找到，因为任何n的参数1将返回结果1。
如果我们有两个可用的助手函数，__int_as_float()将ITU-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可以通过well-known optimizations of integer division by constant divisor减少为一个移位或乘加移位。一些可行的示例是: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的m = 2n * m时才是精确的。否则，与真实的数学结果相比，该近似将提供高估。通过稍微减小偏移量，我们可以将最大相对误差大约减半，从而导致低估和高估之间达到平衡。通过对最佳偏移量进行二分查找来轻松实现这一点，该最佳偏移量使用区间[1、2n]中的所有浮点数作为测试用例。这样做，我们发现: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's fast inverse square root的初始部分。基于观察原始原始偏移量与为最小化最大相对误差而优化的最终偏移量之间的差异，我可以制定一种用于二次校正的启发式方法，从而得出最终的优化偏移量值。
但是，我想知道是否可以通过某些闭式公式确定最佳偏移，以使相对误差的最大绝对值max(|(approx(x，n)-x1 / n)/ x1 /对于[1,2n)中的所有
x，n |)被最小化。为了便于说明，我们可以限制为binary32(IEEE-754单精度)数字。我知道，一般来说，没有关于极大极小值逼近的闭式解，但是我的印象是，对于第n个根等代数函数的多项式逼近，确实存在闭式解。在这种情况下，我们有(逐段)线性近似。
最佳答案
这是一些 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)。从实验上看，最大发生在
x为2的幂时。在这种情况下，有效项为零和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而言。希望以后还会有更多。关于algorithm - 优化的 `rootn(x, n)`低精度近似值，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/32042673/