# -*- coding: ISO-8859-1 -*-
#
#
# Copyright (C) 2005-2006 Michael Schindler <m-schindler@users.sourceforge.net>
# Copyright (C) 2006 Andr Wobst <wobsta@users.sourceforge.net>
#
# This file is part of PyX (http://pyx.sourceforge.net/).
#
# PyX is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# PyX is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with PyX; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
import math, types
# try:
# import Numeric, LinearAlgebra
# _has_numeric = 1
# except:
# _has_numeric = 0
def sign(x):
"""sign of x, i.e. +1 or -1; returns 1 for x == 0"""
if x >= 0:
return 1
return -1
def asinh(x):
"""Return the arc hyperbolic sine of x."""
return math.log(x+math.sqrt(x*x+1))
def acosh(x):
"""Return the arc hyperbolic cosine of x."""
return math.log(x+math.sqrt(x*x-1))
def _realroots_quadratic(a1, a0):
"""gives the real roots of x**2 + a1 * x + a0 = 0"""
D = a1*a1 - 4*a0
if D < 0:
return []
SD = math.sqrt(D)
return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)]
def _realroots_cubic(a2, a1, a0):
"""gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0"""
# see http://mathworld.wolfram.com/CubicFormula.html for details
Q = (3*a1 - a2*a2) / 9.0
R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0
D = Q*Q*Q + R*R
if D > 0: # one real and two complex roots
SD = math.sqrt(D)
if R + SD >= 0:
S = (R + SD)**(1/3.0)
else:
S = -(-R - SD)**(1/3.0)
if R - SD >= 0:
T = (R - SD)**(1/3.0)
else:
T = -(SD - R)**(1/3.0)
return [S + T - a2/3.0]
elif D == 0:
if Q == 0: # one real root (R==0)
return [-a2/3.0]
else: # two real roots (R>0, Q<0)
S = -math.sqrt(-Q)
return [2*S - a2/3.0, -S - a2/3.0]
else: # three real roots (Q<0)
SQ = math.sqrt(-Q)
arg = R / (SQ**3)
if arg >= 1:
theta = 0
elif arg <= -1:
theta = math.pi
else:
theta = math.acos(R/(SQ**3))
return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)]
def _realroots_quartic(a3, a2, a1, a0):
"""gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0"""
# see http://mathworld.wolfram.com/QuarticEquation.html for details
ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3)
ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and y*y-4*a0 >= 0]
if not ys:
return []
y1 = min(ys)
if a3*y1-2*a1 < 0:
return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) +
_realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))))
else:
return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) +
_realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))))
def realpolyroots(*cs):
"""returns the roots of a polynom with given coefficients
polynomial with coefficients given in cs:
0 = \sum_i cs[i] * x^(len(cs)-i-1)
"""
if not cs:
return [0]
try:
f = 1.0/cs[0]
cs = [f*c for c in cs[1:]]
except ArithmeticError:
return realpolyroots(*cs[1:])
else:
n = len(cs)
if n == 0:
return []
elif n == 1:
return [-cs[0]]
elif n == 2:
return _realroots_quadratic(*cs)
elif n == 3:
return _realroots_cubic(*cs)
elif n == 4:
return _realroots_quartic(*cs)
else:
raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
# def realpolyroots_eigenvalue(*cs):
# # as realpolyroots but using an equivalent eigenvalue problem
# # (this code is currently used for functional tests only)
# if not _has_numeric:
# raise RuntimeError("realpolyroots_eigenvalue depends on Numeric")
# if not cs:
# return [0]
# try:
# f = 1.0/cs[0]
# cs = [f*c for c in cs[1:]]
# except ArithmeticError:
# return realpolyroots_eigenvalue(*cs[1:])
# else:
# if not cs:
# return []
# n = len(cs)
# a = Numeric.zeros((n, n), Numeric.Float)
# for i in range(n-1):
# a[i+1][i] = 1
# for i in range(n):
# a[0][i] = -cs[i]
# rs = []
# for r in LinearAlgebra.eigenvalues(a):
# if type(r) == types.ComplexType:
# if not r.imag:
# rs.append(r.real)
# else:
# rs.append(r)
# return rs
#
|