2 # -*- coding: ISO-8859-1 -*-
5 # Copyright (C) 2005-2006 Michael Schindler <m-schindler@users.sourceforge.net>
6 # Copyright (C) 2006 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
27 # import Numeric, LinearAlgebra
34 """sign of x, i.e. +1 or -1; returns 1 for x == 0"""
41 """Return the arc hyperbolic sine of x."""
42 return math
.log(x
+math
.sqrt(x
*x
+1))
46 """Return the arc hyperbolic cosine of x."""
47 return math
.log(x
+math
.sqrt(x
*x
-1))
50 def _realroots_quadratic(a1
, a0
):
51 """gives the real roots of x**2 + a1 * x + a0 = 0"""
56 return [0.5 * (-a1
+ SD
), 0.5 * (-a1
- SD
)]
59 def _realroots_cubic(a2
, a1
, a0
):
60 """gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0"""
61 # see http://mathworld.wolfram.com/CubicFormula.html for details
63 Q
= (3*a1
- a2
*a2
) / 9.0
64 R
= (9*a2
*a1
- 27*a0
- 2*a2
*a2
*a2
) / 54.0
67 if D
> 0: # one real and two complex roots
72 S
= -(-R
- SD
)**(1/3.0)
76 T
= -(SD
- R
)**(1/3.0)
77 return [S
+ T
- a2
/3.0]
79 if Q
== 0: # one real root (R==0)
81 else: # two real roots (R>0, Q<0)
83 return [2*S
- a2
/3.0, -S
- a2
/3.0]
84 else: # three real roots (Q<0)
92 theta
= math
.acos(R
/(SQ
**3))
93 return [2 * SQ
* math
.cos((theta
+ 2*2*i
*math
.pi
)/3.0) - a2
/3.0 for i
in range(3)]
96 def _realroots_quartic(a3
, a2
, a1
, a0
):
97 """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0"""
98 # see http://mathworld.wolfram.com/QuarticEquation.html for details
99 ys
= _realroots_cubic(-a2
, a1
*a3
- 4*a0
, 4*a0
*a2
- a1
*a1
- a0
*a3
*a3
)
100 ys
= [y
for y
in ys
if a3
*a3
-4*a2
+4*y
>= 0 and y
*y
-4*a0
>= 0]
105 return (_realroots_quadratic(0.5*(a3
+math
.sqrt(a3
*a3
-4*a2
+4*y1
)), 0.5*(y1
-math
.sqrt(y1
*y1
-4*a0
))) +
106 _realroots_quadratic(0.5*(a3
-math
.sqrt(a3
*a3
-4*a2
+4*y1
)), 0.5*(y1
+math
.sqrt(y1
*y1
-4*a0
))))
108 return (_realroots_quadratic(0.5*(a3
+math
.sqrt(a3
*a3
-4*a2
+4*y1
)), 0.5*(y1
+math
.sqrt(y1
*y1
-4*a0
))) +
109 _realroots_quadratic(0.5*(a3
-math
.sqrt(a3
*a3
-4*a2
+4*y1
)), 0.5*(y1
-math
.sqrt(y1
*y1
-4*a0
))))
112 def realpolyroots(*cs
):
113 """returns the roots of a polynom with given coefficients
115 polynomial with coefficients given in cs:
116 0 = \sum_i cs[i] * x^(len(cs)-i-1)
122 cs
= [f
*c
for c
in cs
[1:]]
123 except ArithmeticError:
124 return realpolyroots(*cs
[1:])
132 return _realroots_quadratic(*cs
)
134 return _realroots_cubic(*cs
)
136 return _realroots_quartic(*cs
)
138 raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
141 # def realpolyroots_eigenvalue(*cs):
142 # # as realpolyroots but using an equivalent eigenvalue problem
143 # # (this code is currently used for functional tests only)
144 # if not _has_numeric:
145 # raise RuntimeError("realpolyroots_eigenvalue depends on Numeric")
150 # cs = [f*c for c in cs[1:]]
151 # except ArithmeticError:
152 # return realpolyroots_eigenvalue(*cs[1:])
157 # a = Numeric.zeros((n, n), Numeric.Float)
158 # for i in range(n-1):
163 # for r in LinearAlgebra.eigenvalues(a):
164 # if type(r) == types.ComplexType: