<p dir="ltr">Actually, GCC implements 128-bit floats in software and provides them as __float128; there are also quad-precision versions of the usual functions. The Intel compiler provides this as well, I think, but I don't think Microsoft compilers do. A portable quad-precision library might be less painful.</p>
<p dir="ltr">The cleanest way to add extended precision to numpy is by adding a C-implemented dtype. This can be done in an extension module; see the quaternion and half-precision modules online.</p>
<p dir="ltr">Anne</p>
<br><div class="gmail_quote"><div dir="ltr">On Fri, Dec 11, 2015, 16:46¬†Charles R Harris <<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Dec 11, 2015 at 6:25 AM, Thomas Baruchel <span dir="ltr"><<a href="mailto:baruchel@gmx.com" target="_blank">baruchel@gmx.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">From time to time it is asked on forums how to extend precision of computation on Numpy array. The most common answer<br>
given to this question is: use the dtype=object with some arbitrary precision module like mpmath or gmpy.<br>
See <a href="http://stackoverflow.com/questions/6876377/numpy-arbitrary-precision-linear-algebra" rel="noreferrer" target="_blank">http://stackoverflow.com/questions/6876377/numpy-arbitrary-precision-linear-algebra</a> or <a href="http://stackoverflow.com/questions/21165745/precision-loss-numpy-mpmath" rel="noreferrer" target="_blank">http://stackoverflow.com/questions/21165745/precision-loss-numpy-mpmath</a> or <a href="http://stackoverflow.com/questions/15307589/numpy-array-with-mpz-mpfr-values" rel="noreferrer" target="_blank">http://stackoverflow.com/questions/15307589/numpy-array-with-mpz-mpfr-values</a><br>
<br>
While this is obviously the most relevant answer for many users because it will allow them to use Numpy arrays exactly<br>
as they would have used them with native types, the wrong thing is that from some point of view "true" vectorization<br>
will be lost.<br>
<br>
With years I got very familiar with the extended double-double type which has (for usual architectures) about 32 accurate<br>
digits with faster arithmetic than "arbitrary precision types". I even used it for research purpose in number theory and<br>
I got convinced that it is a very wonderful type as long as such precision is suitable.<br>
<br>
I often implemented it partially under Numpy, most of the time by trying to vectorize at a low-level the libqd library.<br>
<br>
But I recently thought that a very nice and portable way of implementing it under Numpy would be to use the existing layer<br>
of vectorization on floats for computing the arithmetic operations by "columns containing half of the numbers" rather than<br>
by "full numbers". As a proof of concept I wrote the following file: <a href="https://gist.github.com/baruchel/c86ed748939534d8910d" rel="noreferrer" target="_blank">https://gist.github.com/baruchel/c86ed748939534d8910d</a><br>
<br>
I converted and vectorized the Algol 60 codes from <a href="http://szmoore.net/ipdf/documents/references/dekker1971afloating.pdf" rel="noreferrer" target="_blank">http://szmoore.net/ipdf/documents/references/dekker1971afloating.pdf</a><br>
(Dekker, 1971).<br>
<br>
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<br>
50 times faster than MPmath. It should be even faster for some other operations since I had to create another np.ones<br>
array for testing this type because inversion isn't implemented here (which could of course be done). You can run this file by yourself<br>
(maybe you will have to discard mpmath or gmpy if you don't have it).<br>
<br>
I would like to discuss about the way to make available something related to that.<br>
<br>
a) Would it be relevant to include that in Numpy ? (I would think to some "contribution"-tool rather than including it in<br>
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<br>
to perform several arithmetic operations by knowing that they can't use cos/sin/etc. on this type; in other words, I am not<br>
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<br>
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<br>
be very happy to code it in some cleaner way)?<br>
<br>
b) Do you think such attempt should remain something external to Numpy itself and be released on my Github account without being<br>
integrated to Numpy?<br></blockquote><div><br></div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>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 :(.<br><br></div><div>Chuck <br></div></div></div></div>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org" target="_blank">NumPy-Discussion@scipy.org</a><br>
<a href="https://mail.scipy.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div>