平方根的C语言实现(二) —— 手算平方根的原理

1. 版权申明：本文为博主窗户(Colin Cai)原创，欢迎转帖。如要转贴，必须注明原文网址

3. 　　http://www.cnblogs.com/Colin-Cai/p/7220506.html

5. 　　作者：窗户

7. 　　QQ：6679072

9. 　　E-mail：6679072@qq.com

　　一个函数从数学上来说可以有无数个函数列收敛于这个函数，那么程序逼近实现来说可以有无数种算法，平方根自然也不例外。

　　不知道有多少人还记得手算平方根，那是满足每次在结果上添加一位，也就是按位逼近运算结果的唯一算法。至于数学上如何证明这个唯一性我就不说了，数学证明不会有那么多人有兴趣。按位逼近更加适合手算，举个大家更熟悉的例子，那就是手算除法。我这里就采用按位逼近的手算方法。

　　要说手算平方根，原理其实非常简单，

　　一是平方根函数是严格单调增函数，

　　二就是以下这个恒等式满足

　　(a*N+b)2 ≡ (a*N)2 + 2*a*b*N + b2

　　　　　　≡ (a*N)2 + b * ((a*N) * 2 + b)

　　我们实例操作一次平方根笔算，来解释一下。

　　我们来求5499025的平方根。

　　先将5499025两位两位从低往高排，为

　　5 49 90 25

　　2*2<5<3*3

　　所以最高位为2，

　　然后我们再来看549的平方根，

　　我们假设549的平方根的整数部分是2*10+b，则根据之前的恒等式，N在这里等于10，a在这里等于2，有

　　549 >=(2*10)2 + b * ((2*10) * 2 + b)

　　整理一下，149 >= b * (40 + b) 

　　3 * 43 < 149 < 4 * 44

　　所以b=3，

　　549的平方根整数部分是23，

　　再假设54990的平方根整数部分为23*10+b,

　　则

　　54990 >= (23*10)2 + b * ((23*10) * 2 + b)

　　整理一下，2090 >= b * (460 + b),

　　464 * 4  < 2090 < 465 * 5

　　所以b=4,

　　54990的平方根整数部分为234，

　　最后再来看5499025的平方根的整数部分，假设为234 * 10 + b,

　　则

　　5499025 >=  (234*10)2 + b * ((234*10) * 2 + b)

　　整理一下， 23425 >= b * (4680 + b)

　　而5 * 4685 = 23425, 等式成立，

　　所以最终我们要求的平方根是2345。当然，小数位其实一样可以用这种方法继续算下去。

　　手算平方根就是如上这样从高位一步步往地位推的过程，写成式子的形式大致如下：

　　 2    3     4     5

　　-------------------

　  | 5   49   90   25

  2 | 4

　　-------------------

      | 1  49

 43 | 1  29　　　　　　——当前算出了2，2*10*2 = 40

        -------------------

         |  20   90

  464 |  18   56　　　　——当前算出了23，23*10*2 = 460

        -------------------

           |  2   34   25

 4685  |  2   34   25　　——当前算出了234，234*10*2 = 4680

        -------------------

                            0

　　当然，如何写不重要，知道过程便可以继续我们的这个设计。接下去我们要去利用之前的这个算法，改装一下，来进行二进制的开平方。

　　二进制的每一位不是1就是0，这样在每次往前推一位的时候就相对简单。

　　举个例子，我们来算121的平方根，也就是二进制下1111001的平方根。

　　两位两位从低位往高位排

       1   0    1     1

　　------------------

     | 1  11  10  01

  1 | 1

　　------------------

　　  | 11

  100 |   0                  ——当前上面算出了1，1右移动两位为100

　　------------------

　　  | 11  10

1001 | 10  01           ——当前上面算出了10，1右移动两位为1000

　　------------------

            | 1 01  01

 10101 | 1 01 01     ——当前上面算出了101，1右移动两位为10100

　　------------------

                        0

每往右边推1位，下面的除数就是上面当前算出来的二进制的数右移两位再加1或者加0

之后，我们就可以用构建利用此算法的平方根了。

转载于:https://www.cnblogs.com/Colin-Cai/p/7220506.html

整数开平方算法:

本算法只采用移位、加减法、判断和循环实现，因为它不需要浮点运算，也不需要乘除运算，因此可以很方便地运用到各种芯片上去。

我们先来看看10进制下是如何手工计算开方的。
先看下面两个算式，
x = 10*p + q  (1)
公式(1)左右平方之后得：
x^2 = 100*p^2 + 20pq + q^2 (2)
现在假设我们知道x^2和p，希望求出q来，求出了q也就求出了x^2的开方x了。
我们把公式(2)改写为如下格式：
q = (x^2 - 100*p^2)/(20*p+q) (3)
这个算式左右都有q，因此无法直接计算出q来，因此手工的开方算法和手工除法算法一样有一步需要猜值。

我们来一个手工计算的例子：计算1234567890的开方

首先我们把这个数两位两位一组分开，计算出最高位为3。也就是(3)中的p，最下面一行的334为余数，也就是公式(3)中的(x^2 - 100*p^2)近似值

       3
    ---------------
    | 12 34 56 78 90
       9
    ---------------
    |  3 34
下面我们要找到一个0-9的数q使它最接近满足公式(3)。我们先把p乘以20写在334左

后面没有了，因为需要付费才能看2023.12.6 

文章知识点与官方知识档案匹配，可进一步学习相关知识
算法技能树首页概览55984 人正在系统学习中
最低0.47元
————————————————
版权声明：本文为CSDN博主「nhconch」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/nhconch/article/details/7402046

 一起装逼！开平方的七种算法

sqrt()函数，是绝大部分语言支持的常用函数，它实现的是开方运算；开方运算最早是在我国魏晋时数学家刘徽所著的《九章算术》被提及。今天写了几个函数加上国外大神的几个神级程序带大家领略sqrt的神奇之处。

1、古人算法（暴力法）

原理：从0开始0.00001,000002...一个一个试，直到找到x的平方根，代码如下：

1. public class APIsqrt {

3.     static double baoliSqrt(double x) {

5.         final double _JINGDU = 1e-6;

6.         double i;

7.         for (i = 0; Math.abs(x - i * i) > _JINGDU; i += _JINGDU)

8.             ;

9.         return i;

10.     }

12.     public static void main(String[] args) {

13.         double x = 3;

14.         double root = baoliSqrt(x);

15.         System.out.println(root);

16.     }

测试结果：

1.7320509999476947

2、牛顿迭代法

计算机科班出身的童鞋可能首先会想到的是《数值分析》中的牛顿迭代法求平方根。原理是：随意选一个数比如说8（牛顿迭代原理中说这个是初始条件，原则上可以是任意数，查看公式中，它对应Xn给定的初始值，因为7种算法都是以3为例求开平方的，下面的过程需要再检验，是否除以2，否则明显不对，下一篇文中和C语言代码中都明确是除以2），要求根号3，我们可以这么算：

(8 + 3/8) = 4.1875

(4.1875 + 3/4.1875) = 2.4519

(2.4519 + 3/2.4519) = 1.837

(1.837 + 3/1.837) = 1.735

上面的等式，明显是不相等的，但是等式左边除以2，就相等了，可是为什么算3的开方，要用到8呢？2023.12.6

做了4步基本算出了近似值了，这种迭代的方式就是传说中的牛顿迭代法了，代码如下：

1. public class APIsqrt {

3.     static double newtonSqrt(double x) {

5.         if (x < 0) {

6.             System.out.println("负数没事开什么方");

7.             return -1;

8.         }

9.         if (x == 0)

10.             return 0;

12.         double _avg = x;

13.         double last_avg = Double.MAX_VALUE;

14.         final double _JINGDU = 1e-6;

16.         while (Math.abs(_avg - last_avg) > _JINGDU) {

17.             last_avg = _avg;

18.             _avg = (_avg + x / _avg) / 2;

19.         }

20.         return _avg;

21.     }

23.     public static void main(String[] args) {

24.         double x = 3;

25.         double root = newtonSqrt(x);

26.         System.out.println(root);

27.     }

28. }

测试结果：

1.7320508075688772

3、暴力-牛顿综合法

原理：还是以根号3为例，先用暴力法讲根号3逼近到1.7，然后再利用上述的牛顿迭代法。虽然没有用牛顿迭代好，但是也为我们提供一种思路。代码如下：

1. public class APIsqrt {

3.     static double baoliAndNewTonSqrt(double x) {

5.         if (x < 0) {

6.             System.out.println("负数没事开什么方");

7.             return -1;

8.         }

9.         if (x == 0)

10.             return 0;

12.         double i = 0;

13.         double _avg;

14.         double last_avg = Double.MAX_VALUE;

15.         for (i = 0; i*i < x; i += 0.1);

16.         _avg = i;

18.         final double _JINGDU = 1e-6;

20.         while (Math.abs(_avg - last_avg) > _JINGDU) {

21.             last_avg = _avg;

22.             _avg = (_avg + x / _avg) / 2;

23.         }

24.         return _avg;

25.     }

27.     public static void main(String[] args) {

28.         double x = 3;

29.         double root = baoliAndNewTonSqrt(x);

30.         System.out.println(root);

31.     }

32. }

测试结果：

1.7320508075689423

4、二分开方法

原理：还是以3举例：

(0+3)/2 = 1.5, 1.5^2 = 2.25, 2.25 < 3;

(1.5+3)/2 = 2.25, 2.25^2 = 5.0625, 5.0625 > 3;

(1.5+2.25)/2 = 1.875, 1.875^2 = 3.515625; 3.515625>3;

直到前后两次平均值只差小于自定义精度为止，代码如下：

1. public class APIsqrt {

3.     static double erfenSqrt(double x) {

5.         if (x < 0) {

6.             System.out.println("负数没事开什么方");

7.             return -1;

8.         }

9.         if (x == 0)

10.             return 0;

12.         final double _JINGDU = 1e-6;

13.         double _low = 0;

14.         double _high = x;

15.         double _mid = Double.MAX_VALUE;

16.         double last_mid = Double.MIN_VALUE;

18.         while (Math.abs(_mid - last_mid) > _JINGDU) {

20.             last_mid = _mid;

21.             _mid = (_low + _high) / 2;

22.             if (_mid * _mid > x)

23.                 _high = _mid;

24.             if (_mid * _mid < x)

25.                 _low = _mid;

27.         }

28.         return _mid;

30.     }

32.     public static void main(String[] args) {

33.         double x = 3;

34.         double root = erfenSqrt(x);

35.         System.out.println(root);

36.     }

37. }

测试结果：

1.732051134109497

5、计算 (int)(sqrt(x))算法

PS:此算法非博主所写

原理：空间换时间，细节请大家自行探究，代码如下：

1. public class APIsqrt2 {

2.     final static int[] table = { 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53,

3.             55, 57, 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84,

4.             86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104,

5.             106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 119,

6.             120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,

7.             133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144,

8.             145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155,

9.             156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166,

10.             167, 167, 168, 169, 170, 170, 171, 172, 173, 173, 174, 175, 176,

11.             176, 177, 178, 178, 179, 180, 181, 181, 182, 183, 183, 184, 185,

12.             185, 186, 187, 187, 188, 189, 189, 190, 191, 192, 192, 193, 193,

13.             194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, 202,

14.             203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210,

15.             211, 211, 212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218,

16.             218, 219, 219, 220, 221, 221, 222, 222, 223, 224, 224, 225, 225,

17.             226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 231, 232, 232,

18.             233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, 239, 240,

19.             240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246,

20.             247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,

21.             253, 254, 254, 255 };

23.     /**

24.      * A faster replacement for (int)(java.lang.Math.sqrt(x)). Completely

25.      * accurate for x < 2147483648 (i.e. 2^31)...

26.      */

27.     static int sqrt(int x) {

28.         int xn;

30.         if (x >= 0x10000) {

31.             if (x >= 0x1000000) {

32.                 if (x >= 0x10000000) {

33.                     if (x >= 0x40000000) {

34.                         xn = table[x >> 24] << 8;

35.                     } else {

36.                         xn = table[x >> 22] << 7;

37.                     }

38.                 } else {

39.                     if (x >= 0x4000000) {

40.                         xn = table[x >> 20] << 6;

41.                     } else {

42.                         xn = table[x >> 18] << 5;

43.                     }

44.                 }

46.                 xn = (xn + 1 + (x / xn)) >> 1;

47.                 xn = (xn + 1 + (x / xn)) >> 1;

48.                 return ((xn * xn) > x) ? --xn : xn;

49.             } else {

50.                 if (x >= 0x100000) {

51.                     if (x >= 0x400000) {

52.                         xn = table[x >> 16] << 4;

53.                     } else {

54.                         xn = table[x >> 14] << 3;

55.                     }

56.                 } else {

57.                     if (x >= 0x40000) {

58.                         xn = table[x >> 12] << 2;

59.                     } else {

60.                         xn = table[x >> 10] << 1;

61.                     }

62.                 }

64.                 xn = (xn + 1 + (x / xn)) >> 1;

66.                 return ((xn * xn) > x) ? --xn : xn;

67.             }

68.         } else {

69.             if (x >= 0x100) {

70.                 if (x >= 0x1000) {

71.                     if (x >= 0x4000) {

72.                         xn = (table[x >> 8]) + 1;

73.                     } else {

74.                         xn = (table[x >> 6] >> 1) + 1;

75.                     }

76.                 } else {

77.                     if (x >= 0x400) {

78.                         xn = (table[x >> 4] >> 2) + 1;

79.                     } else {

80.                         xn = (table[x >> 2] >> 3) + 1;

81.                     }

82.                 }

84.                 return ((xn * xn) > x) ? --xn : xn;

85.             } else {

86.                 if (x >= 0) {

87.                     return table[x] >> 4;

88.                 }

89.             }

90.         }

92.         return -1;

93.     }

94.     public static void main(String[] args){

95.         System.out.println(sqrt(65));

97.     }

98. }

测试结果：8

6、最快的sqrt算法

PS:此算法非博主所写

这个算法很有名，大家可能也见过，作者是开发游戏的，图形算法中经常用到sqrt，作者才写了一个神级算法，和他那神秘的0x5f3759df，代码如下

1. #include <math.h>

2. float InvSqrt(float x)

3. {

4.  float xhalf = 0.5f*x;

5.  int i = *(int*)&x; // get bits for floating VALUE

6.  i = 0x5f375a86- (i>>1); // gives initial guess y0

7.  x = *(float*)&i; // convert bits BACK to float

8.  x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy

9.  return x;

10. }

12. int main()

13. {

14.   printf("%lf",1/InvSqrt(3));

16.    return 0;

17. }

测试结果：

感兴趣的朋友可以参考http://wenku.baidu.com/view/a0174fa20029bd64783e2cc0.html  是作者解释这个算法的14页论文《Fast Inverse Square Root》

7、一个与算法6相似的算法

PS:此算法非博主所写

代码如下：

1. #include <math.h>

2. float SquareRootFloat(float number) {

3.     long i;

4.     float x, y;

5.     const float f = 1.5F;

7.     x = number * 0.5F;

8.     y  = number;

9.     i  = * ( long * ) &y;

10.     i  = 0x5f3759df - ( i >> 1 );

11.     y  = * ( float * ) &i;

12.     y  = y * ( f - ( x * y * y ) );

13.     y  = y * ( f - ( x * y * y ) );

14.     return number * y;

15. }

17. int main()

18. {

19.   printf("%f",SquareRootFloat(3));

21.    return 0;

22. }

测试结果：

-----------------------

公众号：五分钟学算法（ID：CXYxiaowu）

博客：www.cxyxiaowu.com

知乎：程序员吴师兄

小程序：点击前往五分钟学算法小程序

一个正在学习算法的人，致力于将算法讲清楚！

一起装逼！开平方的七种算法-CSDN博客