Polynomials And Rational Functions

Introduction

The ratfun module provides classes for defining polynomial and rational function (ratio of two polynomials) objects. These objects can be used in arithmetic expressions and evaluated at a particular point.

Features

Limitations

Installation

The ratfun module is distributed as part of the rpncalc package. See the installation instructions for that package. If you don't want the entire rpncalc package, install the clnum package and copy the ratfun module to anywhere on your Python path.

Examples

These examples assume the following imports.

from ratfun import Polynomial as poly, RationalFunction as rat
from clnum import mpq

If you are using the module from the RPN calculator package, the import is slightly modified.

from rpncalc.ratfun import Polynomial as poly, RationalFunction as rat
from clnum import mpq

Object Creation

>>> print poly(0)
0
>>> print poly(1)
1
>>> print poly(1+1j)
(1+1j)
>>> print poly(1,2,3)
3*x^2+2*x+1
>>> # (coefficient, exponent) pairs form.  Pairs can be in any order.
>>> print poly((mpq(1,3),3), (mpq(1,2),2))
1/3*x^3+1/2*x^2
>>> print poly([(mpq(1,3),3), (mpq(1,2),2)]) # Sequence of pairs is accepted
1/3*x^3+1/2*x^2
>>> print rat(0)
0
>>> print rat(0.5)
0.5
>>> print rat(1+1j)
(1+1j)
>>> print rat(1,2)
1/2
>>> print rat((-1,1),(1,1))
x-1
---
x+1
>>> print rat((1,-2,1),(1,2,1))
x^2-2*x+1
---------
x^2+2*x+1
>>> print rat((1,2,1),(2,2))  # Reduced to lowest terms and normalized.
1/2*x+1/2

Arithmetic

>>> x, x2, x3 = poly([(1,1)]), poly([(1,2)]), poly([(1,3)])
>>> print x, x2, x3
x x^2 x^3
>>> print x, x**2, x**3
x x^2 x^3
>>> print .5*x + 1.3*x2 + (1+1j)*x3 + 1
(1+1j)*x^3+(1.3+0j)*x^2+(0.5+0j)*x+1
>>> print poly(1,0.5,1.3,1+1j)
(1+1j)*x^3+1.3*x^2+0.5*x+1
>>> print .5*x + 1.3*x2 + (1+1j)*x3 + 1 - poly(1,0.5,1.3,1+1j)
0
>>> print .5*x + mpq(1,3)*x2
1/3*x^2+0.5*x
>>> a = 1*x3 + 3*x2 + 3*x + 1
>>> b = poly(1,1)
>>> print a,b
x^3+3*x^2+3*x+1 x+1
>>> print a/b
x^2+2*x+1
>>> print a%b
0
>>> q,r = divmod(a,b)
>>> print q,r
x^2+2*x+1 0
>>> q,r = divmod(poly(-1,3,-3,1), poly(1,2,1))
>>> print q,r
x-5 12*x+4
>>> print q*poly(1,2,1) + r
x^3-3*x^2+3*x-1
>>> z = rat((-1,1),(1,1))
>>> print z
x-1
---
x+1
>>> print z*z
x^2-2*x+1
---------
x^2+2*x+1
>>> print z**3
x^3-3*x^2+3*x-1
---------------
x^3+3*x^2+3*x+1
>>> print z**-3
x^3+3*x^2+3*x+1
---------------
x^3-3*x^2+3*x-1
>>> print z**3 * z**-3
1
>>> p = poly(1,2,1)
>>> print p
x^2+2*x+1
>>> print p+z
x^3+3*x^2+4*x
-------------
     x+1
>>> print z-p
-x^3-3*x^2-2*x-2
----------------
      x+1
>>> print p*z
x^2-1
>>> rat.screenWidth = 0 # Suppress multi-line string format.
>>> print p/z
(x^3+3*x^2+3*x+1) / (x-1)
>>> rat.screenWidth = 80
>>> print z/p
      x-1
---------------
x^3+3*x^2+3*x+1

Methods

Polynomials

>>> p = poly(1,2,1)
>>> print p
x^2+2*x+1
>>> print p.coef(0), p.coef(1), p.coef(2), p.coef(3), p.coef(4)
1 2 1 0 0
>>> print p.deg    # Degree of the polynomial
2
>>> print p.deriv  # Derivative of the polynomial
2*x+2
>>> print p.integ  # Integral of the polynomial
1/3*x^3+x^2+x
>>> print list(p.sample(0,2,5)) # Evaluate at 5 points in the interval [0..2]
[(0.0, 1.0), (0.5, 2.25), (1.0, 4.0), (1.5, 6.25), (2.0, 9.0)]
>>> print p.float()   # Convert the coefficients to float
x^2+2*x+1
>>> print p.complex() # Convert the coefficients to complex
x^2+(2+0j)*x+1
>>> print list(poly(1,0,0,0,2).coefAsPairs())
[(1, 0), (2, 4)]

Rational Functions

>>> z = rat((-1,1),(1,1))
>>> print z
x-1
---
x+1
>>> print z.numer  # Numerator polynomial
x-1
>>> print z.denom  # Denominator polynomial
x+1
>>> print z.deriv  # Derivative of the rational function
    2
---------
x^2+2*x+1
>>> print list(z.sample(0,2,5)) # Evaluate at 5 points in the interval [0..2]
[(0.0, -1.0), (0.5, -0.33333333333333331), (1.0, 0.0), (1.5, 0.20000000000000001), (2.0, 0.33333333333333331)]
>>> print z.float()   # Convert the coefficients to float
x-1
---
x+1
>>> print z.complex() # Convert the coefficients to complex
x-1
---
x+1

Extended Example

For a more complex example, let's compute a rational function approximation to the natural log function. This expansion can be found in a number of math handbooks. Note that the expansion is done by evaluating a polynomial at a rational function.

>>> z = rat((-1,1),(1,1))
>>> print z
x-1
---
x+1
>>> p = poly([(mpq(1,n), n) for n in xrange(1,8,2)])
>>> print p
1/7*x^7+1/5*x^5+1/3*x^3+x
>>> aln = 2*p(z)
>>> print aln
352/105*x^7+112/15*x^6+112/5*x^5-112/5*x^2-112/15*x-352/105
-----------------------------------------------------------
        x^7+7*x^6+21*x^5+35*x^4+35*x^3+21*x^2+7*x+1
>>> print float(aln(1))  # Evaluate a rational function by calling it.
0.0
>>> print float(aln(2))
0.693134757332
>>> from math import log
>>> err = [abs(y - log(x)) for x,y in aln.float().sample(1,2,100)]
>>> print len(err)
100
>>> print err[0], err[50], err[-1]
2.42861286637e-17 1.26519925481e-07 1.24232276572e-05
>>> print max(err)
1.24232276572e-05

Roots of Polynomials

There are a number of discussions of polynomial root finders on the Internet (see Numerical Recipes or Wikipedia). Laguerre's method is used in the ratfun root finder because it handles ill-conditioned polynomials, complex roots, and multiple roots (convergence is only linear in this case). Unlike the version in Numerical Recipes, this implementation uses high precision complex arithmetic to avoid non-converging limit cycles.

Since the Wilkinson polynomial is ill-conditioned, we will use it as an example.

>>> from ratfun import polyFromRoots, sortRoots
>>> p = polyFromRoots(range(1,21))
>>> roots = p.roots()
>>> for r in roots: print r,
1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0

In the example above, all of the roots were found exactly. Now use the root finder in Numeric on the same polynomial.

>>> import MLab
>>> roots = sortRoots(MLab.roots(p.float()._coef[::-1]))
>>> for r in roots: print r
1.0
2.00000000004
2.99999999789
4.00000004383
4.99999953503
6.00000250568
6.99999869761
7.9999128424
9.0007574115
9.99630988166
11.0127366347
11.9694308696
13.0605362206
13.9178177587
15.0885739494
15.9278470944
17.0383522601
17.9846032298
19.0034961547
19.9996249124
>>> print max([abs(((i+1)-r)/(i+1)) for i,r in enumerate(roots)])
0.0059049299583

Investigate the source of the ill-conditioned behavior by evaluating the polynomial near the largest root. This demonstrates that the polynomial oscillates with very large excursions between the roots.

>>> print len(str(20**20)) # Number of digits needed
27
>>> prec = 37 # Add another 10 digits
>>> dx = mpf('1e-9', prec)
>>> print p(20+dx), p(20-dx)
1.2164510084039714746e8 -1.21645099977266853944e8
>>> dx = mpf('1e-18', prec)
>>> print p(20+dx), p(20-dx)
0.12164510045845919576 -0.121645100372629827505

Now change one of the coefficients by a small amount and observe what happens to the roots.

>>> p1 = p + poly([(mpq(-1,2**23),19)])
>>> print p.coef(19), float(p1.coef(19))
-210 -210.000000119
>>> roots = p1.roots(eps=1e-20)
>>> for r in roots: print r
1.0
2.0000000000000000098
2.999999999999805233
4.000000000261023189
4.9999999275515379094
6.0000069439522957073
6.999697233936013949
8.007267603450376855
8.917250248517070495
(10.095266145129963366-0.643500903863603576j)
(10.095266145129963366+0.643500903863603576j)
(11.79363388107943398-1.6523297281609322825j)
(11.79363388107943398+1.6523297281609322825j)
(13.992358137235671092-2.518830069630272286j)
(13.992358137235671092+2.518830069630272286j)
(16.730737466090704483-2.81262489427003927j)
(16.730737466090704483+2.81262489427003927j)
(19.502439400493681724-1.9403303466644795428j)
(19.502439400493681724+1.9403303466644795428j)
(20.846908101482256915+5.7070363363468749757e-199j)
>>> print max([float(abs(p1(r))) for r in roots])
0.131058942236

Multiple Roots

A polynomial can have a root repeated multiple times. Let p=(x - r)m, then p is a polynomial with a single root r of multiplicity m. This type of polynomial causes problems for iterative root-finding methods like Laguerre's method because the graph of the polynomial is flat in the region surrounding the root. The higher the multiplicity, the flatter the graph. As a consequence of this, Laguerre's method converges only linearly when approaching a multiple root. So what can be done to make the convergence faster?

If you take the derivative of a polynomial with a root of multiplicity greater than one, you will find that the multiplicity is reduced by one. Computing the greatest common divisor of the derivative polynomial with the original will identify all the extra common factors which can then be removed by dividing into the original polynomial. The resulting polynomial will contain only unique roots. This polynomial can then be solved using a standard root finder. The following example shows this reduction.

>>> p = poly(1,1)**5
>>> print p
x^5+5*x^4+10*x^3+10*x^2+5*x+1
>>> q = gcd(p, p.deriv)
>>> print p/q
1/5*x+1/5
>>> print (p/q).roots()
[mpf('-1.0',17)]

The method uniqueRoots performs the above operations and returns a list of the unique roots. It is often possible to find the roots exactly with this technique.

>>> p = poly(-mpq(1,3),1)**5 * poly(-mpq(1,7),1)**7 *  poly(-mpq(1,11),1)**10
>>> print p.deg
22
>>> for r in p.uniqueRoots(): print mpq(r),
1/11 1/7 1/3

If the coefficients of the polynomial are not exact, uniqueRoots may not be able to determine the multiplicity. In this case, it drops back to using the roots method on the original polynomial. Consequently, you loose nothing but the time it takes to compute a derivative and a gcd.

The multiplicity of each root can be determined by repeatedly dividing by the monomial factor until there is a remainder. However, doing this reliably with approximate roots is a difficult problem so there is no method provided to determine the multiplicities.

SourceForge.net Logo