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-2004 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
  25 # XXX fallback for Numeric (eigenvalue computation) to be implemented along
  26 # know algorithms (like from numerical recipes)
  27
  28 import Numeric, LinearAlgebra
  29
  30 def realpolyroots(coeffs, epsilon=1e-5, polish=0):
  31
  32     """returns the roots of a polynom with given coefficients
  33
  34     This helper routine uses the package Numeric to find the roots
  35     of the polynomial with coefficients given in coeffs:
  36       0 = \sum_{i=0}^N x^i coeffs[i]
  37     The solution is found via an equivalent eigenvalue problem
  38     """
  39
  40     try:
  41         1.0 / coeffs[-1]
  42     except:
  43         roots = realpolyroots(coeffs[:-1], epsilon=epsilon)
  44     else:
  45         # divide the coefficients by the maximal order
  46         # this is to be done in the matrix, anyhow and
  47         # makes checking more stable
  48         coeffs = [coeff/coeffs[-1] for coeff in coeffs]
  49
  50         N = len(coeffs) - 1
  51         # build the Matrix of the polynomial problem
  52         mat = Numeric.zeros((N, N), Numeric.Float)
  53         for i in range(N-1):
  54             mat[i+1][i] = 1
  55         for i in range(N):
  56             mat[0][i] = -coeffs[N-1-i]
  57         # find the eigenvalues of the matrix (== the roots of the polynomial)
  58         roots = [complex(root) for root in LinearAlgebra.eigenvalues(mat)]
  59         # take only the real roots
  60         roots = [root.real for root in roots if -epsilon < root.imag < epsilon]
  61
  62         # polish the roots with Newton-Raphson
  63         if polish:
  64             def polish(root, epsilon):
  65                 polynom = 2*epsilon
  66                 while abs(polynom) > epsilon:
  67                     polynom = coeffs[N]*root + coeffs[N-1]
  68                     poprime = coeffs[N]*N
  69                     for i in range(N-2,-1,-1):
  70                         polynom = polynom*root + coeffs[i]
  71                         poprime = poprime*root + coeffs[i+1]*(i+1)
  72                     root -= polynom / poprime
  73                 return root
  74
  75             roots = [polish(root, epsilon) for root in roots]
  76
  77         # # check if the roots are really roots!
  78         # for root in roots:
  79         #     polynom = coeffs[N]*root + coeffs[N-1]
  80         #     for i in range(N-2,-1,-1):
  81         #         polynom = polynom*root + coeffs[i]
  82         #     if abs(polynom) > epsilon:
  83         #         raise Exception("value %f instead of 0" % polynom)
  84
  85     return roots
  86