平方根表及算法
玛丽莲梦兔
960次浏览
2021年02月01日 07:29
最佳经验
本文由作者推荐
淘宝双十一-怎样发邮件
近期看张开川的那个
Ruby
学习文档的
7.2
节的知识就从网上查 询了一下，分享
上来
好久没用到过平方根之类的了，其实平方根算法也不怎么复杂的。
（武汉火麒麟）
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
平方根
97
计算方法一：
我们用
a
来表示
A
的平方根，方程
x-a=0
的解就为
A
的平方根
a。两边平方后
有：
x*x-2ax+A=0
，因为
x
不等于0
，两边除以
x
有：
x-2a+A/x=0
、
a=(x+A/x)/2
所以你只需设置 一个约等于
(x+A/x)/2
的初始值，代入上面公式，可以得到一个
更加近似的值 。
再将它代入，
又可以得到一个更加精确的值
……
依此方法，
最后< br>得到一个足够精度的
(x+A/x)/2
的值即为
A
的平方根值。
真的是这样吗
?
假设我们代入的值
x
﹤
a
由于这里考虑
a
﹥
0
故：
x*x
﹤
a*a
即
x
﹤
A/x
(x+A/x)/2
﹥
(x+x)/2
即（
x+A/x)/2>x
即当代入的
x
﹤
a
时
(x+A/x)/2
的值将比
x
大。
同样可以证明当代入的
x
﹥
a
时
(x+A/x)/2
的值将比
x
小。
这样随着计算次数的 增
加，
(x+A/x)/2
的值就越来越接近
a
的值了。
如
:
计算
sqrt(5)
设初值为
x
=
2
第一次计算：
(2+5/2)/2=2.25
第二次计算：
(2.25+5/2.25)/2=2.236111
第三次计算：
(2.236111+5/2.236111)/2=2.236068
这三步所得的结果和
5
的平方根值相差已经小于
0.001
了。
计算方法二：
我们可以使用二分法来计算平方根。
设
f(x)=x*x
-
A
同样设置
a
为
A
的平方根，哪么
a< br>就是
f(x)=0
的根。
你可以先找两个正值
m,n
使
f(m)<0,f(n)>0
根据函数的单调性
,a
就在区间
(m,n)
间。
然后计算
(m+n)/2,
计算
f((m+n)/2),
如果它大于零
,
那么
a
就在区间
(m,(m+n)/2)
之
间。
小于零
,
就在
((m+n)/2,n)
之间
,
如果等于零
,
那么
(m+n)/2
当然就是
a
。这样重复几
次
,
你可以把
a
存在的范围一步步缩小
,
在最后足 够精确的区间内随便取一个值
,
它就约等于
a
。
计算方法三：
以上的方法都不是很直接，在上世纪
80
年代的初中 数学书上，都还在介绍一种
比较直接的计算方法：
(1)
如求
54 756
的算术平方根时先由个位向左两位两位地定位
:
定位为
5,47,56 ,
接
着象一般除法那样列出除式
.
(2)
先从最高位用最 大平方数试商
:
最大平方数不超过
5
的是
2,
得商后
,
除式
5-4
后
得
1
。把商
2
写上除式 上。
(3)
加上下一位的数
:
得
147
。
(4)
用
20
去乘商后去试商
147:2×
20=40 < br>这
40
可试商为
3,
那就把试商的
3
加上
4
0
去除
147
。得
147÷
43=3,
把
3
写上除式上。这时
147-129=18
。
(5)
加上下一位的数
:
得
1856
。
(6)
用
20
去乘商后去试商
1856:23×
20=460 这
460
可试商为
4,
那就把试商的
4
加到
4 60
去除
1856
。得
4,
把
4
写上除式上。这时
1856-1856=0,
无余数啦。
(7)
这时除式上的商是< br>234,
即是
54756
的平方根。
哪么这种计算方法是怎 么得来的呢？查找了好久都没有找到答案。
静下心来仔细
分平方根的计算过程，后来的步骤都有
20
乘以也有的商再加上预计的商乘上预
计的商。设也有的商为
a
预 计的商为
b
就是
(20*a+b)*b
即
20ab+b*b
。而实质上
预计的商是平方根中已有的商的后一位数字，
平方根实际为
10a+b再乘以
10
的
N
次方（
N
为整数），这里我们可以简化 为平方根为
10a+b
（因为乘
10
的
N
次
方只影 响平方的小数点位置，对数字计算没有影响）。
这下终于明白了
,
设
a
为
A
的平方根的前
n
位，
b
为
A的平方根的
n
位后面的数
字，哪么（
10a+b)
就是
A
的平方根。有：
(10a+b)(10a+b)=100a*a+20ab+b*b=
A
变形后：
(20a+b)b=A-100a*a
上面的计算 中第一次商
2
，然后从结果中减
4
实质就是
A-100a*a
第二次再预计商
3
再减去
(20*2+3)*3
实质就是：
A-100a*a-20ab-b*b
即：
A-(10a+b)(10a+b)
此时
10a+b
看作为新的已有商
a
，再求下一个
b< br>值。这样就可以一位一位地进
行平方根的求解了。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
快速平方根算法
Algorithm 2008-04-19 16:17:38
在
3D
图形编程中，经常要求平方根或平方根的倒数，例如：求向量的长度 或将向量归
一化。
C
数学函数库中的
sqrt
具有理想的精度，但对 于
3D
游戏程式来说速度太慢。我们
希望能够在保证足够的精度的同时，进一步提高速 度。
Carmack
在
QUAKE3
中使用了下面的算 法，它第一次在公众场合出现的时候，几乎震住了
所有的人。据说该算法其实并不是
Carma ck
发明的，它真正的作者是
Nvidia
的
Gary
Tarolli
（未经证实）。
//
//
计算参数
x
的平方根的倒数
//
float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0
×
5f3759df - (i >> 1); //
计算第一个近似根
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x); //
牛顿迭代法
return x;
}
该算法的本质其实就是牛顿迭代法（
Newton-Raphson
Method< br>，简称
NR
），而
NR
的基础
则是泰勒级数（
Tay lor Series
）。
NR
是一种求方程的近似根的方法。首先要估计一个
与方程的根比较靠近的数值，
然后根据公式推算下一个更加近似的数值，
不断重复直到
可以获得满意的精度。其公式如下：
函数：
y=f(x)
其一阶导数为：
y
’
=f
’
(x)
则方程：
f(x)=0
的第
n+1
个近似根为
x[n+1] = x[n] - f(x[n]) / f
’
(x[n])
NR
最关键的地方在于估计第一个近似根。如果该近似根与真根足够靠近的话，那么 只
需要少数几次迭代，就可以得到满意的解。
现在回过头来看看如何利用 牛顿法来解决我们的问题。求平方根的倒数，实际就是求
方程
1/(x^2)-a=0
的解。将该方程按牛顿迭代法的公式展开为：
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
将
1/2
放到括号里面，就得到了上面那个函数的倒数第二行。
接着，我们要设法估计第一个近似根。这也是上面的函数最神奇的地方。它通过某种
方法算出了 一个与真根非常接近的近似根，
因此它只需要使用一次迭代过程就获得了较
满意的解。它是怎样 做到的呢？所有的奥妙就在于这一行：
i = 0x5f3759df - (i >> 1); //
计算第一个近似根
超级莫名其妙的语句，
不是吗？但仔细想一下的话，
还是可以理解的。
我们知道，
IEEE
标准下，
float
类型的数据在
32
位系统上是这样表示的（大体来说就是这样，但 省略了
很多细节，有兴趣可以
BAIDU
）：
bits
：
31 30
…
0
31
：符号位
30-23
：共
8
位，保存指数（
E
）
22-0
：共
23
位，保存尾数（
M
）
所以，
32
位的浮点数用十进制实数表示就是：
M*2^E
。开根然 后倒数就是：
M^(-1/2)*2^(-E/2)
。现在就十分清晰了。语句
i> >1
其工作就是将指数除以
2
，实现
2^(E/2)
的部分。
而前面用一个常数减去它，
目的就是得到
M^(1/2)
同时反转所有指数的
符号。
至于那个
0
×
5f3759df
，呃 ，我只能说，的确是一个超级的
Magic Number
。
那个
Magic
Number
是可以推导出来的，但我并不打算在这里讨论 ，因为实在太繁琐了。
简单来说，其原理如下：因为
IEEE
的浮点数中，尾数
M
省略了最前面的
1
，所以实际的尾
数是
1+M
。如果你 在大学上数学课没有打瞌睡的话，那么当你看到
(1+M)^(-1/2)
这样的
形式 时，
应该会马上联想的到它的泰勒级数展开，
而该展开式的第一项就是常数。
下面给出简单的推导过程：
对于实数
R>0
，假设其在
IEEE
的浮点表示中，
指数为
E
，尾数为
M
，则：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
将
(1+M)^(-1/2)
按泰勒级数展开，取第一项，得：
原式
= (1-M/2) * 2^(-E/2)
= 2^(-E/2) - (M/2) * 2^(-E/2)
如果不考虑指数的符号的话，
(M/2)*2^(E/2)
正是
(R>>1)
，
而在
IEEE
表示中，指数的符号只需简单地加上一个偏移即可，
而式子的前半部分刚好是个常数，所以原式可以转化为：
原式
= C - (M/2)*2^(E/2) = C - (R>>1)
，其中
C
为常数
所以只需要解方程：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
= C - (R>>1)
求出令到相对误差最小的
C
值就可以了
上面的推导过程只是我个人的理解，并未得到证实。而
Chris Lomont
则在 他的论文中
详细讨论了最后那个方程的解法，并尝试在实际的机器上寻找最佳的常数
C
。有兴趣的
朋友可以在文末找到他的论文的链接。
所以，所谓的
Magic
Number
，并不是从
N
元宇 宙的某个星系由于时空扭曲而掉到地球上
的，
而是几百年前就有的数学理论。
只要熟悉
NR
和泰勒级数，
你我同样有能力作出类似
的优化。
在
上有人做过测试，该函数的相对误差约为
0.177585%
，速度比
C
标准库
的
sqrt
提高超过
20%
。如 果增加一次迭代过程，相对误差可以降低到
e-
004
的级数，但