8000 Options for implementing a quadruple-precision dtype · Issue #14574 · numpy/numpy · GitHub
[go: up one dir, main page]
Skip to content

Options for implementing a quadruple-precision dtype #14574

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
2sn opened this issue Sep 22, 2019 · 52 comments
Open

Options for implementing a quadruple-precision dtype #14574

2sn opened this issue Sep 22, 2019 · 52 comments

Comments

@2sn
Copy link
Contributor
2sn commented Sep 22, 2019

Reproducing code example:

import numpy as np
np.exp(8j * np.arctan(np.array(1, dtype=np.float16))).imag                                                              
np.exp(8j * np.arctan(np.array(1, dtype=np.float32))).imag                                                              
np.exp(8j * np.arctan(np.array(1, dtype=np.float64))).imag                                                              
np.exp(8j * np.arctan(np.array(1, dtype=np.float128))).imag

Output:

-0.0019353059714989746
1.7484556000744883e-07
-2.4492935982947064e-16
 1.0033115225336664047e-19

The last one seems inappropriately inexact, more what you would expect from a float80 rather than a float128. (Yes, the exp then operates on comples256 but all variants I tried get the same accuracy.) Is the quad precision internally just using the Intel built-in 80 bit floats rather than real 128 bit arithmetics? If so, why waste the space and not provide a float80 data type? Alignment issues would not be a lot worse than for float16, which is provided.

Maybe, since I cannot find a reference to it at https://docs.scipy.org/doc/numpy/user/basics.types.html it is only experimental at this time?

It would be good if the documentation could elaborate on precision of float128/complex256

I think libraries should be available to do 128 bit precision arithmetics where not supported by hardware, to have most transparent and as expected behaviour for users, even if slow.

Numpy/Python version information:

numpy: 1.17.2
python: 3.7.4
@charris
Copy link
Member
charris commented Sep 22, 2019 •
edited
Loading

AFAIK, the only common platforms that support quad precision are ARM64, POWER9, and SPARC. The rest either make long double and double the same (MSVC) or make long double extended precision (GCC). What platform are you on?

@2sn
Copy link
Contributor Author
2sn commented Sep 23, 2019

Fedora Linux 30, Intel Xeon (current generation). Use gcc, installed version is 9.2.1.

@seberg seberg added the 57 - Close? Issues which may be closable unless discussion continued label Sep 25, 2019
@seberg
Copy link
Member
seberg commented Sep 25, 2019

I feel like we can probably close this? float128 is a misnomer (since it sounds like it was IEEE 128bit float which it is not). But the precision is thus exactly what we have. Deprecating the name might be possible, but has other issues with the file format. We have gh-10288 open about that.

@2sn
Copy link
Contributor Author
2sn commented Sep 25, 2019 via email

@2sn
Copy link
Contributor Author
2sn commented Sep 25, 2019

OK, in gfortran (Version 9.2.1, "GNU Fortran (GCC) 9.2.1 20190827 (Red Hat 9.2.1-1)") I get from

      program test
      print*,aimag(exp(atan(1.q0)*(0.q8,8.q0)))
      end program test

the output

  -1.73436202602475620495940880520867045E-0034

which is the precision it should be.

I would be really good if NumPy could use/implement that as well. It should not be worse than fortran.

@charris
Copy link
Member
charris commented Sep 25, 2019

We need to worry about other compilers than gcc.

@seberg
Copy link
Member
seberg commented Sep 25, 2019

Wrting an extension type outside numpy is certainly doable right now in principle (it will have some limitation, but for almost everything should be fine). For within nump, there would need to be simple fairly general support.

There was a discussion at some point whether we can just switch to higher precision, real IEEE745 quad. It seems to me like the largest issue could be reading old pickle/npy's (although even now they are not compatible between systems)

@charris
Copy link
Member
charris commented Sep 25, 2019

Several year ago I looked into quad precision libraries and found a couple, but none of them had a NumPy compatible license.

@carlkl
8000 Copy link
Member
carlkl commented Jan 4, 2021

Sleef - see also SLEEF: A Portable Vectorized Library of C Standard Mathematical Functions and Sleef.org - recently published additional support for a software implementation of a IEEE754 __float128 alike datatype with implementations of C99 trigonometric functions (both scalar and vectorized) under the permissive boost license. The initial issue with some implementation details can be found here: Add a vectorized quadruple precision math library.

Sleef supports different CPU architectures and compiler toolchains - see Sleef.org and has a numpy compatible licence.

Another interesting bit is the MPLAPACK, the multiprecision linear algebra package library (formerly known as MPACK). This library can be compiled with __float128 as backend.

In summary, the necessary open source components are available to provide full support for quadrupole precision.

As usual, integrating all of this together into numpy will require some effort.

@charris
Copy link
Member
charris commented Jan 4, 2021

@carlkl Thanks for the heads up, SLEEF looks interesting. @Qiyu8 @seiko2plus Ping.

@shibatch
Copy link
shibatch commented Jan 5, 2021

Hello,

As for SLEEF, there are a couple of things to pay attention.

  1. Beware of this bug : https://bugs.llvm.org/show_bug.cgi?id=47665 It seems that there are other bugs on clang, but I cannot reliably reproduce the problems. GCC seems okay, but gcc 10.2 cannot compile the quad library with SVE support due to ICE.
  2. Arithmetic functions (add, sub, mul, and div) do not give perfectly correctly-rounded results. They have very small amount of error. Please tell me if this is problematic.

@Qiyu8
Copy link
Member
Qiyu8 commented Jan 5, 2021

@shibatch I didn't see your scalar version quadruple precision implementation, the vectorized version has an error bound like 0.5000000001 ULP, which is unbearable in numpy's precision standard, maybe scalar version can overcame this flaw.

@shibatch
Copy link
shibatch commented Jan 5, 2021

It is not so hard to just implement a correctly-rounded version of arithmetic functions. The problem is that there is a tradeoff between speed and accuracy. If 0.5000000001 ULP is not bearable, I will add correctly-rounded but slower functions.

@carlkl
Copy link
Member
carlkl commented Jan 5, 2021 •
edited
Loading

Perhaps arithmetic functions used for the trigonometric functions do not necessarily have to be correctly-rounded, only the arithmetic functions exposed to the user for IEEE 754 conformity.
In this context, the question arises as to whether GNU libquadmath arithmetic functions are correctly-rounded or not.

@shibatch
Copy link
shibatch commented Jan 5, 2021

Arithmetic functions are implemented within glibc, and they are correctly rounded.

@carlkl
Copy link
Member
carlkl commented Jan 5, 2021

Paul Zimmermann (INRIA) recently compared the accuracy of several math libraries: accuracy.pdf. This includes a comparison of quadrupole precision between glibc and icc. See also his talk How Slow is Quadruple Precision on Variable Precision in Mathematical and Scientific Computing.

@shibatch
Copy link
shibatch commented Jan 5, 2021

Quadruple precision arithmetic functions are slow because packing to/unpacking from IEEE 754 format are pretty complicated. If computation speed is important, double-double is much faster. SLEEF quad library is designed in a way that the frontend and backend are separated. The frontend does packing and unpacking. Computation by the backend is all carried out with triple-double operations. Thus, it is easy to add another frontend for a different format.

@carlkl
Copy link
Member
carlkl commented Jan 5, 2021

I want to point out my comment #14574 (comment) is motivated by my desire to discuss the need for quadrupole precision in Numpy. It seems to me, that adding quadrupole precision to numpy is feasible in forseable future with the help off Sleef and MPLAPACK for a considerable number of target systems.

Personally I would like to have the possibility to combine the efficiency and elegance of numpy with enhanced precision. On the other side I know that this is a small but important niche within computational science. Usually Fortran or C++ are used there.

I'm not a numpy maintainer, so I have no idea if the effort of integration into numpy pays of. This is the question that has to be discussed first IMHO.

@seberg
Copy link
Member
seberg commented Jan 5, 2021

@carlkl I am doing a lot of slow work to allow you to externally write a quadruple precision that integrates seamlessly. Right now, that is not finished and the parts that are, are not public API yet, but it should fall in places for 1.21 far enough that one can start to use it (experimentally).

There could be some small limitations (at least I assume those won't go away quickly), such as having to use dtype=quadruple.Quadruple instead of dtype="quadruple" (i.e. no string). But it will smoothen out, or at least set the stage, the current issues with externally defined dtypes. That said, even now quaternions work fairly well written outside of NumPy.

My current opinion is to start writing such dtypes outside of NumPy first. While it probably makes some things work less well, it avoids having to worry about partial implementation or stranger corners in NumPy.
I admit, there is the problem that doing it right now with the old API means you need to update it soon, OTOH, I doubt that boilerplate is hard to update. If you wait a few months, it would be a (very simple) trial baloon for a new API.

@charris
Copy link
Member
charris commented Jan 5, 2021

There are uses for quad precision apart from it starting to show up with hardware support. I looked into it because I'd like to use it for a new datetime type. I think astropy uses double double for that, they need the extra precision. As a historical note, VAX had a quad precision type back in the 1980's and I found it useful on occasion.

@seberg
Copy link
Member
seberg commented Jan 5, 2021

I am really not as such opposed with adding things to NumPy, I just feel that comes with a lot of baggage (whether real or mental); and also with stronger than typical backcompat guarantees. Also someone not familiar with NumPy can probably dive into an external implementation much easier (and things tend to work better driven by interested parties).

I agree that having a larger than normal inexact type around does seem like a useful thing (I could imagine even just using it for double checking the accuracy of an algorithm). Datetimes are a great reason, but also one of those complex beasts that, to me, seem better started outside NumPy.

If someone wants to start these projects, I would be happy to put them on the NumPy organization but outside the NumPy code base. For DTypes there is the annoyance of private API being much more powerful and less wrinkly, I admit, hopefully that problem will be gone in a few months.

@carlkl
Copy link
Member
carlkl commented Jan 6, 2021

@seberg I am doing a lot of slow work to allow you to externally write a quadruple precision that integrates seamlessly. Right now, that is not finished and the parts that are, are not public API yet, but it should fall in places for 1.21 far enough that one can start to use it (experimentally)
Is this https://github.com/seberg/numpy-dispatch? BTW. I don't expect any workable results in the near future. Therefore I do not plan a fixed schedule, just thinking about the next steps - probably first testing mplapack with sleef as a backend.

I did a quick look at what the developers of https://github.com/moble/quaternion have done, but at a first glance I didn't get the idea how to extend dtypes with quadrupole-precision. Are there tutorials how to extend numpy with new datatypes?

@seberg
Copy link
Member
seberg commented Jan 6, 2021

No, numpy-dispatch is on the array object level. The idea of that is that you can write code that supports both NumPy and e.g. Jax, etc. arrays. What I am talking about here is the work related to: https://numpy.org/neps/nep-0041-improved-dtype-support.html which will allow to write dtypes outside of NumPy with much fewer issues (also allowing much more complex dtypes than numpy can currently handle reasonably).

There is fairly little information on how to write a dtype currently. And the API has to change pretty much completely. @WarrenWeckesser also has similar tries here: https://github.com/WarrenWeckesser/numtypes which may give a better overview and more explanation.
Its a bit verbose and arcane, but overall it should actually be fairly straight forward. You have to create a new dtype object in C and then register it with NumPy (the main step there is to implement the set of "ArrFunctions") after that you can look into casting and ufuncs. (Pretty much all of that API has to be revamped, but the basic puzzle pieces will remain of course.)

8000
@carlkl
Copy link
Member
carlkl commented Jan 6, 2021

Great, https://github.com/WarrenWeckesser/numtypes has a link to https://github.com/stefanv/teaching, which includes an example for extending with quad types.

@carlkl
Copy link
Member
carlkl commented Jan 7, 2021

When I read Hierarchy and abstract classes from NEP 42, I can read out the following. A quadrupole-precision type should be implementing something like this in future:

Abstract Type:  DType -> Numeric -> Floating -> Float128
Concrete Type:  DType -> Float128

@seberg with future numpy-1.21 or higher is what I assume is meant?

@seberg
Copy link
Member
seberg commented Jan 7, 2021

@carlkl I hope 1.21 is what it means, I would have hoped to be closer to having it done, but thats life. There doesn't need to be an abstract type, Float128 can be a subclass of an abstract Floating, but it doesn't necessarily matter too much to begin with, and those abstract dtypes don't even exist yet.

@aragilar
Copy link
Contributor

@carlkl I had an attempt at adding quad precision floats in https://github.com/aragilar/numpy/ (there's a number of branches as I kept needing to rebase the changes, I don't remember which was the latest one). The issue I ran into was the lack of a quad precision implementation that wasn't LGPL (or didn't depend on a LGPL library along the line - there are some very basic quad precision implementations around which have just the data type plus the basic operation—no exp/sin/cos etc.), if SLEEF solves that aspect, then you may be able to reuse some of the dtype logic I wrote in one of those branches (or use as a guide, depending on how much the numpy internals have changed).

@carlkl carlkl mentioned this issue Nov 13, 2021
Closed
@lesshaste
Copy link

Was there any progress on quad precision support for numpy?

@carlkl
Copy link
Member
carlkl commented Mar 31, 2023

About one year age or so @shibatch ceased the development of sleef. Whether this is still valid I do not know. I haven't seen any progress since then. Sleef itself works without major bugs and is used in some libraries today. However, can we rely on a library that is unmaintained? I guess it will be really hard to maintain the sleef code without his help, but that should not stop anyone from judging's for themselves.

There were also other concerns about using sleef for development of numpy:

  • slow performance for quad precision calculations. double-double is somewhat less precise but much faster.
  • sleef doesn't give perfectly correctly-rounded results for quadmath arithmetik functions.

There are some counter arguments as well:

  • unclear double-double licence of D. Bailey's implementation.
  • no perfectly correctly-rounded results for quadmath arithmetik functions: but does it really matter for qud-precision?
  • there is always a trade off between precision and performance.

All in all, I decided for myself not to pursue this development any further.

BTW: there is now a PR ENH: Add support SLEEF for transcendental functions #23068 see also https://github.com/yamadafuyuka/numpy/tree/add_SLEEF not related to quad precision but related to the integration of the vectorized single/double sleef functions.

@seberg
Copy link
Member
seberg commented Mar 31, 2023

If there is a good basis, starting a project in https://github.com/numpy/numpy-user-dtypes might be cool. That includes an example for a DType start that wraps MPF for multi-precision floats. (I wouldn't consider it usable and obviously this is all using as-of-now unstable/experimental API).

@shibatch
Copy link
shibatch commented Apr 1, 2023 via email

@carlkl
Copy link
Member
carlkl commented Apr 2, 2023

@shibatch, this sounds really promising! I won't have time for such a quadmath-numpy project in the next two months either. After that, let's see. And maybe @yamadafuyuka could comment on the topic sleef-for-numpy too.

@shibatch
Copy link
shibatch commented Apr 3, 2023 via email

@carlkl
Copy link
Member
carlkl commented Apr 3, 2023

@shibatch, this question could be raised in one of the numpy meetings. @charris, what do you think about that?

@carlkl
Copy link
Member
carlkl commented Apr 4, 2023

Will the numpy community adopt SLEEF with the current 0.5+alpha precision arithmetic operators?

@shibatch, how big can alpha become in sleef's quadmath elementary arithmetic operations?

@shibatch
Copy link
shibatch commented Apr 4, 2023 via email

@carlkl
Copy link
Member
carlkl commented Apr 4, 2023

@shibatch, thank you for this information. According to the documentation sleef's quadmath has 0.5000000001 ULP for elementary arithmetic and 1.0 ULP for trigonometric functions.

Does this apply to the entire range? I'm curius because for some double/singe precision sleef trigonometric functions ranges are given in the documenation.

@shibatch
Copy link
shibatch commented Apr 4, 2023 via email

@2sn
Copy link
Contributor Author
2sn commented Apr 4, 2023

@shibatch Sounds excellent.

@carlkl
Copy link
Member
carlkl commented Apr 4, 2023 •
edited
Loading

Yes, this is excellent. And it is hard to believe, that sleefs quadmath 0.5000000001 ULP instead of the IEEE specified < 0.5 ULP really matters in practice. But I leave this question to the experts.

I would really like to see sleef as well sleef quadmath in Paul Zimmmermanns accuracy paper:
https://homepages.loria.fr/PZimmermann/papers/#accuracy

@rgommers rgommers removed the 57 - Close? Issues which may be closable unless discussion continued label Jun 3, 2024
@rgommers rgommers changed the title issues with accuracy of quad precision? Options for implementing a quadruple-precision dtype Jun 3, 2024
@carlkl
Copy link
Member
carlkl commented Aug 8, 2024

Adding Quaddtype from @SwayamInSync

@carlkl
Copy link
Member
carlkl commented Feb 26, 2025

Just curious: as numpy/numpy-user-dtypes#98 is merged - is there a reason to not close this issue?

@mhvk
Copy link
Contributor
mhvk commented Feb 27, 2025

Intrigued by the comment above, I tried the QuadPrecDtype, but it is not ready for prime time yet -- at least I get core dumps even on just adding two numbers together - see numpy/numpy-user-dtypes#105

That said, hopefully that is just a relatively simple glitch and it will provide real support!

Anyway, I think closing this issue would require having a working package and a link to it from the numpy docs.

@shibatch
Copy link
shibatch commented Mar 1, 2025

I believe using tlfloat is much easier.
You will get octuple presicion as a bonus.

@lesshaste
Copy link

Is the proposal to have a true 128 bit float type in numpy? If so that is a great idea. Can anyone explain the status of QuadPrecDtype or if the general idea is still going forwards?

@carlkl
Copy link
Member
carlkl commented Mar 1, 2025

Is the proposal to have a true 128 bit float type in numpy?

Yes

Can anyone explain the status of QuadPrecDtype or if the general idea is still going forwards?

@lesshaste, this is a question for @SwayamInSync. However sleef as well as tlfloat are both external projects created by @shibatch.

His suggestion to use tlfloat in favour of sleef is worth considering.

@seberg
Copy link
Member
seberg commented Mar 2, 2025 •
edited
Loading

numpy-user-dtypes was created to show that you can extend NumPy quite easily outside of the proper repository.
We can always considering moving things into NumPy proper eventually, but for a large chunk, there are few downsides for a dtype living outside of NumPy proper.

So I think the status is still as it was: Improvements to or forks of the quad-dtype in numpy-user-dtypes are very welcome and we can help a bit. We can also provide a repository within the numpy org and "bless" e.g. a numpy-float128 or similar package name.

But, I would still not aim for a true float128 in NumPy. Because that comes with a lot of baggage and things to consider. It might grow into that, but living in NumPy raises expectations and compatibility concerns that are possibly not quick to deal with.

EDIT: And I'll note that I actually like to have prominent dtypes defined outside of NumPy (ml_types is maybe becoming that). Having an ecosystem of those will be great and embracing some defined outside of NumPy will help get more ;).

@shibatch
Copy link
shibatch commented Mar 2, 2025

Solving compatibility issues is one of the goals of tlfloat's development, though.
If you decide to use tlfloat, it is easy to use both true float128 and float256.

@mhvk
Copy link
Contributor
mhvk commented Mar 2, 2025

@shibatch - is there a specific reason you recommend using tlfoat over sleef? Just simplicity of use as it avoids some of the vectorization?

@shibatch
Copy link
shibatch commented Mar 2, 2025

My understanding is that the current implementation on the Numpy side is not vectorized.
If that is the case, tlfloat is faster for basic arithmetic operations.
The accuracy problem I explained above is not present in tlfloat.
Also, tlfloat implements mathematical functions that are not implemented in sleef.

@mhvk
Copy link
Contributor
mhvk commented Mar 2, 2025

OK, thanks, that makes sense! I certainly think that having more functions is more important than vectorized speed, at least initially!

@SwayamInSync
Copy link
SwayamInSync commented Mar 2, 2025 •
edited
Loading

Is the proposal to have a true 128 bit float type in numpy?

Yes

Can anyone explain the status of QuadPrecDtype or if the general idea is still going forwards?

@lesshaste, this is a question for @SwayamInSync. However sleef as well as tlfloat are both external projects created by @shibatch.

His suggestion to use tlfloat in favour of sleef is worth considering.

We do have a initial experimental dtype here numpy-user-dtypes/quaddtype (idea is to have a true 128 bit dtype with backward compatibility of long double) It supports major operations and integrating well with NumPy. But there is still lot of room for improvement and expansion

Agreeing to @shibatch and we mention this in the quaddtype report as well, during development, we experimented with TLFloat as a potential SLEEF replacement due to its IEEE 754 compliance, correctly rounded arithmetic, and 1-ULP accuracy for math functions.

The current major plans regarding quaddtype is to integrate it with long double tests and keep iterating to make it a good long double replacement

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
0