魔鬼数字

目录

一，开方的近似计算

1，代码

2，原理

二，开方倒数的近似计算

1，代码

2，原理

三，除法换成乘法

1，5的倍数除5

2，任意数除5

---

进公司之前对于魔鬼数字的理解，是真魔鬼。

进公司之后发现，2也是魔鬼数字。

收录一下真魔鬼的魔鬼数字。

一，开方的近似计算

1，代码

float MagicSqrt(float x)
{
	float xhalf = x / 2;
	int i = *(int*)&x;
	i = 0x1fbd1df5 + (i >> 1);
	return *(float*)&i;
}

int main()
{
	cout << MagicSqrt(2);
	return 0;
}

输出1.47748

如果再加牛顿法进行精度优化：

float MagicSqrt(float x)
{
	float xhalf = x / 2;
	int i = *(int*)&x;
	i = 0x1fbd1df5 + (i >> 1);
	x = *(float*)&i;
	x = x / 2 + xhalf / x;
	x = x / 2 + xhalf / x;
	return x;
}

int main()
{
	cout << MagicSqrt(2);
	return 0;
}

输出1.41421

输入5000输出70.7107，误差不到0.0001

2，原理

根据浮点数的编码方案，可以表示成：

  ，其中s和e表示小数（0到1）和指数（-127到128）

当0<x<1时，，其中c是一个常数。

所以上式简化成

即

如果把浮点数的32位比特按照int来解读，那么绝对值就等于

所以上式化成 

用不同方法估算的c值都在0.05左右，https://www.lomont.org/papers/2003/InvSqrt.pdf

给出了0.0450465是怎么算出来的。

所以

二，开方倒数的近似计算

1，代码

float rsqrt(float number)
{
	float y = number;
	long i = *(long *)&y;
	i = 0x5f3759df - (i >> 1);
	y = *(float *)&i;
	return y;
}

int main()
{
	cout << rsqrt(2);
	return 0;
}

输出0.716215

2，原理

和开方原理类似，最后也是用到

其中c取0.0450465

三，除法换成乘法

1，5的倍数除5

求一个5的倍数除5等于多少，可以换成乘法

int main()
{
	cout << 100*3435973837;
	return 0;
}

输出20

原理很简单：

cout << 3435973837 * 5;  会输出1，发生了溢出截断。

那么这个数是怎么来的呢？也很简单，解不定方程

最小解是

除了5之外，其他所有奇数几乎都可以这么搞。

2，任意数除5

取x=(2^33 + 3)/5，即5x=2的33次方加5

对于任意i32范围的5a+b，0<=b<=4，

有(5a+b)x=(2^33 + 3)a+b(2^33 + 3)/5=2^33 * a + 3a+b(2^33 + 3)/5

其中3a+b(2^33 + 3)/5的绝对值一定是小于2^33的。

所以，(5a+b)x右移33位一定是等于a

gcc64位编译器就是这么算的，任意i32数num除以5可以表示成：

（1）取x=(2^33 + 3)/5=1717986919

（2）计算num*x，结果是个64位的整数，存放在64位寄存器中

（3）右移32位，这是一个高速的算子

（4）右移1位