[Numpy-discussion] mpfr integration into numpy?

[Nathaniel Smith]

On Mon, May 6, 2013 at 8:52 AM, Funky Dev <funkydevnull at gmail.com> wrote:
> Hi,
>
> I've got a project in which it turns out we need much higher precision
> than even __float128 and playing around with a few alternatives mpfr
> seems to be the highest performing possibility.  So I've starting
> writing a cythonized class mpfr_array which provides array-like
> functionality but with mpfr_t as a "primitive" type.
>
> While this is doable it feels very much like re-inventing the wheel as
> I'm basically rewriting some part of numpy's functionality.  So I was
> wondering, as I think this is a very standard use case, is there any
> interest in adding mpfr_t as a "primitive" dtype to numpy?  I know one
> can define new dtype's but I don't want it to by a python object since
> then there will be a large python overhead for all operations (i.e. dot
> products will happen in python, not C).  mpfr_t is a very natural dtype
> to add as its fast, C-based and supports general precision.

You actually can define new dtypes implemented in C code -- in fact,
right now, that's the only way to do it. There are some examples of
how to do this from third-party code here:
  https://github.com/numpy/numpy-dtypes
Since mpfr is GPLed, it couldn't be included in core numpy in any
case, but I don't see any reason you couldn't implement it this way
and distribute it as a third-party extension.

There would be some subtleties to work out regarding memory management
-- I guess mpfr_t objects have a fixed width (which is a requirement
for dtypes), but internally malloc some buffers, so you have to be
careful with copying them etc.? Should be doable with some care.

> On a related note, if this was done, would it automatically work with
> functionality such as numpy.linalg.inv(), etc...?  In principle such
> functions could have been written with macros to be more type-flexible
> (i.e. an ADD(a,b,c) mapping to e.g. a=b+c for floats or to mpfr_add(a,
> b, c) for mpfr_t) but I suspect this is not the case.

No, it wouldn't work automatically. But in numpy 1.8 you can register
new implementations of 'inv' and friends when your library is loaded,
which np.linalg.inv() will then automatically use if it sees your new
dtype.

-n