pyx-checkins

 [PyX-checkins] pyx/pyx helper.py,1.16,1.17 From: Michael Schindler - 2005-09-05 14:04:13 Update of /cvsroot/pyx/pyx/pyx In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv22977 Modified Files: helper.py Log Message: change the coefficients order in realpolyroots Index: helper.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/helper.py,v retrieving revision 1.16 retrieving revision 1.17 diff -C2 -d -r1.16 -r1.17 *** helper.py 1 Sep 2005 08:26:21 -0000 1.16 --- helper.py 5 Sep 2005 14:04:05 -0000 1.17 *************** *** 36,40 **** import Numeric, LinearAlgebra ! def realpolyroots(coeffs, epsilon=1e-5): """returns the roots of a polynom with given coefficients --- 36,40 ---- import Numeric, LinearAlgebra ! def realpolyroots(coeffs, epsilon=1e-5, polish=0): """returns the roots of a polynom with given coefficients *************** *** 42,53 **** This helper routine uses the package Numeric to find the roots of the polynomial with coefficients given in coeffs: ! 0 = \sum_{i=0}^N x^{N-i} coeffs[i] The solution is found via an equivalent eigenvalue problem """ try: ! 1.0 / coeffs[0] except: ! return realpolyroots(coeffs[1:], epsilon=epsilon) else: --- 42,53 ---- This helper routine uses the package Numeric to find the roots of the polynomial with coefficients given in coeffs: ! 0 = \sum_{i=0}^N x^i coeffs[i] The solution is found via an equivalent eigenvalue problem """ try: ! 1.0 / coeffs[-1] except: ! roots = realpolyroots(coeffs[:-1], epsilon=epsilon) else: *************** *** 58,75 **** mat[i+1][i] = 1 for i in range(N): ! mat[0][i] = -coeffs[i+1] / coeffs[0] ! # find the eigenvalues of the matrix (== the zeros of the polynomial) ! zeros = [complex(zero) for zero in LinearAlgebra.eigenvalues(mat)] ! # take only the real zeros ! zeros = [zero.real for zero in zeros if -epsilon < zero.imag < epsilon] ! ## check if the zeros are really zeros! ! #for zero in zeros: ! # p = 0 ! # for i in range(N+1): ! # p += coeffs[i] * zero**(N-i) ! # if abs(p) > epsilon: ! # raise Exception("value %f instead of 0" % p) ! return zeros --- 58,90 ---- mat[i+1][i] = 1 for i in range(N): ! mat[0][i] = -coeffs[N-1-i] / coeffs[N] ! # find the eigenvalues of the matrix (== the roots of the polynomial) ! roots = [complex(root) for root in LinearAlgebra.eigenvalues(mat)] ! # take only the real roots ! roots = [root.real for root in roots if -epsilon < root.imag < epsilon] ! # polish the roots with Newton-Raphson ! if polish: ! def polish(root, epsilon): ! polynom = 2*epsilon ! while abs(polynom) > epsilon: ! polynom = coeffs[N]*root + coeffs[N-1] ! poprime = coeffs[N]*N ! for i in range(N-2,-1,-1): ! polynom = polynom*root + coeffs[i] ! poprime = poprime*root + coeffs[i+1]*(i+1) ! root -= polynom / poprime ! return root ! roots = [polish(root, epsilon) for root in roots] ! ! # # check if the roots are really roots! ! # for root in roots: ! # polynom = coeffs[N]*root + coeffs[N-1] ! # for i in range(N-2,-1,-1): ! # polynom = polynom*root + coeffs[i] ! # if abs(polynom) > epsilon: ! # raise Exception("value %f instead of 0" % polynom) ! ! return roots