【C】不使用sqrt函数计算开平方
不使用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语言表示：
#include <stdio.h>
#include <math.h>
#define ITERATIONTIMES (128) //最大精确迭代次数
double CustomSquare(double number)
{
int i;
if (number == 0.0)
return 0.0;
double result = number / 2.0; //迭代初始数
for (i = 1; i <= ITERATIONTIMES; i++)
{
result = 0.5 * (result + number / result);
}
return result;
}
int main(void)
{
double di;
printf("Input number:");
scanf("%lf", &di);
printf("Result by sqrt function:%1.9f\n", sqrtf(di));
printf("Result by custom sqrt function:%1.9f\n", CustomSquare(di));
}
本帖最后由 乘简 于 2017-9-28 15:40 编辑
网上抄了个代码，不用循环，也可以计算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
return y;
} 其实可以参考下牛顿迭代法 乘简 发表于 2017-9-28 15:38
网上抄了个代码，不用循环，也可以计算sqrt，精度与库函数一样，号称比库函数更快的算法。。。
这个代码中的 0x5f3759df 大有来头，有兴趣可以搜索一下，建议加个关键词 约翰卡马克 乘简 发表于 2017-9-28 15:38
网上抄了个代码，不用循环，也可以计算sqrt，精度与库函数一样，号称比库函数更快的算法。。。
然而你这个计算的不是sqrt，而是rsqrt，也就是reciprocal of the square root——倒数平方根（平方根倒数）。可用于向量的“标准化” 很好的选择 学习一下。
页:
[1]