【C】不使用sqrt函数计算开平方

cyycoish

不使用sqrt函数计算开平方
零、介绍
@元始天尊
和天尊兄第一次见面，他就给我出了一道难题：
不使用Sqrt函数，只能使用+-*/和流程控制，实现sqrt函数。
一、开始
我们先来实验一个数：sqrt(2)


然后设初始值为1，将1的代入公式：

得到x=1.5，将1.5再次代入公式：

得到x = 1.416666666，再次带入
第四次，所得到的x值已经为：x=1.414213562和同等小数精度的sqrt值相等
至此我们已经得到sqrt(2)的迭代公式：

由此，我们来推广sqrt(A),A属于任意实数的迭代公式：

我们来试一下：
假设A=5，求sqrt(A)的值：

反复带入：
第2次sqrt(5) = 2.333333.....
第3次sqrt(5) = 2.238095238
第4次sqrt(5) = 2.236068896
第5次sqrt(5) = 2.236067978
第6次sqrt(5) = 2.236067977
第7次sqrt(5) = 2.236067977
.........
所以，对于任意实数A，其sqrt(A)的迭代公式为：

二、关于迭代精度问题
假设迭代初始值都从1开始，我们分别对2、3、4、5、6、7、8、9进行计算
当迭代结果精度与sqrt函数返回值小数第九位的精度相等，表示我们已经达到目的
开始测试：
对象2
Sqrt(2) ~=1.414213562
1.Sqrt(2) ~=1.5
2. Sqrt(2) ~=1.416666667
3. Sqrt(2) ~=1.414215686
4. Sqrt(2) ~=1.414213562(OK)
对象 3
Sqrt(3) ~=1.732050808
1. Sqrt(3) ~=2
2. Sqrt(3) ~=1.75
3. Sqrt(3) ~=1.732142857
4. Sqrt(3) ~=1.73205081
5. Sqrt(3) ~=1.732050808(OK)
对象 4
Sqrt(4) =2
1. Sqrt(4) ~=2.5
2. Sqrt(4) ~=2.05
3. Sqrt(4) ~=2.000609756
4. Sqrt(4) ~=2.000000093
5. Sqrt(4) ~=2.000000000(OK)
对象 5
Sqrt(5) ~=2.236067977
1. Sqrt(5) ~=3
2. Sqrt(5) ~=2.333333333
3. Sqrt(5) ~=2.238095238
4. Sqrt(5) ~=2.236068896
5. Sqrt(5) ~=2.236067978
6. Sqrt(5) ~=2.236067977(OK)
对象 6
Sqrt(6) ~=2.449489743
1. Sqrt(6) ~=3.5
2 Sqrt(6) ~=2.607142857
3. Sqrt(6) ~=2.454235636
4. Sqrt(6) ~=2.449494372
5. Sqrt(6) ~=2.449489743(OK)
对象 7
Sqrt(7) ~=2.645751311
1. Sqrt(7) ~=4
2. Sqrt(7) ~=2.875
3. Sqrt(7) ~=2.654891304
4.Sqrt(7) ~=2.645767044
5. Sqrt(7) ~=2.645751311(OK)
对象 8
Sqrt(8) ~=2.828427125
1. Sqrt(8) ~=4.5
2. Sqrt(8) ~=3.138888889
3. Sqrt(8) ~=2.843780728
4. Sqrt(8) ~=2.828468572
5. Sqrt(8) ~=2.828427125(OK)
对象 9
Sqrt(9) =3
1.Sqrt(9) ~=5
2. Sqrt(9) ~=3.4
3. Sqrt(9) ~=3.023529412
4. Sqrt(9) ~=3.000091554
5. Sqrt(9) ~=3.000000001
6. Sqrt(9) ~=3.000000000(OK)
结论，可见，要使sqrt函数精确到小数点后第九位，迭代初始数为1，迭代次数5次以上即可。
三、算法的优化
我们说，迭代初始值从1开始还是太慢，计算量太大了。因为，如果我们计算sqrt到小数点后1千位呢？所以，这一小节我们试着探索优化迭代初始值。
那么，经过猜测迭代初始值应该越接近sqrt函数的返回值（也就是x的开平方）越好：
考虑到算法的复杂性要求，我们将x迭代初始值设为(x/2)。因为：对于任一实数r有sqrt(r)恒大于等于r/2（几何原理）。
那么根据迭代公式：

改变迭代初始值为:

我们继续试验：
这里我用5，9，3来进行测试，因为上面测试结果中5和9的最大精确迭代次数都超过了平均值，而数字3则是从先前测试结果大于等于平均值的对象中随机选取的。
对象 5
迭代初始值=5/2
Sqrt(5) ~=2.236067977
1. Sqrt(5) ~=2.25
2. Sqrt(5) ~=2.236111111
3. Sqrt(5) ~=2.236067978
4. Sqrt(5) ~=2.236067977(OK)
对象 9
迭代初始值=9/2
Sqrt(9) =3
1.Sqrt(9) ~=4.5
2. Sqrt(9) ~=3.25
3. Sqrt(9) ~=3.009615385
4. Sqrt(9) ~=3.00001536
5. Sqrt(9) ~=3(OK)
对象 3
迭代初始值=1.5
Sqrt(3) ~=1.732050808
1. Sqrt(3) ~=1.75
2. Sqrt(3) ~=1.732142857
3. Sqrt(3) ~=1.732050801
4. Sqrt(3) ~=1.732050808(OK)
那么根据以上结果，三个测试数最小精确迭代次数均<=先前平均值。
从而证明：迭代初始值设为x/2能够起到一部分减小运算量的效果。
因为测试数据量不够，所以此结果仅仅作为猜想。
四、算法的C语言表示：

1. #include <stdio.h>
   
2. #include <math.h>
   
3. 
   
4. #define ITERATIONTIMES (128) //最大精确迭代次数
   
5. 
   
6. double CustomSquare(double number)
   
7. {
   
8.         int i;
   
9.         if (number == 0.0)
   
10.                 return 0.0;
    
11.         double result = number / 2.0; //迭代初始数
    
12.         for (i = 1; i <= ITERATIONTIMES; i++)
    
13.         {
    
14.                 result = 0.5 * (result + number / result);
    
15.         }
    
16.         return result;
    
17. }
    
18. 
    
19. int main(void)
    
20. {
    
21.         double di;
    
22.         printf("Input number:");
    
23.         scanf("%lf", &di);
    
24.         printf("Result by sqrt        function:%1.9f\n", sqrtf(di));
    
25.         printf("Result by custom sqrt function:%1.9f\n", CustomSquare(di));
    
26. }
    

复制代码

不使用sqrt函数计算开平方.docx (20.05 KB, 下载次数: 15)

2015-9-26 14:37 上传

点击文件名下载附件

乘简

网上抄了个代码，不用循环，也可以计算sqrt，精度与库函数一样，号称比库函数更快的算法。。。

1. float Q_rsqrt( float number )  
   
2. {  
   
3.     long i;  
   
4.     float x2, y;  
   
5.     const float threehalfs = 1.5F;  
   
6.   
   
7.     x2 = number * 0.5F;  
   
8.     y  = number;  
   
9.     i  = * ( long * ) &y;            // evil floating point bit level hacking  
   
10.     i  = 0x5f3759df - ( i >> 1 ); // what the fuck?  
    
11.     y  = * ( float * ) &i;  
    
12.     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration  
    
13.     //      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed  
    
14.   
    
15.     return y;  
    
16. }

复制代码

唐凌

其实可以参考下牛顿迭代法

令狐公子

乘简 发表于 2017-9-28 15:38
网上抄了个代码，不用循环，也可以计算sqrt，精度与库函数一样，号称比库函数更快的算法。。。
[mw_shl_cod ...

这个代码中的 0x5f3759df 大有来头，有兴趣可以搜索一下，建议加个关键词 约翰卡马克

0xAA55

乘简 发表于 2017-9-28 15:38
网上抄了个代码，不用循环，也可以计算sqrt，精度与库函数一样，号称比库函数更快的算法。。。
[mw_shl_cod ...

然而你这个计算的不是sqrt，而是rsqrt，也就是reciprocal of the square root——倒数平方根（平方根倒数）。可用于向量的“标准化”

玫瑰花葬礼

很好的选择

Cer

学习一下。