[Numpy-discussion] Fast vectorized arithmetic with ~32 significant digits under Numpy

Charles R Harris charlesr.harris at gmail.com
Fri Dec 11 12:51:06 EST 2015


On Fri, Dec 11, 2015 at 10:45 AM, Nathaniel Smith <njs at pobox.com> wrote:

> On Dec 11, 2015 7:46 AM, "Charles R Harris" <charlesr.harris at gmail.com>
> wrote:
> >
> >
> >
> > On Fri, Dec 11, 2015 at 6:25 AM, Thomas Baruchel <baruchel at gmx.com>
> wrote:
> >>
> >> From time to time it is asked on forums how to extend precision of
> computation on Numpy array. The most common answer
> >> given to this question is: use the dtype=object with some arbitrary
> precision module like mpmath or gmpy.
> >> See
> http://stackoverflow.com/questions/6876377/numpy-arbitrary-precision-linear-algebra
> or http://stackoverflow.com/questions/21165745/precision-loss-numpy-mpmath
> or
> http://stackoverflow.com/questions/15307589/numpy-array-with-mpz-mpfr-values
> >>
> >> While this is obviously the most relevant answer for many users because
> it will allow them to use Numpy arrays exactly
> >> as they would have used them with native types, the wrong thing is that
> from some point of view "true" vectorization
> >> will be lost.
> >>
> >> With years I got very familiar with the extended double-double type
> which has (for usual architectures) about 32 accurate
> >> digits with faster arithmetic than "arbitrary precision types". I even
> used it for research purpose in number theory and
> >> I got convinced that it is a very wonderful type as long as such
> precision is suitable.
> >>
> >> I often implemented it partially under Numpy, most of the time by
> trying to vectorize at a low-level the libqd library.
> >>
> >> But I recently thought that a very nice and portable way of
> implementing it under Numpy would be to use the existing layer
> >> of vectorization on floats for computing the arithmetic operations by
> "columns containing half of the numbers" rather than
> >> by "full numbers". As a proof of concept I wrote the following file:
> https://gist.github.com/baruchel/c86ed748939534d8910d
> >>
> >> I converted and vectorized the Algol 60 codes from
> http://szmoore.net/ipdf/documents/references/dekker1971afloating.pdf
> >> (Dekker, 1971).
> >>
> >> A test is provided at the end; for inverting 100,000 numbers, my type
> is about 3 or 4 times faster than GMPY and almost
> >> 50 times faster than MPmath. It should be even faster for some other
> operations since I had to create another np.ones
> >> array for testing this type because inversion isn't implemented here
> (which could of course be done). You can run this file by yourself
> >> (maybe you will have to discard mpmath or gmpy if you don't have it).
> >>
> >> I would like to discuss about the way to make available something
> related to that.
> >>
> >> a) Would it be relevant to include that in Numpy ? (I would think to
> some "contribution"-tool rather than including it in
> >> the core of Numpy because it would be painful to code all ufuncs; on
> the other hand I am pretty sure that many would be happy
> >> to perform several arithmetic operations by knowing that they can't use
> cos/sin/etc. on this type; in other words, I am not
> >> sure it would be a good idea to embed it as an every-day type but I
> think it would be nice to have it quickly available
> >> in some way). If you agree with that, in which way should I code it
> (the current link only is a "proof of concept"; I would
> >> be very happy to code it in some cleaner way)?
> >>
> >> b) Do you think such attempt should remain something external to Numpy
> itself and be released on my Github account without being
> >> integrated to Numpy?
> >
> >
> > I think astropy does something similar for time and dates. There has
> also been some talk of adding a user type for ieee 128 bit doubles. I've
> looked once for relevant code for the latter and, IIRC, the available
> packages were GPL :(.
>
> You're probably thinking of the __float128 support in gcc, which relies on
> a LGPL (not GPL) runtime support library. (LGPL = any patches to the
> support library itself need to remain open source, but no restrictions are
> imposed on code that merely uses it.)
>
> Still, probably something that should be done outside of numpy itself for
> now.
>

No, there are several other software packages out there. I know of the gcc
version, but was looking for something more portable.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20151211/9a97b54d/attachment.html>


More information about the NumPy-Discussion mailing list