计算机与代数---如何计算sqrt---方法和实现

最新推荐文章于 2021-05-29 00:38:37 发布

ARTELE

最新推荐文章于 2021-05-29 00:38:37 发布

阅读量3.1k

点赞数

分类专栏：计算机与数学文章标签：计算根号数学 sqrt C++实现开方

本文链接：https://blog.csdn.net/ARTELE/article/details/109635244

版权

计算机与数学专栏收录该内容

9 篇文章 5 订阅

订阅专栏

0.简介

开根号可以用之前的pow方法计算，当然，还有一些更便捷的方法。

1.二分法

最常见的开根号方法就是二分法，计算 $\small \sqrt{n}$ ，另 $\small x*x-n<10^{-a}$ ，当n的值大于1的时候，可以让l = 0,r = n，当n<1的时候另l=0,r=1，然后利用如下代码计算出结果。

long double l = 0, r = (n>=1 ? n:1), mid=0;
while (fabs(mid * mid - n)>1e-8)//fabs(l-r)<1e-8
{
	mid = (l + r) / 2;
	if (mid * mid > n) r = mid;
	else l = mid;
}

2.牛顿迭代法

牛顿迭代法通过 $\small y=x^{2}-n$ ，来计算 $\small \sqrt{n}$ 。

首先找到曲线y上一点，一般取 $\small x=n$ ，此时得到点 $\small A(n^{2}-n,n)$ ，在A点做切线，与x轴交于点B，这一步主要就是计算从 $\small x=n$ 开始，然后下一步A点应该怎么走能接近y=0的位置，从上图中可以看出A点要向B方向走才可以，这里为了方便计算，索性直接取下一个A点以 $\small x=B_{x}$ ，如下图。

这样一来，B点就逐渐接近y=0的点，经过多次迭代，使得 $\small B_{x}^{2}-n<10^{-a}$ ，其中 $\small B_{x}$ 表示B的x轴坐标值，a是大于1的比较大的数字，当B点很接近y=0点的时候则找到了函数的解。牛顿迭代法在计算根号的时候，初始点要从函数向下凸的位置开始，因为每次算的切线与x轴交点都是在解的同一侧，只能是越来越逼近，不会导致计算结果不收敛。实际实现如下。

long double x = n,y = -1,k,b;
while (fabs(x*x-n)>1e-6)
{
	y = x * x - n;
	k = 2 * x;
	b = y - k * x;
	x = -(b / k);
}

3.问题

实际上，上面的主要思想领悟了，这个算法实现出来就好了，但是我实际上遇到了问题，当计算9000000000000000000000000000000000000000000000000000000000000000000.0的时候，算法竟然失效了，一直在死循环，经过调试我发现，上述两个方法，在fabs(x*x-n)<1e-6这里出现了问题，当数字很大的时候，浮点的底数计算多少会有一点不准确，一丁点误差，但是这一点误差由于有比较大的指数，所以误差实际会很大很大，导致fabs(x*x-n)<1e-6始终成立。并且在计算0.00000000000000000000000000000000000000000000000000000009的时候，由于数字过小，fabs(x*x-n)<1e-6不成立，在第一次无法进入进入循环，针对数字很大和很小的情况，我将要计算的数字都缩放到0.1-1.x之间，实际是(0.1,9.99999...)之间，这样就不会出现计算很大很小的数字，缩放的倍数都是10倍为基准，所以记录乘或者除了几个10，在最终结果上放大回来就可以了。所以有了如下的修改版。

double mysqrt(long double n)
{
	int tens_count = 0;
	double ten = 1;
	if (n > 1)
	{
		while(n > 1)
		{
			n /= 10;
			tens_count++;
		}
		//对于10倍数的奇偶个数要处理好
		if (tens_count % 2 == 1)
		{
			n *= 10;
			tens_count--;
		}
		while (tens_count > 0)
		{
			ten *= 10;
			tens_count -= 2;
		}
	}
	else if (n < 1)
	{
		while (n < 1)
		{
			n *= 10;
			tens_count++;
		}
		if (tens_count % 2 == 1)
		{
			n /= 10;
			tens_count--;
		}
		while (tens_count > 0)
		{
			ten /= 10;
			tens_count -= 2;
		}
	}
	//long double x = n,y = -1,k,b,t;
	//while (fabs(x*x-n)>1e-6)
	//{
	//	t = y;
	//	y = x * x - n;
	//	k = 2 * x;
	//	b = y - k * x;
	//	x = -(b / k);
	//}
	//return x*ten;
	long double l = 0, r = (n>=1 ? n:1), mid=0;
	while (fabs(mid * mid - n)>1e-8)
	{
		mid = (l + r) / 2;
		if (mid * mid > n) r = mid;
		else l = mid;
	}
	return l*ten;
}