2 # -*- coding: ISO-8859-1 -*-
5 # Copyright (C) 2002-2004 Jörg Lehmann <joergl@users.sourceforge.net>
6 # Copyright (C) 2002-2004 André Wobst <wobsta@users.sourceforge.net>
8 # This file is part of PyX (http://pyx.sourceforge.net/).
10 # PyX is free software; you can redistribute it and/or modify
11 # it under the terms of the GNU General Public License as published by
12 # the Free Software Foundation; either version 2 of the License, or
13 # (at your option) any later version.
15 # PyX is distributed in the hope that it will be useful,
16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 # GNU General Public License for more details.
20 # You should have received a copy of the GNU General Public License
21 # along with PyX; if not, write to the Free Software
22 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
25 # XXX fallback for Numeric (eigenvalue computation) to be implemented along
26 # know algorithms (like from numerical recipes)
28 import Numeric
, LinearAlgebra
30 def realpolyroots(coeffs
, epsilon
=1e-5, polish
=0):
32 """returns the roots of a polynom with given coefficients
34 This helper routine uses the package Numeric to find the roots
35 of the polynomial with coefficients given in coeffs:
36 0 = \sum_{i=0}^N x^i coeffs[i]
37 The solution is found via an equivalent eigenvalue problem
43 roots
= realpolyroots(coeffs
[:-1], epsilon
=epsilon
)
45 # divide the coefficients by the maximal order
46 # this is to be done in the matrix, anyhow and
47 # makes checking more stable
48 coeffs
= [coeff
/coeffs
[-1] for coeff
in coeffs
]
51 # build the Matrix of the polynomial problem
52 mat
= Numeric
.zeros((N
, N
), Numeric
.Float
)
56 mat
[0][i
] = -coeffs
[N
-1-i
]
57 # find the eigenvalues of the matrix (== the roots of the polynomial)
58 roots
= [complex(root
) for root
in LinearAlgebra
.eigenvalues(mat
)]
59 # take only the real roots
60 roots
= [root
.real
for root
in roots
if -epsilon
< root
.imag
< epsilon
]
62 # polish the roots with Newton-Raphson
64 def polish(root
, epsilon
):
66 while abs(polynom
) > epsilon
:
67 polynom
= coeffs
[N
]*root
+ coeffs
[N
-1]
69 for i
in range(N
-2,-1,-1):
70 polynom
= polynom
*root
+ coeffs
[i
]
71 poprime
= poprime
*root
+ coeffs
[i
+1]*(i
+1)
72 root
-= polynom
/ poprime
75 roots
= [polish(root
, epsilon
) for root
in roots
]
77 # # check if the roots are really roots!
79 # polynom = coeffs[N]*root + coeffs[N-1]
80 # for i in range(N-2,-1,-1):
81 # polynom = polynom*root + coeffs[i]
82 # if abs(polynom) > epsilon:
83 # raise Exception("value %f instead of 0" % polynom)