python计算方程式根的方法

''' roots = polyRoots(a).
 Uses Laguerre's method to compute all the roots of
 a[0] + a[1]*x + a[2]*x^2 +...+ a[n]*x^n = 0.
 The roots are returned in the array 'roots',
''' 
from evalPoly import *
from numpy import zeros,complex
from cmath import sqrt
from random import random
def polyRoots(a,tol=1.0e-12):
 def laguerre(a,tol):
 x = random()
 # Starting value (random number)
 n = len(a) - 1
 for i in range(30):
 p,dp,ddp = evalPoly(a,x)
 if abs(p) < tol: return x
 g = dp/p
 h = g*g - ddp/p
 f = sqrt((n - 1)*(n*h - g*g))
 if abs(g + f) > abs(g - f): dx = n/(g + f)
 else: dx = n/(g - f)
 x = x - dx
 if abs(dx) < tol: return x
 print 'Too many iterations'
 def deflPoly(a,root): # Deflates a polynomial
 n = len(a)-1
 b = [(0.0 + 0.0j)]*n
 b[n-1] = a[n]
 for i in range(n-2,-1,-1):
 b[i] = a[i+1] + root*b[i+1]
 return b
 n = len(a) - 1
 roots = zeros((n),dtype=complex)
 for i in range(n):
 x = laguerre(a,tol)
 if abs(x.imag) < tol: x = x.real
 roots[i] = x
 a = deflPoly(a,x)
 return roots
 raw_input("
Press return to exit")