On Mon, May 6, 2013 at 7:57 AM, Nathaniel Smith <njs@pobox.com> wrote:
On Mon, May 6, 2013 at 8:52 AM, Funky Dev <funkydevnull@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.
The rational test class recently merged into numpy is more complete and fixed up than the version in numpy-dtypes, although I can't help thinking it could be streamlined even more.
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.
Chuck