什么“卧槽”的代码

七十二旹

已于 2023-12-05 20:28:35 修改

阅读量191

点赞数

文章标签：算法启发式算法

于 2023-11-07 23:59:22 首次发布

本文链接：https://blog.youkuaiyun.com/m0_72553226/article/details/134278050

版权

牛顿迭代公式

设 $r$ 是 $f(x)=0$ 的根，选取 $x_{0}$ 作为 $r$ 的初始近似值，过点 $(x_{0}$,$f(x_{0}))$ 做曲线 $y=f(x)$ 的切线 $L$ ，

$L:y=f(x_{0})+f'(x_{0})(x-x_{0})$ ，则 $L$ 与 $x$ 轴交点的横坐标 $x_{1}=x_{0}-\frac{f(x_{0})}{f'(x_{0})}$ ，称 $x_{1}$ 为 $r$ 的一次近似值。过点 $(x_{1},f(x_{1}))$ 做曲线 $y=f(x)$ 的切线，并求该切线与x轴交点的横坐标 $x_{2}=x_{1}-\frac{f(x_{1})}{f'(x_{1})}$ ，称 $x_{2}$ 为 $r$ 的二次近似值。重复以上过程，得 $r$ 的近似值序列，其中， $x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}$ 称为 $r$ 的 $n+1$ 次近似值，上式称为牛顿迭代公式。

用C语言写出：

#include<stdio.h>
#include<math.h>
double newton_sqrt(double a);
int main()
{
	double a,s;
	printf("请输入需要计算平方根的数：");
	scanf("%lf", &a);
	s=newton_sqrt(a);
	printf("%lf的平方根近似值为%8.5f", a,s);
	return 0;
}
double newton_sqrt(double a)
{
	double x = 1.0;//初始值
	while (1)
	{
		double next_x = 0.5 * (x + a / x);//牛顿迭代法
		if (fabs(next_x - x) <1e-7)//判断精度
		{
			return next_x;
		}
		x = next_x;
	}
	return 0;
}

但在我网上冲浪时，看到了来自Greg.Walsh格雷格·沃什的平方根倒数快速算法。

Amazing！！！

计算平方根是3D程序的基础。为了简化平方根代码，沃什开始了研究。

他利用《计算机组成原理》中对浮点数的存储方式（408第二章）（IEEE754标准）

$y=(1+\frac{M}{2^{23}})*2^{E-127}$ ,但这复杂程度不亚于牛顿法。

但沃什祭出对数大杀器

$log_{2}[(1+\frac{M}{2^{23}})*2^{E-127}]$

$=log_{2}(1+\frac{M}{2^{23}})=E-127$ 再利用近似值进行 $log_{2}(1+\frac{M}{2^{23}})\approx\frac{M}{2^{23}}$ 转换

$\approx \frac{M}{2^{23}}+E-127$

$=\frac{1}{2^{23}}(2^{23}*E+M)-127$

我们可以清楚看到这是二进制的表示方式

浮点数y,以二进制存储的Y

$log_{2}y=\frac{1}{2^{23}}*Y-127$

这已经很强了，但沃什不满足，他进行了又一次优化。

设a为y的平方根倒数

$log_{2}a=log_{2}y^{-\frac{1}{2}}=-\frac{1}{2}log_{2}y$

把a 和y换为它们的二进制形式A和Y

$\frac{1}{2^{23}}*A-127=-\frac{1}{2}(\frac{1}{2^{23}}*Y-127)$

二进制形式 $A=381*2^{22}-\frac{1}{2}*Y$

$381*2^{22}$ 转换为十六进制为5f400000

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;
}

其中0x5f3759df是对5f400000的修正

其中还有程序员之神约翰·卡马克的这句婉转悠扬的注释what the fuck?