pyx/helper.py

   1 #!/usr/bin/env python
   2 # -*- coding: ISO-8859-1 -*-
   3 #
   4 #
   5 # Copyright (C) 2002-2004 Jörg Lehmann <joergl@users.sourceforge.net>
   6 # Copyright (C) 2002-2005 André Wobst <wobsta@users.sourceforge.net>
   7 # Copyright (C) 2005 Michael Schindler <m-schindler@users.sourceforge.net>
   8 #
   9 # This file is part of PyX (http://pyx.sourceforge.net/).
  10 #
  11 # PyX is free software; you can redistribute it and/or modify
  12 # it under the terms of the GNU General Public License as published by
  13 # the Free Software Foundation; either version 2 of the License, or
  14 # (at your option) any later version.
  15 #
  16 # PyX is distributed in the hope that it will be useful,
  17 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 # GNU General Public License for more details.
  20 #
  21 # You should have received a copy of the GNU General Public License
  22 # along with PyX; if not, write to the Free Software
  23 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  24
  25
  26 def sign(x):
  27     return (x >= 0) and 1 or -1
  28
  29 # XXX fallback for Numeric (eigenvalue computation) to be implemented along
  30 # know algorithms (like from numerical recipes)
  31
  32 import Numeric, LinearAlgebra
  33
  34 def realpolyroots(coeffs, epsilon=1e-5, polish=0):
  35
  36     """returns the roots of a polynom with given coefficients
  37
  38     This helper routine uses the package Numeric to find the roots
  39     of the polynomial with coefficients given in coeffs:
  40       0 = \sum_{i=0}^N x^i coeffs[i]
  41     The solution is found via an equivalent eigenvalue problem
  42     """
  43
  44     try:
  45         1.0 / coeffs[-1]
  46     except:
  47         roots = realpolyroots(coeffs[:-1], epsilon=epsilon)
  48     else:
  49         # divide the coefficients by the maximal order
  50         # this is to be done in the matrix, anyhow and
  51         # makes checking more stable
  52         coeffs = [coeff/coeffs[-1] for coeff in coeffs]
  53
  54         N = len(coeffs) - 1
  55         # build the Matrix of the polynomial problem
  56         mat = Numeric.zeros((N, N), Numeric.Float)
  57         for i in range(N-1):
  58             mat[i+1][i] = 1
  59         for i in range(N):
  60             mat[0][i] = -coeffs[N-1-i]
  61         # find the eigenvalues of the matrix (== the roots of the polynomial)
  62         roots = [complex(root) for root in LinearAlgebra.eigenvalues(mat)]
  63         # take only the real roots
  64         roots = [root.real for root in roots if -epsilon < root.imag < epsilon]
  65
  66         # polish the roots with Newton-Raphson
  67         if polish:
  68             def polish(root, epsilon):
  69                 polynom = 2*epsilon
  70                 while abs(polynom) > epsilon:
  71                     polynom = coeffs[N]*root + coeffs[N-1]
  72                     poprime = coeffs[N]*N
  73                     for i in range(N-2,-1,-1):
  74                         polynom = polynom*root + coeffs[i]
  75                         poprime = poprime*root + coeffs[i+1]*(i+1)
  76                     root -= polynom / poprime
  77                 return root
  78
  79             roots = [polish(root, epsilon) for root in roots]
  80
  81         # # check if the roots are really roots!
  82         # for root in roots:
  83         #     polynom = coeffs[N]*root + coeffs[N-1]
  84         #     for i in range(N-2,-1,-1):
  85         #         polynom = polynom*root + coeffs[i]
  86         #     if abs(polynom) > epsilon:
  87         #         raise Exception("value %f instead of 0" % polynom)
  88
  89     return roots
  90