pyx/mathutils.py

   1 #!/usr/bin/env python
   2 # -*- coding: ISO-8859-1 -*-
   3 #
   4 #
   5 # Copyright (C) 2005-2006 Michael Schindler <m-schindler@users.sourceforge.net>
   6 # Copyright (C) 2006 André Wobst <wobsta@users.sourceforge.net>
   7 #
   8 # This file is part of PyX (http://pyx.sourceforge.net/).
   9 #
  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.
  14 #
  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.
  19 #
  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
  23
  24 import math, types
  25
  26 # try:
  27 #     import Numeric, LinearAlgebra
  28 #     _has_numeric = 1
  29 # except:
  30 #     _has_numeric = 0
  31
  32
  33 def sign(x):
  34     """sign of x, i.e. +1 or -1; returns 1 for x == 0"""
  35     if x >= 0:
  36         return 1
  37     return -1
  38
  39
  40 def asinh(x):
  41     """Return the arc hyperbolic sine of x."""
  42     return math.log(x+math.sqrt(x*x+1))
  43
  44
  45 def acosh(x):
  46     """Return the arc hyperbolic cosine of x."""
  47     return math.log(x+math.sqrt(x*x-1))
  48
  49
  50 def _realroots_quadratic(a1, a0):
  51     """gives the real roots of x**2 + a1 * x + a0 = 0"""
  52     D = a1*a1 - 4*a0
  53     if D < 0:
  54         return []
  55     SD = math.sqrt(D)
  56     return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)]
  57
  58
  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
  62
  63     Q = (3*a1 - a2*a2) / 9.0
  64     R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0
  65     D = Q*Q*Q + R*R
  66
  67     if D > 0:         # one real and two complex roots
  68         SD = math.sqrt(D)
  69         if R + SD >= 0:
  70             S = (R + SD)**(1/3.0)
  71         else:
  72             S = -(-R - SD)**(1/3.0)
  73         if R - SD >= 0:
  74             T = (R - SD)**(1/3.0)
  75         else:
  76             T = -(SD - R)**(1/3.0)
  77         return [S + T - a2/3.0]
  78     elif D == 0:
  79         if Q == 0:    # one real root (R==0)
  80             return [-a2/3.0]
  81         else:         # two real roots (R>0, Q<0)
  82             S = -math.sqrt(-Q)
  83             return [2*S - a2/3.0, -S - a2/3.0]
  84     else:             # three real roots (Q<0)
  85         SQ = math.sqrt(-Q)
  86         arg = R / (SQ**3)
  87         if arg >= 1:
  88             theta = 0
  89         elif arg <= -1:
  90             theta = math.pi
  91         else:
  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)]
  94
  95
  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]
 101     if not ys:
 102         return []
 103     y1 = min(ys)
 104     if a3*y1-2*a1 < 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))))
 107     else:
 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))))
 110
 111
 112 def realpolyroots(*cs):
 113     """returns the roots of a polynom with given coefficients
 114
 115     polynomial with coefficients given in cs:
 116       0 = \sum_i  cs[i] * x^(len(cs)-i-1)
 117     """
 118     if not cs:
 119         return [0]
 120     try:
 121         f = 1.0/cs[0]
 122         cs = [f*c for c in cs[1:]]
 123     except ArithmeticError:
 124         return realpolyroots(*cs[1:])
 125     else:
 126         n = len(cs)
 127         if n == 0:
 128             return []
 129         elif n == 1:
 130             return [-cs[0]]
 131         elif n == 2:
 132             return _realroots_quadratic(*cs)
 133         elif n == 3:
 134             return _realroots_cubic(*cs)
 135         elif n == 4:
 136             return _realroots_quartic(*cs)
 137         else:
 138             raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
 139
 140
 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")
 146 #     if not cs:
 147 #         return [0]
 148 #     try:
 149 #         f = 1.0/cs[0]
 150 #         cs = [f*c for c in cs[1:]]
 151 #     except ArithmeticError:
 152 #         return realpolyroots_eigenvalue(*cs[1:])
 153 #     else:
 154 #         if not cs:
 155 #             return []
 156 #         n = len(cs)
 157 #         a = Numeric.zeros((n, n), Numeric.Float)
 158 #         for i in range(n-1):
 159 #             a[i+1][i] = 1
 160 #         for i in range(n):
 161 #             a[0][i] = -cs[i]
 162 #         rs = []
 163 #         for r in LinearAlgebra.eigenvalues(a):
 164 #             if type(r) == types.ComplexType:
 165 #                 if not r.imag:
 166 #                     rs.append(r.real)
 167 #             else:
 168 #                 rs.append(r)
 169 #         return rs
 170 #