Sqrt-如何计算根号2

二分法

二分法：通过二分的方式，不断计算比较，推出根号2的一个近似值。

[1, 2]
[1, 1.5]
....

缺点：性能不高

牛顿迭代法

y = x^2 - 2 的图像

其实，计算根号2就是计算与 X 轴的交点。

切线是曲线的逼近

曲线上任一一点 (a, a^2-2) 处的切线为 y=2ax-a^2-2

这条切线和 X 轴的交点的横坐标为 0.5*(a + a/2);

根据这个结论，我们可以去逼近结果。

逼近的过程

取 x=2

做直线 x=2，和曲线相交于点 A(2, 2);

过 A 做切线 L1, 带入上面的 0.5*(a + a/2)，可知与 X 轴相交于 (1.5, 0);

取 x=1.5

同理。

做直线 x=1.5，和曲线相交于点 B(1.5, 0.25);

过 B 做切线 L2, 带入上面的 0.5*(a + a/2)，不断逼近结果，直到达到我们的精度为止。

实现源码

js - 牛顿迭代法

//求n的算术平方根，参数n不能为负数
function sqrt(n) {
    //当n>=1时，从n开始迭代；
    //当n<1时，从1开始迭代
    let res = n >= 1 ? n : 1;
    while(res * res - n > 1e-8)
        res = 0.5 * (res + n / res);
    return res;
}

c - 更高效的计算方式

java 的 Math.sqrt() 实际调用的是 c 实现。

• java

/**
 * Returns the correctly rounded positive square root of a
 * {@code double} value.
 * Special cases:
 * <ul><li>If the argument is NaN or less than zero, then the result
 * is NaN.
 * <li>If the argument is positive infinity, then the result is positive
 * infinity.
 * <li>If the argument is positive zero or negative zero, then the
 * result is the same as the argument.</ul>
 * Otherwise, the result is the {@code double} value closest to
 * the true mathematical square root of the argument value.
 *
 * @param   a   a value.
 * @return  the positive square root of {@code a}.
 */
public static native double sqrt(double a);

• 雷神III 中代码

这段代码的作用就是求number的平方根，并且返回它的倒数。

经过测试，它的效率比上述牛顿法程序要快几十倍。

也比c++标准库的sqrt()函数要快好几倍。

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
	#ifndef Q3_VM #
	ifdef __linux__ assert( !isnan(y) ); // bk010122 - FPE?
	#endif
	#endif return y; 
}

0x5f3759df 是什么？

很多文章写到上面就停止了，可是 0x5f3759df 才正是精妙所在。

到底为什么是这个数字呢？

这个在左耳朵耗子的专栏中讲解过，本文就把这个过程按自己的理解复述一遍。

开始之前，我们需要一些计算机基础知识。

结论

一个小数的计算方式是下面的公式：

(-1)^S * (1+M/2^23) * 2^(E-127)

下面是 IEEE 754 标准的简介，可以跳过，直接进入 简化浮点数公式

计算机的浮点数表示

在IEEE标准中，浮点数在内存中的表示是将特定长度的连续字节的所有二进制位按特定长度划分为符号域，指数域和尾数域三个连续域。

float

float 类型在内存中占用的位数为: 1+8+23=32bits

double

double 类型在内存中占用的位数为: 1+11+52=64bits

第一位s代表符号位，1代表负数，0代表正数。

第二个域是指数域，对于单精度float类型，指数域有８位，可以表示　0-255个指数值。指数值规定了小数点的位置，小数点的移动代表了所表示数值的大小。但是，指数可以为正数，也可以为负数。为了处理负指数的情况，实际的指数值按要求需要加上一个偏差（Bias）值作为保存在指数域中的值，单精度数的偏差 值为 -127，而双精度double类型的偏差值为 -1023。比如，单精度指数域中的64 则表示实际的指数值 -63。 偏差的引入使得对于单精度数，实际可以表达的指数值的范围就变成-127 到 128 之间（包含两端）。我们不久还将看到，实际的指数值-127（保存为 全 0）以及 +128（保存为全1）保留用作特殊值的处理。这样，实际可以表达的有效指数范围就在 -126 和 +127 之间。

第三个域为尾数域，其中单精度数为 23 位长，双精度数为 52 位长。比如一个单精度尾数域中的值为: 00001001000101010101000, 第二个域中的指数值则规定了小数点在尾数串中的位置，默认情况下小数点位于尾数串首位之前。　　

例子

比如指数值为 -1,则该float数即为：.000001001000101010101000,如果为+1,则该float 数值为:0.0001001000101010101000。

我们知道引入浮点数的目的在于用尽可能少的位数表示既高精度又大范围的实数，其中的范围大小是由指数域位长确定的，而尾数域的长度则确定了所能表示实数的精度，所以double比float数的精度更高，范围更大，相应的也就占用更多的内存。　

刚才我们介绍的对尾数域中的值的解释并不能实现这个精度最大化的目标，因为在尾数串第一个”1”之前还有４个”0”,这4个”0”实际上是多余的，因为我们把小数点向前移动时，前端的＂０＂是自动添加的，所以可以把这４个“0”删除，然后尾数域末端多出４个位来表示更高精度的数值。

也就是说尾数的第一位一定是＂１＂，那么既然第一位一定是＂１＂，那么我们也就没有必要把它存储在尾数域中，而是直接默认尾数为1.xxxx…xxx的形式。

尾数的首位从小数点后开始。那么上面的例子所表示的尾数就是：

1.00001001000101010101000。　

用23位表示了24位的信息 (小数点不占位置).

表示

一个规格化的32位浮点数ｘ的真值为：

x=(−1)^s × (1.M) × 2^E−127

一个规格化的64位浮点数ｘ的真值为：

x=(−1)^s × (1.M) × 2^E−1023

下面举一个32位单精度浮点数-3.75表示的例子帮助理解：

(1) 首先转化为2进制表示

−3.75=−(2+1+1/2+1/4)=−1.111×2^1

(2) 整理符号位并进行规格化表示

−1.111×21=(−1)(1)×(1+0.1110 0000 0000 0000 0000 000)×2^1

(3) 进行阶码的移码处理

(−1)^(1)×(1+0.1110 0000 0000 0000 0000 000)×2^1

=(−1)^(1)×(1+0.1110 0000 0000 0000 0000 000)×2^128−127

于是，符号位S=1，尾数M为1110 0000 0000 0000 0000 000阶码E为128_10=1000 0000_2,

则最终的32位单精度浮点数为 1 1110 0000 0000 0000 0000 000 1000 0000

简化浮点数公式

(-1)^S * (1+M/2^23) * 2^(E-127) 这个公式比较复杂，我们简化一下。

令 m=(M/2^23)，e=(E-127)

∵ 符号位在 y=x^(-1/2) 的两端都是 0(正数)，可以省去。 ∴ 公式简化为 (1+m) * 2^e

上面的算式是从一个 32 bits 二进制计算机推导而来，整形算式为 M + E*2^23

案例

比如，0.015 的 32bits 的二进制是：00111100011101011100001010001111，也就是整型的：

7717519 + 120 * 2^23
= 1014350479
= 0X3C75C28F

平方根倒数的推导

平方根公式：y = x^(-1/2)

两边同时取以2为基数的对数：

log2(y) = -1/2 * log2(x)

我们在计算浮点数，所以将 x、y 都使用 (1+m) * 2^e 替换掉。

得到：

TODO…

参考资料

• 牛顿迭代

如何求根号2

梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现

牛顿迭代法(Newton’s Method)

如何通俗易懂地讲解牛顿迭代法求开方？

• 计算机浮点数表示

计算机浮点数规格化表示

浅谈计算机中浮点数的表达方法（IEEE 754）

• 二分法
• 牛顿迭代法
  • 切线是曲线的逼近
  • 逼近的过程
    • 取 x=2
    • 取 x=1.5
• 实现源码
  • js - 牛顿迭代法
  • c - 更高效的计算方式
• 0x5f3759df 是什么？
  • 结论
• 计算机的浮点数表示
  • float
  • double
  • 例子
  • 表示
• 简化浮点数公式
  • 案例
• 平方根倒数的推导
• 参考资料