On Fri, Dec 11, 2015 at 10:45 AM, Nathaniel Smith <n...@pobox.com> wrote:
> On Dec 11, 2015 7:46 AM, "Charles R Harris" <charlesr.har...@gmail.com> > wrote: > > > > > > > > On Fri, Dec 11, 2015 at 6:25 AM, Thomas Baruchel <baruc...@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
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
