找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 7842|回复: 6

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

[复制链接]

85

主题

175

回帖

3990

积分

用户组: 超级版主

No. 418

UID
418
精华
14
威望
53 点
宅币
1974 个
贡献
1582 次
宅之契约
0 份
在线时间
252 小时
注册时间
2014-8-9
发表于 2015-9-26 14:31:33 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有账号?立即注册→加入我们

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

然后设初始值为1,将1的代入公式:
2.JPG
得到x=1.5,将1.5再次代入公式:
3.JPG
得到x = 1.416666666,再次带入
第四次,所得到的x值已经为:x=1.414213562和同等小数精度的sqrt值相等
至此我们已经得到sqrt(2)的迭代公式:
4.JPG
由此,我们来推广sqrt(A),A属于任意实数的迭代公式:
5.JPG
我们来试一下:
假设A=5,求sqrt(A)的值:
6.JPG
反复带入:
第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)的迭代公式为:
6.JPG
二、关于迭代精度问题
假设迭代初始值都从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(几何原理)。
那么根据迭代公式:
7.JPG
改变迭代初始值为:
8.JPG
我们继续试验:
这里我用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. #define ITERATIONTIMES (128) //最大精确迭代次数

  4. double CustomSquare(double number)
  5. {
  6.         int i;
  7.         if (number == 0.0)
  8.                 return 0.0;
  9.         double result = number / 2.0; //迭代初始数
  10.         for (i = 1; i <= ITERATIONTIMES; i++)
  11.         {
  12.                 result = 0.5 * (result + number / result);
  13.         }
  14.         return result;
  15. }

  16. int main(void)
  17. {
  18.         double di;
  19.         printf("Input number:");
  20.         scanf("%lf", &di);
  21.         printf("Result by sqrt        function:%1.9f\n", sqrtf(di));
  22.         printf("Result by custom sqrt function:%1.9f\n", CustomSquare(di));
  23. }
复制代码

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

本帖被以下淘专辑推荐:

In the beginning I was not the best.
And the world was also not the best.
But I still know that I am who I am.
Because I think that it is good.
I have been working hard.
I have been keeping growth with the world.
And it was so.
回复

使用道具 举报

13

主题

49

回帖

513

积分

用户组: 大·技术宅

UID
2285
精华
0
威望
39 点
宅币
362 个
贡献
11 次
宅之契约
0 份
在线时间
36 小时
注册时间
2017-2-25
发表于 2017-9-28 15:38:45 | 显示全部楼层
本帖最后由 乘简 于 2017-9-28 15:40 编辑

网上抄了个代码,不用循环,也可以计算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. }
复制代码
回复 赞! 1 靠! 0

使用道具 举报

65

主题

113

回帖

1万

积分

用户组: 超级版主

OS与VM研究学者

UID
1043
精华
35
威望
789 点
宅币
8246 个
贡献
1094 次
宅之契约
0 份
在线时间
2053 小时
注册时间
2015-8-15
发表于 2017-10-16 10:17:56 | 显示全部楼层
其实可以参考下牛顿迭代法
回复 赞! 靠!

使用道具 举报

0

主题

1

回帖

15

积分

用户组: 初·技术宅

UID
3846
精华
0
威望
2 点
宅币
10 个
贡献
0 次
宅之契约
0 份
在线时间
0 小时
注册时间
2018-5-15
发表于 2018-5-15 00:44:42 | 显示全部楼层
乘简 发表于 2017-9-28 15:38
网上抄了个代码,不用循环,也可以计算sqrt,精度与库函数一样,号称比库函数更快的算法。。。
[mw_shl_cod ...


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

使用道具 举报

1109

主题

1649

回帖

7万

积分

用户组: 管理员

一只技术宅

UID
1
精华
244
威望
743 点
宅币
24180 个
贡献
46222 次
宅之契约
0 份
在线时间
2294 小时
注册时间
2014-1-26
发表于 2018-11-14 18:18:02 | 显示全部楼层
乘简 发表于 2017-9-28 15:38
网上抄了个代码,不用循环,也可以计算sqrt,精度与库函数一样,号称比库函数更快的算法。。。
[mw_shl_cod ...

然而你这个计算的不是sqrt,而是rsqrt,也就是reciprocal of the square root——倒数平方根(平方根倒数)。可用于向量的“标准化”
回复 赞! 靠!

使用道具 举报

0

主题

43

回帖

60

积分

用户组: 小·技术宅

UID
4536
精华
0
威望
2 点
宅币
12 个
贡献
0 次
宅之契约
0 份
在线时间
6 小时
注册时间
2018-12-6
发表于 2018-12-6 10:55:04 | 显示全部楼层
很好的选择
回复 赞! 靠!

使用道具 举报

1

主题

32

回帖

124

积分

用户组: 小·技术宅

UID
5683
精华
0
威望
6 点
宅币
79 个
贡献
0 次
宅之契约
0 份
在线时间
6 小时
注册时间
2020-3-9
发表于 2020-3-10 11:27:58 | 显示全部楼层
学习一下。
回复

使用道具 举报

QQ|Archiver|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )|网站地图

GMT+8, 2024-3-28 19:06 , Processed in 0.049712 second(s), 37 queries , Gzip On.

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表