計算機與代數---如何計算sqrt---方法和實作-有解無憂

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;
}