[Numpy-discussion] Re: status of long double support and what to do about it

[Scott Ransom]

On 11/17/22 7:13 PM, Charles R Harris wrote:



On Thu, Nov 17, 2022 at 3:15 PM Ralf Gommers > wrote:


Hi all,

We have to do something about long double support. This is something I 
wanted to propose a long
time ago already, and moving build systems has resurfaced the pain yet 
again.

This is not a full proposal yet, but the start of a discussion and gradual 
plan of attack.

The problem
---
The main problem is that long double support is *extremely* painful to 
maintain, probably far
more than justified. I could write a very long story about that, but 
instead I'll just
illustrate with some of the key points:

(1) `long double` is the main reason why we're having such a hard time with 
building wheels on
Windows, for SciPy in particular. This is because MSVC makes long double 
64-bit, and Mingw-w64
defaults to 80-bit. So we have to deal with Mingw-w64 toolchains, proposed 
compiler patches,
etc. This alone has been a massive time sink. A couple of threads:
https://github.com/numpy/numpy/issues/20348 


https://discuss.scientific-python.org/t/releasing-or-not-32-bit-windows-wheels/282


The first issue linked above is one of the key ones, with a lot of detail 
about the fundamental
problems with `long double`. The Scientific Python thread focused more on 
Fortran, however that
Fortran + Windows problem is at least partly the fault of `long double`. 
And Fortran may be
rejuvenated with LFortran and fortran-lang.org ; 
`long double` is a
dead end.

(2) `long double` is not a well-defined format. We support 9 specific 
binary representations,
and have a ton of code floating around to check for those, and manually 
fiddle with individual
bits in long double numbers. Part of that is the immediate pain point for 
me right now: in the
configure stage of the build we consume object files produced by the 
compiler and parse them,
matching some known bit patterns. This check is so weird that it's the only 
one that I cannot
implement in Meson (short of porting the hundreds of lines of Python code 
for it to C), see
https://github.com/mesonbuild/meson/issues/11068
 for details. To get an 
idea of the complexity
here:

https://github.com/numpy/numpy/blob/9e144f7c1598221510d49d8c6b79c66dc000edf6/numpy/core/setup_common.py#L264-L434
 


https://github.com/numpy/numpy/blob/9e144f7c1598221510d49d8c6b79c66dc000edf6/numpy/core/src/npymath/npy_math_private.h#L179-L484
 


https://github.com/numpy/numpy/blob/main/numpy/core/src/npymath/npy_math_complex.c.src#L598-L619



https://github.com/numpy/numpy/blob/main/numpy/core/src/multiarray/dragon4.c#L2480-L3052


Typically `long double` has multiple branches, and requires more code than 
float/double.

(3) We spend a lot of time dealing with issues and PR to keep `long double` 
working, as well as
dealing with hard-to-diagnose build issues. Which sometimes even stops 
people from
building/contributing, especially on Windows. Some recent examples:
https://github.com/numpy/numpy/pull/20360 

https://github.com/numpy/numpy/pull/18536 

https://github.com/numpy/numpy/pull/21813 

https://github.com/numpy/numpy/pull/22405 

https://github.com/numpy/numpy/pull/19950 

https://github.com/numpy/numpy/pull/18330/commits/aa9fd3c7cb

https://github.com/scipy/scipy/issues/16769 

https://github.com/numpy/numpy/issues/14574 


(4) `long double` isn't all that useful. On both Windows and macOS `long 
double` is 64-bit,
which means it is just a poor alias to `double`. So it does literally 
nothing for the majority
of our users, except confuse them and take up extra memory. On Linux, `long 
double` is 80-bit
precision, which means it doesn't do all that much there either, just a 
modest bump in precision.

Let me also note that it's not just the user-v

[Numpy-discussion] Re: status of long double support and what to do about it

[Scott Ransom]

On 11/17/22 8:53 PM, Charles R Harris wrote:



On Thu, Nov 17, 2022 at 6:30 PM Scott Ransom mailto:sran...@nrao.edu>> wrote:



On 11/17/22 7:13 PM, Charles R Harris wrote:
 >
 >
 > On Thu, Nov 17, 2022 at 3:15 PM Ralf Gommers mailto:ralf.gomm...@gmail.com>
 > <mailto:ralf.gomm...@gmail.com <mailto:ralf.gomm...@gmail.com>>> wrote:
 >
 >     Hi all,
 >
 >     We have to do something about long double support. This is something 
I wanted to propose
a long
 >     time ago already, and moving build systems has resurfaced the pain 
yet again.
 >
 >     This is not a full proposal yet, but the start of a discussion and 
gradual plan of attack.

 > I would agree that extended precision is pretty useless, IIRC, it was 
mostly intended as an
accurate
 > way to produce double precision results. That idea was eventually 
dropped as not very useful.
I'd
 > happily do away with subnormal doubles as well, they were another not 
very useful idea. And
strictly
 > speaking, we should not support IBM double-double either, it is not in 
the IEEE standard.
 >
 > That said, I would like to have a quad precision type. That precision is 
useful for some
things, and
 > I have a dream that someday it can be used for a time type. 
Unfortunately, last time I looked
 > around, none of the available implementations had a NumPy compatible 
license.
 >
 > The tricky thing here is to not break downstream projects, but that may 
be unavoidable. I
suspect
 > the fallout will not be that bad.
 >
 > Chuck

A quick response from one of the leaders of a team that requires 80bit 
extended precision for
astronomical work...

"extended precision is pretty useless" unless you need it. And the 
high-precision pulsar timing
community needs it. Standard double precision (64-bit) values do not 
contain enough precision
for us
to pass relative astronomical times via a single float without extended 
precision (the precision
ends up being at the ~1 microsec level over decades of time differences, 
and we need it at the
~1-10ns level) nor can we store the measured spin frequencies (or do 
calculations on them) of our
millisecond pulsars with enough precision. Those spin frequencies can have 
16-17 digits of base-10
precision (i.e. we measure them to that precision). This is why we use 
80-bit floats (usually via
Linux, but also on non X1 Mac hardware if you use the correct compilers) 
extensively.

Numpy is a key component of the PINT software to do high-precision pulsar 
timing, and we use it
partly *because* it has long double support (with 80-bit extended 
precision):
https://github.com/nanograv/PINT <https://github.com/nanograv/PINT>
And see the published paper here, particularly Sec 3.3.1 and footnote #42:
https://ui.adsabs.harvard.edu/abs/2021ApJ...911...45L/abstract
<https://ui.adsabs.harvard.edu/abs/2021ApJ...911...45L/abstract>

Going to software quad precision would certainly work, but it would 
definitely make things much
slower for our matrix and vector math.

We would definitely love to see a solution for this that allows us to get 
the extra precision we
need on other platforms besides Intel/AMD64+Linux (primarily), but giving 
up extended precision on
those platforms would *definitely* hurt. I can tell you that the pulsar 
community would definitely
be against option "B". And I suspect that there are other users out there 
as well.

Scott
NANOGrav Chair
www.nanograv.org <http://www.nanograv.org>



Pulsar timing is one reason I wanted a quad precision time type. I thought Astropy was using a self 
implemented double-double type to work around that?


That is correct. For non-compute-intensive time calculations, Astropy as a Time object that 
internally uses two 64-bit floats. We use it, and it works great for high precision timekeeping over 
astronomical times.


*However*, it ain't fast. So you can't do fast matrix/vector math on time differences where your 
precision exceeds a single 64-bit float. That's exactly where we are with extended precision for our 
pulsar timing work.


Scott

--
Scott M. RansomAddress:  NRAO
Phone:  (434) 296-0320   520 Edgemont Rd.
email:  sran...@nrao.edu Charlottesville, VA 22903 USA
GPG Fingerprint: A40A 94F2 3F48 4136 3AC4  9598 92D5 25CB 22A6 7B65
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: status of long double support and what to do about it

[Scott Ransom]

Hi Mike,

Good to hear from you!

On 11/17/22 10:36 PM, Mike McCarty wrote:
Hi Scott - Thanks for providing input! Do you have a minimal example that shows the type of 
calculations you would like to be faster with extended precision? Just so we are clear on the goal 
for this use case.


The majority of what we do is have a vector of extended precision times that we do various kinds of 
algebra on, using many NumPy routines and vectorized custom functions.


We don't usually have to do too many LAPACK-style heavy-duty matrix operations 
in extended precision.


Would you or some of your colleagues be open to helping maintain a library 
that adds the 80-bit
extended precision dtype into NumPy? This would be a variation of Ralf's "option 
A."

The NumPy community has resources and assistance for onboarding new contributors, and I'm happy to 
provide additional assistance.


There are possibly a couple of us who might be willing to have a look at that. It is important 
enough for what we do that we really need a solution.


Scott

--
Scott M. RansomAddress:  NRAO
Phone:  (434) 296-0320   520 Edgemont Rd.
email:  sran...@nrao.edu Charlottesville, VA 22903 USA
GPG Fingerprint: A40A 94F2 3F48 4136 3AC4  9598 92D5 25CB 22A6 7B65
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: status of long double support and what to do about it

[Scott Ransom]

On 12/14/22 3:01 PM, Carl Kleffner wrote:

Hi Scott,

may I ask you which kind of vector / matrix operation in extended precision (np.longdouble) is 
supported in 'pint' ? It can't be backed by the underlying blas library as extended precision 
functions are not supported there.


Carl


Hi Carl,

The vast majority of them are simply a vector of np.longdoubles (which are usually times or time 
differences) which are operated on by regular Numpy math functions to compute new quantities.  For 
example, polynomial computations base on time, many trigonometric operations on the times for orbit 
calculations, things like that.


We don't really do big longdouble matrix operations.

Scott

--
Scott M. RansomAddress:  NRAO
Phone:  (434) 296-0320   520 Edgemont Rd.
email:  sran...@nrao.edu Charlottesville, VA 22903 USA
GPG Fingerprint: A40A 94F2 3F48 4136 3AC4  9598 92D5 25CB 22A6 7B65
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com