求平方根
数学里面有很多操作是互逆的,正向操作简单,但是逆向操作很复杂。平方的逆操作是开平方根,这个开方操作计算起来并不容易。费了我们先人很多功夫。
在中学的课本中,会讲到竖式计算法。它的基本思路是这样的:
比如1156是四位数,不难看出它的平方根是个两位数,且十位上的数字是3.于是问题就转换成:怎样求出它的个位数a?根据两数和的平方公式,可以得到
1156=(30+a)^2=30^2+2*30a+a^2,所以 1156-30^2=2×30a+a^2,即 256=(30×2+a)a,
这就是说, a是这样一个数,它与30*2的和,再乘以它本身,等于256.
为便于求得a,可用下面的方法进行计算:
3是平方根的十位数.将1156-30^2,得到余数 256,把这个余数试除以30×2,得4(如果未除尽则取整数位).由于4与30×2的和64,与4的积正好等于256,开方正好开尽,4就是所求的个位数a.于是得到 1156=34^2。对于不能开尽的数,通过补充00来一步步接近。
这个最早出现在中国古代数学著作《九章算术》(约2000年前)的“少广”章。这是中国古代数学的杰出成就,系统地介绍了中国古代的数学体系。但是该书的整个体系是偏向于实用的,没有原理证明,与以《几何原本》为代表的公理逻辑体系形成鲜明对比,所以流于表面解决问题,后劲不足,最后没有发展出现代科学。这是我们在了解中国古代科技成就的时候需要明了的。
上面的开平方运算非常繁,我们现在很少用到。感谢古人的不辞辛劳,通过手工计算给了我们一个庞大的开方表。技术的进步都是一点一滴的,通过无数代人接力,历经千年万年,人类整体智慧就会到达一个令人惊叹的地步。所以,我们永远不要忘记先人的呕心沥血。
在计算机算法里,我们其实不是用的上述办法,那实在太麻烦了,我们会用几种别的方法,比如二分法,如牛顿逼近法(牛顿真的是大神,学着学着就会碰到他)。
二分法是一种渐渐逼近的办法,比如要求10的平方根,观察一下不难得知这个数肯定在0到10之间,我们取中间值5来算一下,为25,比10大,所以这个数一定比5小而比0大,我们再次取中间值2.5,这样一直算下去逐步逼近真数。
这个办法很笨很笨,通过很多次循环一点点接近。这里再次看到了计算机的笨和快,因为计算机跑得快,所以笨方法还很有效。
我们用下面程序实现上面的算法:
print ("Square root calculator")
x = 10
a=0
b=x
c=(a+b)/2
while c**2-x>0.0001 or c**2-x<-0.0001:
if c**2>x:
a=a
b=c
if c**2<x:
a=c
b=b
c=(a+b)/2
print (c)
稍微解释一下,x为给定的原值,a为平方根的下限,b为平方根的上限,c为猜测的值,初始为中间值c=(a+b)/2,我们用这个值逼近真值。
While 的循环条件为精度要求。如果还没有达到精度要求,就判断试值的平方大于小时小于原值,如果大于,说明这个试值猜测得大了,那么下一个猜测得值应该在下限与当前试值中间,所以我么你重新划定上下限,这就是a=a,b=c 两句的含义,反之,我们用a=c,b=b划定上下限。之后再用新的上下限的中间值当成下一个试值。一直这么循环下去,最后得到符合精度要求的结果值。
要记住一点,其实我们得不到真值的,得到的永远是近似值。
自然,我们审视上面的程序的时候,也会发现问题。有了以前的基础,估计你也能知道几个明显的问题,第一个就是x=10,这个写死了,不好,应该是从外部接受一个值。除了这个呢?还有没有问题?有的,你看到了,这个程序没有处理负数。
但是,这个程序最大的问题还不在这里,而是有一个大bug。如果给定的这个值小于1,这个程序会陷于无穷循环中,结束不了了。这是因为x<1的时候,x的平方根是大于x的,落在的0到x区间之外。所以我们要改进一下我们的算法,判断一下x,如果x<1,就把上下限设定为x到1之间。修改后的程序如下:
print ("Square root calculator")
x = float(input("enter one number:"))
if x>1:
a=0
b=x
elif x>=0:
a=x
b=1
else:
print("error: negative number.")
exit()
c=(a+b)/2
while c**2-x>0.0001 or c**2-x<-0.0001:
if c**2>x:
a=a
b=c
if c**2<x:
a=c
b=b
c=(a+b)/2
print (c)
修改后的程序,原值由外部输入,判断了大于1大于0小于0三种情况。其中用到了exit()这个Python提供的系统函数,用于退出程序执行。这个办法不是很好,以后我们会讲到怎么更好地处理。
有些用心的人会想知道我们究竟用了多少次循环才得到了这个逼近的数,这里与精度与关,你可以在循环代码里面加上一个计数器,看看结果。
我自己试验的结果,四位精度下,计算10的平方根,循环18次,十位精度下是35次。
有一些人特别好奇,究竟多少位精度合适呢?我们计算圆周率Pi的时候经常也会有这样的疑惑。这个与对计算结果的要求有关了,反正呐,我看过一个说明,如果Pi取小数点后35位,那么计算出来的太阳系的周长的值的误差只有一个原子那么一点。
有时候会觉得数学真的很神奇,让我们坐在桌子前用一支铅笔居然能算出浩瀚的宇宙。
还有一些严谨的人在使用算法的时候老是心里犯嘀咕:怎么知道这个算法本身是正确的呢?一般的做法是多测试几个数。但是数是无穷的,测试再多也不能表明算法是正确的,算法需要证明才能成为真理。大家放心,这些算法都是得到了证明的。我这里不给出证明,这需要用到高等数学知识,你们可以从下面的表达感受一下:
求解 f(x)=x2−n 的零点,在 [0,n] 的区间上,f(x)=x2−n 自然单调递增,可证二分法必然收敛。
这已经到了初等数学知识的边界了(中国的数学到此为止)。我们需要下车,再换乘高等数学(近代的数学)的列车,也有一批人会永久下车。
在计算平方根的几种算法中,还有一种用得比较广泛,就是牛顿逼近法,这个算法更加高效。
它用到了迭代,公式为Xi+1=Xi - (Xi2 - n) / (2Xi) =(Xi + n/Xi) / 2。有了这个迭代公式,程序就好编写了:
n=float(input("enter one number:"))
if n>=0:
y=n/2
while y**2-n>0.000001 or y**2-n<-0.000001:
y=(y+n/y)/2
print(y)
else:
print("Error: negative number")
解释一下,初始值取得是n/2,不是取的n,因为n的平方根不会超过n/2。并且这个方法不需要单独考虑小于1的情况,它自己通过迭代会收敛到真值。
程序的主体就是那个while循环实现迭代公式,y=(y+n/y)/2。
那么我们究竟用了多少次迭代呢?我自己试验的结果,计算10的平方根,十位精度下是5次。
看到牛顿逼近算法的力量了吧?程序更加简洁,效率更高。
这个算法其实是一个推论,先就有数学上的严格的表达和证明。它使用函数 f(x) 的泰勒级数的前面几项来寻找方程 f(x)=0 的根。该法效率高,证明它是平方收敛的,该法应用广,并不限于求解实数的平方根,相反求解实数的平方根只是其一个具体的应用而已。
对于求解实数平方根的函数 f(x)=x2−n, 自然其根的迭代公式为:
xn+1=xn−f(xn)/f′(xn)=xn−f(xn)/2xn
我们不深究理论,但是从这里我们可以感受到牛顿的伟大和近代数学的魅力,它帮助我们精确分析时空和运动。有一首诗写道:
自然隐藏在幽暗之中,
上帝说:
让牛顿降生吧。
于是有了光明。
求阶乘的平方根
这个题目有点趣味,组合了前面两个题目。
虽然内心有点疑惑,但是同古人同样不辞辛劳的你把上面两个题目的程序拼在一起:
n=int(input("enter a number:"))
result = 1
while n>1 :
result = result * n
n = n - 1
print (result)
n=result
if n>=0:
y=n/2
while y**2-n>0.0001 or y**2-n<-0.0001:
y=(y+n/y)/2
print(y)
else:
print("Error: negative number")
运行结果,输入10,结果输出
3628800
1904.9409439665096
没有错。
虽然如此,你心里还是觉得不对劲:计算机不应该这么笨吧?它不能把以前写的程序用某种方式复用一下吗?
你的感觉是对的,有办法的。而且这个办法在程序出现的时候就有了,这就是Ada提出来的子过程的概念。我们可以把一段程序封装成一个单独的子过程,接受输入,计算处理后,输出结果,这个子过程可以被别的程序反复调用。
Python里面是这么定义子过程的:
def functionname(p1,p2,...):
解释一下,def关键字代表这里要定义一个子过程,后面是子过程的名字,接下来的括号是参数(输入值,数学上叫自变量)说明,后面过程代码段。用return命令输出结果。
按照这个定义,我们先改写阶乘程序(factorial.py):
def factorial(n):
result = 1
while n>1 :
result = result * n
n = n - 1
return result
print (factorial(10))
程序里先用def定义一个叫做factorial的子过程,计算后return result。
使用这个子过程的代码很简单,直接调用factorial(10)。由于子程序是给人调用的,所以最后提供出来的时候要把print()语句去掉,否则会网屏幕上输出结果。
同样照此办理,我们改写牛顿逼近算法程序为(sqrt.py):
def sqrt(n):
if n>=0:
y=n/2
while y**2-n>0.000001 or y**2-n<-0.000001:
y=(y+n/y)/2
return y
else:
return -1
print (sqrt(10))
测试运行,表明新创建的sqrt()子过程也是对的。
注意这个子过程的实现,else里面是return -1。以前我们对负数是屏幕输出错误。但是子过程的目的是供别的程序调用,所以子过程自己不要这么处理,返回一个永远不可能的值(-1),错误处理交给调用者。
以前调用exit退出程序,那个办法不好,到现在我们看到了新的处理方法。
好,我们手头现在有了两个子过程,我们这个题目就可以直接使用了。
一开头,你大概会这么编写:
n=int(input("enter a number:"))
x = factorial(n)
y = sqrt(x)
print(y)
运行一下,出错了:
NameError: name 'factorial' is not defined
计算机不知道这个factorial()子过程,同理也不知道sqrt()子过程。
可是,我们明明定义了啊。问题出在程序空间,本题的程序不能识别以前的题目编写的程序,为了识别,重用以前的程序,必须import一下。程序如下:
import factorial
import sqrt
n=int(input("enter a number:"))
x = factorial.factorial(n)
y = sqrt.sqrt(x)
print(y)
程序开始的两行import了子过程程序文件,factorial.py和sqrt.py,这样把以前编写的程序引入了。调用的时候不直接写子过程名字,而是在前面加了程序文件名,如factorial.factorial(),表示使用factorial这个程序里的factorial子过程。
更加好的做法,是把相关的子过程写在一个程序文件中,作为一个函数包供大家调用。
我们可以这么提供一个maths程序(maths.py):
#this is a mathematics library
#square root function
def sqrt(n):
if n>=0:
y=n/2
while y**2-n>0.000001 or y**2-n<-0.000001:
y=(y+n/y)/2
return y
else:
return -1
#factorial function
def factorial(n):
result = 1
while n>1 :
result = result * n
n = n - 1
return result
利用这个maths数学包,我们把本题的程序改写为:
import maths
n=int(input("enter a number:"))
x = maths.factorial(n)
y = maths.sqrt(x)
print(y)
现在我们看到了,如何通过子过程复用和简化问题。实际上,一个大一点的程序就是仔细分解子程序,最后通过某种结构将子程序组合在一起,积木式进行程序构建。最微型的原子子过程,组合成功能大一点的次小子过程,再进一步组合成大一些的子过程,这么一级级分解组合。看起来程序像一台机器,有很多组成部分,每一部分又分成许多部件,每一个部件又由许多零件组成。这样就可以由专门的人负责专门的子过程了,大家的工作通过调用的参数和输出关联在一起,这种组织方式叫“软件工程”。