这个 sqrt 近似内联汇编函数是如何工作的?
Posted
技术标签:
【中文标题】这个 sqrt 近似内联汇编函数是如何工作的?【英文标题】:How does this sqrt approximation inline assembly function work? 【发布时间】:2017-01-21 22:54:06 【问题描述】:通读3D游戏编程大师的技巧,我发现了这个用内联汇编编写的排序函数:
inline float FastSqrt(float Value)
float Result;
_asm
mov eax, Value
sub eax, 0x3F800000
sar eax, 1
add eax, 0x3F800000
mov Result, eax
return(Result);
这是实际平方根的近似值,但准确度足以满足我的需要。
这实际上是如何工作的?这个神奇的0x3F800000
值是什么?我们如何通过减法、旋转和加法获得平方根?
这是它在 C/C++ 代码中的样子:
inline float FastSqrt_C(float Value)
float Result;
long Magic = *((long *)&Value);
Magic -= 0x3F800000;
Magic >>= 1;
Magic += 0x3F800000;
Result = *((float *)&Magic);
return(Result);
【问题讨论】:
0x3F800000 是 1.0 的 32 位浮点表示 有趣,这就是为什么当我将参数 Value 更改为 int 时我猜我会得到不正确的结果?看起来这个函数只适用于浮点数? 更重要的是,这是指数偏差。所以它消除了偏差,将指数减半,然后将偏差加回去。它也与尾数有点混淆。 @vexe 所涉及的一般原理可以看this question 【参考方案1】:很多人指出0x3f800000
是1.0
的表示。虽然这是真的,但它与计算的方式无关。要理解它,您需要知道非负浮点数是如何存储的。 f = (1+m)*2^x
,0 <= m < 1
和 m
是尾数,x
是指数。另请注意,x
存储有偏差,因此二进制文件中的实际内容是x+127
。 32 位值由符号位(在我们的例子中为零)、8 位指数存储 x+127
和最后 23 位尾数 m
组成。 (见wikipedia article)。
应用一些基本的数学,
sqrt(f) = sqrt((1+m)*2^x)
= sqrt(1+m)*sqrt(2^x)
= sqrt(1+m)*2^(x/2)
因此,作为粗略的近似,我们需要将指数减半,但由于存在偏差,我们不能只做x/2
,我们需要(x-127)/2 + 127
。这个127
移动到适当的位位置是神奇的0x3f800000
。
除以 2 是通过右移一位来实现的。由于这对整个浮点数起作用,因此它对尾数也有副作用。
首先,假设原始指数是偶数。然后移出的最低有效位为零。因此,尾数也减半,所以我们最终得到的是:sqrt(f) = (1+m/2)*2^(x/2)
。我们得到了正确的指数,但尾数是 (1+m/2)
而不是 sqrt(1+m)
。这方面的最大相对误差是(1.5 - sqrt(2))/sqrt(2) ~ 6%
,如果m
几乎是1
意味着f
接近,但小于2
的奇数幂,就会发生这种情况。以f=7.99
为例。该公式为我们提供了关于2.998
而不是2.827
确实有6%
的错误。
现在,如果指数是奇数,那么最低有效位将是1
,当移入尾数时会导致增加一半。因此,我们得到sqrt(f) = (1.5+m/2)*2^((x-1)/2)
。最大的错误实际上是m=0
,这将是(1.5/sqrt(2)-sqrt(1))/sqrt(1)
,它又是6%
。这发生在从上面接近 2 的奇数次幂的数字上。
如果输入值恰好接近 2 的奇数幂,则这两种情况相结合意味着最差的不准确度约为 6%。对于 2 的偶数次方,结果是准确的。
【讨论】:
【参考方案2】:0x3F800000 in float 是 1。这是因为浮点数的存储方式。您可以在https://gregstoll.dyndns.org/~gregstoll/floattohex/ 看到一个视觉表示。
这是一个很好的近似值,我相信 sqrt。这源于一个游戏 Quake for inverse sqrt (https://en.wikipedia.org/wiki/Fast_inverse_square_root#Aliasing_from_floating_point_to_integer_and_back)。
【讨论】:
如果你画出 y=(x+1)/2 和 y=sqrt(x) 的图,你会发现当 x 在 [1,2] 中时它们很接近。所以我猜这是该区间内值的近似值。 @Roadowl 它确实不计算(x+1)/2
您链接的***文章讨论了一种与问题中提供的代码完全不同的近似平方根的方法。【参考方案3】:
下面是这个机制的一个例子:
FastSqrt(4.0) == 2.0
4.0 to hex -> 0x40800000
0x40800000 - 0x3f800000 = 0x1000000
0x1000000 to binary -> 00000001 00000000 00000000 00000000
shift toward the lsb (sar) -> 00000000 10000000 00000000 00000000
00000000 10000000 00000000 00000000 back to hex -> 0x00800000
0x00800000 + 0x3f800000 = 0x40000000
0x40000000 to dec -> 2.0
【讨论】:
如果您还显示每个步骤的指数/尾数字段,而不是只是整个 binary32 位模式,这会更好。【参考方案4】:浮点数 f = (1 + m)* [2^(e+127)],其中 m 为尾数部分,e 为指数部分。
因此:sqrt(f) = (f)^(1/2) = ((1 + m)* [2^(e+127)] )^(1/2)
-> ((1 + m)* [2^(e+127)] )^(1/2) = (1 + m)^(1/2) * 2^((e + 127)/ 2)
在指数部分,2^((e + 127)/2):
2^((e + 127)/2) = 2^( (e-127/2) + 127)
因此,在浮动表示中, 它是 (e - 0x3F800000) /2 + 0x3F800000
尾数部分,(1 + m)^(1/2):
从二项级数公式,(1 + x)^r = 1 + rx + (r(r - 1)/2)*(x^2) +....
因此,(1 + m)^(1/2) 等于 (1 + m/2 - (m^2)/8 + ...) 它大约等于 1 + m/2(一阶的典型近似值) 因此,尾数部分应除以2。
但是,尾数和指数组合为 A 数,右移除以指数和尾数 BOTH。
要评估错误,您可以考虑二项式级数的第二项 - (m^2)/8。
因为 m 总是小于 1,我将 m 替换为 0.9999 (0.5 + 0.25 + 0.125 + ...)
(m^2)/8 = 0.12497500125,这是最坏的情况。
【讨论】:
以上是关于这个 sqrt 近似内联汇编函数是如何工作的?的主要内容,如果未能解决你的问题,请参考以下文章
实用技能分享,充分利用内联函数,内联汇编,内部函数和嵌入式汇编提升代码执行效率和便捷性(2021-12-17)
实用技能分享,充分利用内联函数,内联汇编,内部函数和嵌入式汇编提升代码执行效率和便捷性(2021-12-17)