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         theta = math.acos(R/(SQ**3))
  87         return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)]
  88
  89
  90 def _realroots_quartic(a3, a2, a1, a0):
  91     """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0"""
  92     # see http://mathworld.wolfram.com/QuarticEquation.html for details
  93     ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3)
  94     ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and  y*y-4*a0 >= 0]
  95     if not ys:
  96         return []
  97     y1 = min(ys)
  98     if a3*y1-2*a1 < 0:
  99         return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) +
 100                 _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))))
 101     else:
 102         return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) +
 103                 _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))))
 104
 105
 106 def realpolyroots(*cs):
 107     """returns the roots of a polynom with given coefficients
 108
 109     polynomial with coefficients given in cs:
 110       0 = \sum_i  cs[i] * x^(len(cs)-i-1)
 111     """
 112     if not cs:
 113         return [0]
 114     try:
 115         f = 1.0/cs[0]
 116         cs = [f*c for c in cs[1:]]
 117     except ArithmeticError:
 118         return realpolyroots(*cs[1:])
 119     else:
 120         n = len(cs)
 121         if n == 0:
 122             return []
 123         elif n == 1:
 124             return [-cs[0]]
 125         elif n == 2:
 126             return _realroots_quadratic(*cs)
 127         elif n == 3:
 128             return _realroots_cubic(*cs)
 129         elif n == 4:
 130             return _realroots_quartic(*cs)
 131         else:
 132             raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
 133
 134
 135 # def realpolyroots_eigenvalue(*cs):
 136 #     # as realpolyroots but using an equivalent eigenvalue problem
 137 #     # (this code is currently used for functional tests only)
 138 #     if not _has_numeric:
 139 #         raise RuntimeError("realpolyroots_eigenvalue depends on Numeric")
 140 #     if not cs:
 141 #         return [0]
 142 #     try:
 143 #         f = 1.0/cs[0]
 144 #         cs = [f*c for c in cs[1:]]
 145 #     except ArithmeticError:
 146 #         return realpolyroots_eigenvalue(*cs[1:])
 147 #     else:
 148 #         if not cs:
 149 #             return []
 150 #         n = len(cs)
 151 #         a = Numeric.zeros((n, n), Numeric.Float)
 152 #         for i in range(n-1):
 153 #             a[i+1][i] = 1
 154 #         for i in range(n):
 155 #             a[0][i] = -cs[i]
 156 #         rs = []
 157 #         for r in LinearAlgebra.eigenvalues(a):
 158 #             if type(r) == types.ComplexType:
 159 #                 if not r.imag:
 160 #                     rs.append(r.real)
 161 #             else:
 162 #                 rs.append(r)
 163 #         return rs
 164 #