pyx/mathutils.py

   1 # -*- coding: ISO-8859-1 -*-
   2 #
   3 #
   4 # Copyright (C) 2005-2006 Michael Schindler <m-schindler@users.sourceforge.net>
   5 # Copyright (C) 2006 André Wobst <wobsta@users.sourceforge.net>
   6 #
   7 # This file is part of PyX (http://pyx.sourceforge.net/).
   8 #
   9 # PyX is free software; you can redistribute it and/or modify
  10 # it under the terms of the GNU General Public License as published by
  11 # the Free Software Foundation; either version 2 of the License, or
  12 # (at your option) any later version.
  13 #
  14 # PyX is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU General Public License for more details.
  18 #
  19 # You should have received a copy of the GNU General Public License
  20 # along with PyX; if not, write to the Free Software
  21 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  22
  23 import math, types
  24
  25 # try:
  26 #     import Numeric, LinearAlgebra
  27 #     _has_numeric = 1
  28 # except:
  29 #     _has_numeric = 0
  30
  31
  32 def sign(x):
  33     """sign of x, i.e. +1 or -1; returns 1 for x == 0"""
  34     if x >= 0:
  35         return 1
  36     return -1
  37
  38
  39 def asinh(x):
  40     """Return the arc hyperbolic sine of x."""
  41     return math.log(x+math.sqrt(x*x+1))
  42
  43
  44 def acosh(x):
  45     """Return the arc hyperbolic cosine of x."""
  46     return math.log(x+math.sqrt(x*x-1))
  47
  48
  49 def _realroots_quadratic(a1, a0):
  50     """gives the real roots of x**2 + a1 * x + a0 = 0"""
  51     D = a1*a1 - 4*a0
  52     if D < 0:
  53         return []
  54     SD = math.sqrt(D)
  55     return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)]
  56
  57
  58 def _realroots_cubic(a2, a1, a0):
  59     """gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0"""
  60     # see http://mathworld.wolfram.com/CubicFormula.html for details
  61
  62     Q = (3*a1 - a2*a2) / 9.0
  63     R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0
  64     D = Q*Q*Q + R*R
  65
  66     if D > 0:         # one real and two complex roots
  67         SD = math.sqrt(D)
  68         if R + SD >= 0:
  69             S = (R + SD)**(1/3.0)
  70         else:
  71             S = -(-R - SD)**(1/3.0)
  72         if R - SD >= 0:
  73             T = (R - SD)**(1/3.0)
  74         else:
  75             T = -(SD - R)**(1/3.0)
  76         return [S + T - a2/3.0]
  77     elif D == 0:
  78         if Q == 0:    # one real root (R==0)
  79             return [-a2/3.0]
  80         else:         # two real roots (R>0, Q<0)
  81             S = -math.sqrt(-Q)
  82             return [2*S - a2/3.0, -S - a2/3.0]
  83     else:             # three real roots (Q<0)
  84         SQ = math.sqrt(-Q)
  85         arg = R / (SQ**3)
  86         if arg >= 1:
  87             theta = 0
  88         elif arg <= -1:
  89             theta = math.pi
  90         else:
  91             theta = math.acos(R/(SQ**3))
  92         return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)]
  93
  94
  95 def _realroots_quartic(a3, a2, a1, a0):
  96     """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0"""
  97     # see http://mathworld.wolfram.com/QuarticEquation.html for details
  98     ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3)
  99     ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and  y*y-4*a0 >= 0]
 100     if not ys:
 101         return []
 102     y1 = min(ys)
 103     if a3*y1-2*a1 < 0:
 104         return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) +
 105                 _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))))
 106     else:
 107         return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) +
 108                 _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))))
 109
 110
 111 def realpolyroots(*cs):
 112     """returns the roots of a polynom with given coefficients
 113
 114     polynomial with coefficients given in cs:
 115       0 = \sum_i  cs[i] * x^(len(cs)-i-1)
 116     """
 117     if not cs:
 118         return [0]
 119     try:
 120         f = 1.0/cs[0]
 121         cs = [f*c for c in cs[1:]]
 122     except ArithmeticError:
 123         return realpolyroots(*cs[1:])
 124     else:
 125         n = len(cs)
 126         if n == 0:
 127             return []
 128         elif n == 1:
 129             return [-cs[0]]
 130         elif n == 2:
 131             return _realroots_quadratic(*cs)
 132         elif n == 3:
 133             return _realroots_cubic(*cs)
 134         elif n == 4:
 135             return _realroots_quartic(*cs)
 136         else:
 137             raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
 138
 139
 140 # def realpolyroots_eigenvalue(*cs):
 141 #     # as realpolyroots but using an equivalent eigenvalue problem
 142 #     # (this code is currently used for functional tests only)
 143 #     if not _has_numeric:
 144 #         raise RuntimeError("realpolyroots_eigenvalue depends on Numeric")
 145 #     if not cs:
 146 #         return [0]
 147 #     try:
 148 #         f = 1.0/cs[0]
 149 #         cs = [f*c for c in cs[1:]]
 150 #     except ArithmeticError:
 151 #         return realpolyroots_eigenvalue(*cs[1:])
 152 #     else:
 153 #         if not cs:
 154 #             return []
 155 #         n = len(cs)
 156 #         a = Numeric.zeros((n, n), Numeric.Float)
 157 #         for i in range(n-1):
 158 #             a[i+1][i] = 1
 159 #         for i in range(n):
 160 #             a[0][i] = -cs[i]
 161 #         rs = []
 162 #         for r in LinearAlgebra.eigenvalues(a):
 163 #             if type(r) == types.ComplexType:
 164 #                 if not r.imag:
 165 #                     rs.append(r.real)
 166 #             else:
 167 #                 rs.append(r)
 168 #         return rs
 169 #