[Numpy-discussion] Idea: fractional polynomial class
Pauli Virtanen
pav at iki.fi
Fri Apr 24 13:09:06 EDT 2009
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
More information about the NumPy-Discussion
mailing list