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.
Create polynomials and rational functions with an arbitrary number of coefficients. The coefficients can be any numeric type (conversion to complex is the test of a numeric type).
Support the usual arithmetic operations (+, -, *, /, **). Polynomials also support %, and divmod.
Combine polynomials, rational functions, and ordinary numbers together in arithmetic expressions.
Automatically handles rational and extended precision floating point coefficients. This includes converting integers to rationals when a division is performed.
Evaluate a polynomial or rational function by calling it.
Find the roots of a polynomial to any desired precision. Using the default tolerance, the method converges exactly even for the ill-conditioned Wilkinson polynomial of degree 100.
Only supports polynomials and rational functions in a single variable.
Mixing Python float and complex numbers with polynomials or rational functions that have mpq or cmpq coefficients cause exceptions. There are methods to convert all of the coefficients to an appropriate type to overcome this problem.
Operating on high degree rational functions can get very slow due to the GCD calculations which are not optimized.
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
>>> 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
>>> 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
>>> 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)]
>>> 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
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
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
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.