[Numpy-discussion] Idea: fractional polynomial class

[Pauli Virtanen]

Hi,

Fri, 24 Apr 2009 12:52:33 -0400, Michael S. Gilbert kirjoitti:
> I've been working with numpy's poly1d class recently, and it would be
> very useful to me if the class worked with fractions instead of floats
> (since I'm encountering quantities that often mostly cancel out, which
> lead to instabilities in my algorithms; hence it would be useful to use
> fractions that don't lose precision). Since python 2.6 introduced a
> "fractions" class, it should be fairly straightforward to
> implement/reimplement the polynomial class with support for fractions.
> 
> I'm curious as to whether there would be any interest from the project
> on this and whether you all would be open to patches/code contributions.

Definitely, provided fractions support the usual Python operators.
(I object to having a specialized 'fractional' polynomial class, though.)

The changes required are probably very small: only `polyint` needs some 
fixes.

So, thanks to duck typing, it mostly works out of the box already. For 
example, I can use mpmath's arbitrary precision numbers in poly1d:

>>> import numpy as np
>>> from mpmath import mp, mpf
>>> p = np.poly1d([mpf('1.0'), mpf('2.0'), mpf('3.0')])
>>> p(1/mpf('3'))
mpf('3.777777777777777777777777777777777777777773')
>>> p * mpf('1.9999999999999999999999999')
poly1d([1.9999999999999999999999999, 3.9999999999999999999999998,
       5.9999999999999999999999997], dtype=object)
>>> p*p
poly1d([1.0, 4.0, 10.0, 12.0, 9.0], dtype=object)
>>> np.polyder(p*p, 3)
poly1d([24.0, 24.0], dtype=object)

Polyder uses only exact integers, so it does not lose precision.

The only thing that does not work is integration, since it seems to 
force-cast to float:

>>> np.polyint(p)
poly1d([ 0.33333333,  1.        ,  3.        ,  0.        ])

This is a bug in itself, since poly1d could well have complex 
coefficients.

-- 
Pauli Virtanen