平方根倒数计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

这是在《雷神之锤 III 竞技场》源代码中的一个函数，功能是计算 平方根倒数.

在 3D 游戏场景的计算机图形学中，要求取照明和投影的光照与反射效果，就需要对一个曲面上的多个点坐标计算平方根倒数，即

浮点数的表示

浮点数在线转换工具：IEEE-754 Floating Point Converter

在 C 语言中，浮点数表示采用 IEEE 754 标准。对一个 32 位二进制浮点数，从高到低分成三部分

  • [31:31]：表示符号位(Sign, S)；
  • [30:23]：表示指数(Exponent, E)；
  • [22:0]：表示尾数(Mantissa, M)；

以小数 为例，

  • 非负数，符号位为 ；
  • 8 比特的指数位表示 区间内的整数，为了表示小于 的浮点数，指数需要取负数，所以将取值修改在 这个区间；因为 ，所以 ；
  • 23 比特偏移量，取值范围在区间 ，将 这个区间分成了 段， 表示从 开始的第 段，，

最终得到 32 位二进制转换为浮点数的计算公式为

其中 ，得到 的二进制表示为

0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 1 0 0 0 0 1 1

带入公式得到近似值为

平方根倒数计算

的定义域与值域均为，，于是浮点数的计算公式可以简化为 该 32 比特二进制数表示的整型为 。 两边同时取对数，得到 将 替换为 ， 令 ，则 整理得

其中 为表示浮点数 的 32 比特二进制的整型， 为表示浮点数 的 32 比特二进制的整型， 是一个常数。

文首程序中 i = * ( long * ) &y; 首先将浮点数的 32 比特二进制转换为整型，然后 i = 0x5f3759df - ( i >> 1 ) 与刚刚推导出的公式对应，。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

对于一个输入的浮点数 ，首先得到首次近似值 ，然后代入牛顿迭代公式 。

对于 ，，则 对应代码 y = y * ( threehalfs - ( x2 * y * y ) );.

只搞懂了这段代码的数学原理，还是不知道 这个数从何而来……

参考资料