【快乐数学】雷神III中平方根倒数程序分析

获取代码

代码在github上
git clone https://github.com/raspberrypi/quake3.git
代码目录

在这款游戏中有一个经典的计算平方根的倒数的代码

float Q_rsqrt( float number )//函数 求平方根的倒数

即：
$$y=\frac{1}{\sqrt{x}}$$

前言

在游戏中 模拟物理现象，例如计算光反射之类的常会用到
例如：雷神III的里面代码

code/renderer/tr_shade_calc.c

计算反射 光的alpha值 不重要 可以略过 goto 正文

/*
** RB_CalcSpecularAlpha
**
** Calculates specular coefficient and places it in the alpha channel
*/
vec3_t lightOrigin = { -960, 1980, 96 };		// FIXME: track dynamically

void RB_CalcSpecularAlpha( unsigned char *alphas ) {
	int			i;
	float		*v, *normal;
	vec3_t		viewer,  reflected;
	float		l, d;
	int			b;
	vec3_t		lightDir;
	int			numVertexes;

	v = tess.xyz[0];
	normal = tess.normal[0];

	alphas += 3;

	numVertexes = tess.numVertexes;
	for (i = 0 ; i < numVertexes ; i++, v += 4, normal += 4, alphas += 4) {
		float ilength;

		VectorSubtract( lightOrigin, v, lightDir );
//		ilength = Q_rsqrt( DotProduct( lightDir, lightDir ) );
		VectorNormalizeFast( lightDir );

		// calculate the specular color
		d = DotProduct (normal, lightDir);
//		d *= ilength;

		// we don't optimize for the d < 0 case since this tends to
		// cause visual artifacts such as faceted "snapping"
		reflected[0] = normal[0]*2*d - lightDir[0];
		reflected[1] = normal[1]*2*d - lightDir[1];
		reflected[2] = normal[2]*2*d - lightDir[2];

		VectorSubtract (backEnd.or.viewOrigin, v, viewer);
		ilength = Q_rsqrt( DotProduct( viewer, viewer ) );
		l = DotProduct (reflected, viewer);
		l *= ilength;

		if (l < 0) {
			b = 0;
		} else {
			l = l*l;
			l = l*l;
			b = l * 255;
			if (b > 255) {
				b = 255;
			}
		}

		*alphas = b;
	}
}

正文：

代码：
code\tools\lcc\src\c.h

.....
typedef union {
	float f;
	int i;
	unsigned int ui;
} floatint_t
.....

code/qcommon/q_math.c

/*
** float q_rsqrt( float number )
*/
float Q_rsqrt( float number )
{
	floatint_t t;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	t.f  = number;
	t.i  = 0x5f3759df - ( t.i >> 1 );               // what the fuck?
	y  = t.f;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

浮点数float存储值为I，浮点数真实值为val
$$\begin{gathered} I={(-1)}^S+2^{23}\times E+M \\ E=e+127 \\ M=m\times 2^{23} \\ val=(-1{)}^s\times (1+m)\times 2^e \end{gathered}$$
S是1bit （第31位）符号位 1/0
E是8bit （第23-30位）范围-128~127 即 正3需要加127 表示指数位
M是23bit（第0-22位）表示小数位
E与e以及M与m关系
$$\begin{gathered} E=e+B(B=127) \\ M=L\times m(L=2^{23}) \end{gathered}$$

举个栗子：
8.25在内存中存储的方式
$$\begin{gathered} 由val=(1+m)\times 2^e \\ 得到8.25=1000.01B=1.00001\times 2^{3}B(B表示是二进制) \\ 即： \\ e=3 \\ m=0.00001B(去掉最前面的1因为默认是1) \\ S=0 \\ E=127+3=10000010B \\ M=00001000000000000000000B \\ I={(-1)}^S+2^{23}\times E+M(将值填充进去) \\ I=0/10000010/00001000000000000000000B \\ I=01000001000001000000000000000000B \\ I=0X41040000=0X41040000 \end{gathered}$$
即 8.25在内存中存储的结构为

推到
$$\begin{gathered} y=\frac{1}{\sqrt{x}} \\ {\log }_{2}y={\log }_2{x}^{-\frac{1}{2}} \\ \because y\ge 0x\ge 0 \\ 且 \\ y=(1+m_y)\times 2^{e_y} \\ x=(1+m_x)\times 2^{e_x} \\ \therefore \\ {\log }_2((1+m_y)\times 2^{e_y})=-\frac{1}{2}\times ({\log }_2((1+m_x)\times 2^{e_x})) \\ lo{g}_2(1+m_y)+e_y=-\frac{1}{2}\times (lo{g}_2(1+m_x)+e_x)(公式1) \end{gathered}$$
研究公式：
$$\begin{gathered} f_{(x)}=lo{g}_2(1+x)xϵ[0,1] \\ 与之拟合的函数 \\ {f}_{(x)}=x+b \end{gathered}$$
公式1：
$$\begin{gathered} lo{g}_2(1+m_y)+e_y=-\frac{1}{2}\times (lo{g}_2(1+m_x)+e_x) \\ b+m_y+e_y=-\frac{1}{2}(b+m_x+e_x) \\ \because e=E-B且m=\frac{M}{L} \\ b+\frac{M_y}{L}+E_y-B=-\frac{1}{2}\times (b+\frac{M_x}{L}+E_x-B) \end{gathered}$$
整理公式得到
$$\begin{gathered} M_y+L\times E_y=\frac{3}{2}\times (B-b)-\frac{1}{2}\times (M_x+L\times E_x) \\ 即 \\ {I}_y=k-(I_x>>2) \end{gathered}$$

迭代：

验证部分拟合

利用windows自带的计算器进行画图得到图像
$$\begin{gathered} f_{(x)}=lo{g}_2(1+x)xϵ[0,1] \\ {f}_{(x)}=x+bxϵ[0,1] \end{gathered}$$
设置b为值0.0450615 经验

拟合情况很近似
关于b的取值0.0450165的值

https://zhuanlan.zhihu.com/p/79699768

验证浮点数

#include <stdio.h>
#include <stdlib.h>



int main(int argc, char * argv[])
{
    float  n = 8.25;

    printf("0x%x \n", *(unsigned long *)&n);
    
    return 0;
}

扩展发散：

求平方根

李永乐老师的教授的求平方根的办法：

$$\sqrt{s}=a+\frac{b}{2a+\frac{b}{2a+\frac{b}{2a+\frac{b}{2a+...}}}}$$

#include <stdio.h>
#include <stdlib.h>

#define SPEC 10 //迭代10次的精度

float q_sqrt(float num)
{
    int  i = 0;
    int  lim = (int) num/2;
    
    for(i = 0; i < lim ; i++)
    {
    
        if( i * i > num)
        {
            i --;
            break;
        }
    }

    int     a = i;
    float   b = num - i*i ;
    float   t = b / (2*a) ;
    
    for (i = 0; i < SPEC; i++)
    {
        t = b/(2*a + t);
    }
    
    return a+t;
}


int main(int argc, char * argv[])
{
    float ret = 0.0;
    
    ret = q_sqrt(3.0);

    printf("%f \n", ret);

    return 0;
}

雷神之锤