[Numpy-discussion] NEP 50: deprecating np.find_common_type

[Inessa Pawson]

As part of the work outlined in NEP 50, the function np.find_common_type
will be deprecated soon after the release of NumPy 1.24. To provide
feedback on this decision, please respond in the Conversation section of
the pull request:  https://github.com/numpy/numpy/pull/22539.

-- 
Cheers,
Inessa

Inessa Pawson
Contributor Experience Lead | NumPy
https://numpy.org/
GitHub: inessapawson
___
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] status of long double support and what to do about it

[Ralf Gommers]

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-visible dtypes that we have to
consider; long double types are also baked into the libnpymath static
library that we ship with NumPy. That's a thing we have to do something
about anyway (shipping static libraries is not the best idea, see
https://github.com/numpy/numpy/issues/20880). We just have to make sure to
not forget about it when thinking about solutions here.


Potential solutions
---

(A) The ideal solution would be to have a proper, portable quad-precision
(128 bits of precision) dtype. It's now possible to write that externally,
after all the work that Sebastian and others have put into the dtype
infrastructure. The dtype itself already exists (
https://github.com/peytondmurray/quaddtype, maybe there are more
implementations floating around). It just need the people who have an
actual need for it to drive that. It's still a significant amount of work,
so I'll not go into this one more right now.

(B) Full deprecation and removal of all `long double` support from NumPy
(and SciPy), irrespective of whether the quad dtype comes to life.

Given the current state, I'm personally convinced that that is easily
justified. However, I know some folks are going to be hesitant, given that
we don't know how many remaining users we have or what use cases they have.
So let's see if we can find more granular solutions (note: these are ideas,
not all fully researched solutions that 

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

[Ilhan Polat]

Thanks Ralf, this indeed strikes a chord also from SciPy point of view.

In my limited opinion, it's not just the dtype but also everything else
that this dtype is going to be used in defines its proper implementation.
Upstream (lapack, fft, so on) and downstream (practically everything else)
does not have generic support on all platforms hence we are providing
convenience to some and taking it away afterwards by not providing the
necessary machinery essentially making it a storage format (similar but
less problematic situation with half) with some valuable and convenient
arithmetics on top.

Hence I'd lean on (B) but I'm curious how we can pull off (C) because that
is going to be tricky for Cythonization etc. which would, probably in
practice, force downcasting to double to do anything with it. Or maybe I'm
missing the context of its utilization.



On Thu, Nov 17, 2022 at 11:12 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-visible dtypes that we have
> to consider; long double types are also baked into the libnpymath static
> library that we ship with NumPy. That's a thing we have to do something
> about anyway (shipping static libraries is not the best idea, see
> https://github.com/numpy/numpy/issues/20880). We just have to make sure
> to not forget about it when thinking about solu

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

[Charles R Harris]

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-visible dtypes that we have
> to consider; long double types are also baked into the libnpymath static
> library that we ship with NumPy. That's a thing we have to do something
> about anyway (shipping static libraries is not the best idea, see
> https://github.com/numpy/numpy/issues/20880). We just have to make sure
> to not forget about it when thinking about solutions here.
>
>
> Potential solutions
> ---
>
> (A) The ideal solution would be to have a proper, portable quad-precision
> (128 bits of precision) dtype. It's now possible to write that externally,
> after all the work that Sebastian and others have put into the dtype
> infrastructure. The dtype itself already exists (
> https://github.com/peytondmurray/quaddtype, maybe there are more
> implementations floating around). It just need the people who have an
> actual need for it to drive that. It's still a significant amount of work,
> so I'll not go into this one more right now.
>
> (B) Full deprecation and removal of all `long double` support from NumPy
> (and SciPy), irrespective of whether the quad dtype comes to life.
>
> Given the current state, I'm personally convinced that that is easily
> justified. However, I know some folks are

[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

[Stephan Hoyer]

On Thu, Nov 17, 2022 at 5:29 PM Scott Ransom  wrote:

> 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
> And see the published paper here, particularly Sec 3.3.1 and footnote #42:
> 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.
>

Hi Scott,

Thanks for sharing your feedback!

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."

Best,
Stephan


>
> Scott
> NANOGrav Chair
> www.nanograv.org
>
>
> --
> 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: sho...@gmail.com
>
___
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

[Charles R Harris]

On Thu, Nov 17, 2022 at 6:30 PM Scott Ransom  wrote:

>
>
> 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.
> 
> > 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
> And see the published paper here, particularly Sec 3.3.1 and footnote #42:
> 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
>
>
>
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?

Chuck
___
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

[James Tocknell]

Around 2015-2016 I created a fork of numpy with IEEE 128bit float
support, but I never upstreamed it because I couldn't see a way to
make it not depend on GCC/libquadmath (the licensing issue Chuck
brought up). I wonder if it's possible now with the dtype changes to
have different dtypes on different platforms (it didn't appear to be
the case when I looked then), which would avoid the need to distribute
third-party libraries to cover certain platforms.

Regards
James

On Fri, 18 Nov 2022 at 12:28, Scott Ransom  wrote:
>
>
>
> 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 
> > 

[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>
 > >> 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 
And see the published paper here, particularly Sec 3.3.1 and footnote #42:
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 



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

[Mike McCarty]

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.

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.


Ralf - Are there other use cases that you have already identified?

Best,
Mike



Mike McCarty
Software Engineering Manager
NVIDIA

From: Scott Ransom 
Sent: Thursday, November 17, 2022 9:04 PM
To: numpy-discussion@python.org 
Subject: [Numpy-discussion] Re: status of long double support and what to do 
about it

External email: Use caution opening links or attachments


On 11/17/22 8:53 PM, Charles R Harris wrote:
>
>
> On Thu, Nov 17, 2022 at 6:30 PM Scott Ransom  > wrote:
>
>
>
> 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.
> 
>  > 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://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnanograv%2FPINT&data=05%7C01%7Cmmccarty%40nvidia.com%7C6bc649fe7b1047fe14d108dac909c126%7C43083d15727340c1b7db39efd9ccc17a%7C0%7C0%7C638043340985084225%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=VgjqmdHsQmp2K%2Bqil47D2JD9h5zWIHbUJD%2Fle8MqxgM%3D&reserved=0
>  
> 
> And see the published paper here, particularly Sec 3.3.1 and footnote #42:
> 
> https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fui.adsabs.harvard.edu%2Fabs%2F2021ApJ...911...45L%2Fabstract&data=05%7C01%7Cmmccarty%40nvidia.com%7C6bc649fe7b1047fe14d108dac909c126%7C43083d15727340c1b7db