Investigate use with numpy-quaddtype

@SwayamInSync

Hi!

I'm one of the NumPy maintainers. I'm also an engineer @Quansight-Labs. Last year as part of our internship program @SwayamInSync implemented a new 128 bit floating point data type for NumPy. See his blog post for more detail. You can also see the code here.

Our plan with the quaddtype is to provide a way for users of numpy.longdouble to have access to both a fast hardware 80 bit float as well as a float128 implementation provided by the SLEEF library. You can choose which "backend" you want to use at runtime.

We still have a bit more work to do to improve the packaging of numpy-quaddtype, but I'm also curious if what we have now is generally useful for pulsar timing work.

@SwayamInSync

Thank you for your comments. We have been struggling with this issue for a while (see #1170). I did see the work you and @SwayamInSync did, and was very interested in it, but when I looked at the page I admit I was scared off by "These dtypes are not meant for real-world use yet." I did check back occasionally but didn't see any dramatic changes. But if you think it's ready for testing I will do that.

OK, I just started to test this but couldn't even get the suggested commands at https://github.com/numpy/numpy-user-dtypes/tree/main/quaddtype to run, or the packaged tests. I assume I must be doing something wrong since it seems there are others who have used this correctly. I can file an issue over on https://github.com/numpy/numpy-user-dtypes/tree/main/quaddtype if that's preferred, or if you have any suggestions here I would be grateful.

On a M2 Mac with x64 conda:

created new environment with python=3.12
installed numpy, as well as astropy and a few others through conda
pip install numpy-quaddtype
Then:

In [1]: from numpy_quaddtype import QuadPrecDType, QuadPrecision
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 from numpy_quaddtype import QuadPrecDType, QuadPrecision

File ~/pythonpackages/numpy-user-dtypes/quaddtype/numpy_quaddtype/__init__.py:1
----> 1 from ._quaddtype_main import (
      2     QuadPrecision,
      3     QuadPrecDType,
      4     is_longdouble_128,
      5     get_sleef_constant,
      6     set_num_threads,
      7     get_num_threads,
      8     get_quadblas_version
      9 )
     11 __all__ = [
     12     'QuadPrecision', 'QuadPrecDType', 'SleefQuadPrecision', 'LongDoubleQuadPrecision',
     13     'SleefQuadPrecDType', 'LongDoubleQuadPrecDType', 'is_longdouble_128', 
   (...)     18     'set_num_threads', 'get_num_threads', 'get_quadblas_version'
     19 ]
     21 def SleefQuadPrecision(value):

ModuleNotFoundError: No module named 'numpy_quaddtype._quaddtype_main'

Or if I run the tests in the GitHub package, I get 779 failures (don't want to list them all). Is there something wrong with this installation procedure?

@dlakaplan

Hey @dlakaplan
You can follow the steps present here https://github.com/numpy/numpy-user-dtypes/tree/main/quaddtype#installation-from-source and please install from the source, pip version is pretty old (we will update it with the upcomong v1.0 release)

The incorrect in above step it seems from the output is that you installed the very older version from pip but since within the folder the python is searching APIs from the current project source

Looking at the issue you mentioned, it seems we have pretty much to be work as the alternative for xprec/ld mentioned in issue.
Let me know I can help in building a patch with numpy-quaddtype to test for that issue.

Thanks for the response. I am attempting to install from source but it's not entirely straightforward. So far I can get sleef installed but still having trouble with quaddtype itself.

For our user base having a package installation (pip and/or conda) is important, but it sounds like that's in your plans. Will that package include sleef as well?

Will that package include sleef as well?

Yes even the version uploaded currently of PyPI comes with included sleef. So we will be continuing with this it's just current source code has multiple better APIs and support so that's why I am recommending it.

For our user base having a package installation (pip and/or conda) is important, but it sounds like that's in your plans.

There's a package on conda-forge too: https://anaconda.org/conda-forge/numpy_quaddtype

@SwayamInSync

@SwayamInSync did a bunch of work the last couple weeks to simplify the packaging situation. Maybe you can try again?

Indeed, I could install it easily and run the sample code. Thanks! It will take longer to see how well this works for us. I tried to just look at the reported finfo but it seems like that's not available, and that you were aware of this (numpy/numpy#27231). So we may need to adjust our tests, but that should be doable.

Any suggestions for how to quickly test for the available precision without finfo? We had a number of tests that used finfo.eps.

any suggestions for how to quickly test for the available precision without finfo? We had a number of tests that used finfo.eps.

We internally package this within the numpy_quaddtype, it should be as simple as

In [1]: import numpy_quaddtype as npq

In [2]: npq.precision
Out[2]: 33

In [3]: npq.epsilon
Out[3]: QuadPrecision('1.925929944387235853055977942584927e-034', backend='sleef')

We're aware that np.finfo doesn't know about user dtypes, it's on our list of things to work on...

OK, I finally got around to installing this on a machine under the arm platform. I used a fresh conda install and a fresh environment, with numpy-quaddtype installed through pip. But when I try some of the basic tests I get:

>>> import numpy_quaddtype as npq
OMP: Error numpy/numpy-user-dtypes#15: Initializing libomp.dylib, but found libomp.dylib already initialized.
OMP: Hint This means that multiple copies of the OpenMP runtime have been linked into the program. That is dangerous, since it can degrade performance or cause incorrect results. The best thing to do is to ensure that only a single OpenMP runtime is linked into the process, e.g. by avoiding static linking of the OpenMP runtime in any library. As an unsafe, unsupported, undocumented workaround you can set the environment variable KMP_DUPLICATE_LIB_OK=TRUE to allow the program to continue to execute, but that may cause crashes or silently produce incorrect results. For more information, please see http://openmp.llvm.org/

I don't know if this is something specific to the machine or to how numpy-quaddtype was packaged. Any ideas?

I've never seem this before but IMO getting an error like this on import after installing a wheel is bu definition a bug in the quaddtype packaging. Would you mind filing an issue on the numpy-user-dtypes repo about this? It would also help to get information about your system and how you set up your Python environment so we can reproduce the issue.

OK, thanks. numpy/numpy-quaddtype#22

OK, moving on now. One possible way to implement this is to monkey-patch np.longdouble to QuadPrecDType(backend="sleef"). But having a callable function of the dtype would help with that. Is that something that's in the works? I don't know the tradeoffs.

Long-term, we probably would need to write a NEP to add the dtype to numpy itself. Out of curiosity, why is monkeypatching numpy needed on your end? Is there a lot of code that uses string type codes or something like that?

@ngoldbaum

Hi, @ngoldbaum and @SwayamInSync One issue I'm encountering now is that we store our data in astropy.table.Table objects and that works fine, but it fails when we try to stack two tables:

import numpy as np
from numpy_quaddtype import QuadPrecDType, QuadPrecision
from astropy.table import Table, Column, vstack

x1 = np.array(np.arange(10), dtype=QuadPrecDType)
t1 = Table([Column(x1, name="x")])
print(f"x.dtype={t1['x'].dtype}")
x2 = np.array(np.arange(10), dtype=QuadPrecDType)
t2 = Table([Column(x2, name="x")])
t = vstack([t1, t2])

so the last operation fails with:

    515     if data is None:
--> 516         self_data = np.zeros((length,) + shape, dtype=dtype)
    517     elif isinstance(data, BaseColumn) and hasattr(data, "_name"):
    518         # When unpickling a MaskedColumn, ``data`` will be a bare
    519         # BaseColumn with none of the expected attributes.  In this case
    520         # do NOT execute this block which initializes from ``data``
    521         # attributes.
    522         self_data = np.array(data.data, dtype=dtype, copy=copy)

TypeError: data type "QuadPrecDType(backend='sleef')" not understood

I don't know why this is happening here, since I can clearly create new arrays with that dtype. A workaround may be possible, but in case there was something obvious/subtle that I was missing, I thought I'd ask. Note that this is with astropy 8.0.

@mhvk

I'm not sure offhand.

@mhvk do you happen to have any insight about the interaction between astropy and quaddtype described in the comment above? Hopefully Marten has a bit more context about this.

Interesting conversation! I've been meaning to try the quaddtype for astropy's Time (which currently uses two doubles), so nice to see the progress.

Directly related to here: The error message suggests it is not related to astropy since it is just trying to create a zero-filled array with np.zeros... (would check, but my office machine seems to old for numpy-quaddtype, I get an error on import, RuntimeError: Invalid slot with value 1031 passed in.)

@mhvk

Hi, @mhvk Nice to hear from you. I understand your point about this not being astropy related, but something like np.zeros((40,),dtype=QuadPrecDType) works just fine. So I don't know if somehow it's not seeing the import of that package. In my example I am creating arrays in that way to start.

@dlakaplan

Ah I missed this was just using numpy.

@dlakaplan if you can make a more minimal example that triggers this behavior it might help to figure out the issue.

Unfortunately I cannot figure out a more minimal example. That is my problem. Array creation on its own works fine, but somehow when this is done from within the table module it does not work. So I don't know if it's a subtle namespace issue or something else. If you have ideas about things to try I can do so.

Can you at least share the full traceback in e.g. a github gist? It's not clear to me where the snippet of your traceback comes from. Astropy internals?

It might also help to put a breakpoint right above the problematic np.zeros call to see exactly what dtype is (e.g. its type). Maybe the repr is fooling us?

Ah its the issue from astropy, they are storing the dtype as a string and then passing that string to create zeros, basically they are doing something as following

In [3]: np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")

TypeError: data type "QuadPrecDType(backend='sleef')" not understood

They should just be storing it like a dtype object rather than string

Here this line is the issue https://github.com/astropy/astropy/blob/c0ef822afb1fee1e40f6a980ba6b64f9c1486863/astropy/utils/metadata/utils.py#L59

Generally in numpy string parsing of dtype works, but this is a custom dtype, I am not sure whether there's a way to make custom dtype also get parsed in NumPy but TBH I don't like this method anyways

@SwayamInSync

Ah, yes, I see. It sounds like @SwayamInSync identified the issue, but in case it's helpful the full traceback is:

---> 10 t = vstack([t1, t2])

File ~/miniforge3/envs/nptest2/lib/python3.13/site-packages/astropy/table/operations.py:715, in vstack(tables, join_type, metadata_conflicts)
    712 if len(tables) == 1:
    713     return tables[0]  # no point in stacking a single table
--> 715 out = _vstack(tables, join_type, metadata_conflicts)
    717 # Merge table metadata
    718 _merge_table_meta(out, tables, metadata_conflicts=metadata_conflicts)

File ~/miniforge3/envs/nptest2/lib/python3.13/site-packages/astropy/table/operations.py:1492, in _vstack(arrays, join_type, metadata_conflicts)
   1488     raise NotImplementedError(
   1489         f"vstack unavailable for mixin column type(s): {col_cls.__name__}"
   1490     )
   1491 try:
-> 1492     col = col_cls.info.new_like(cols, n_rows, metadata_conflicts, out_name)
   1493 except metadata.MergeConflictError as err:
   1494     # Beautify the error message when we are trying to merge columns with incompatible
   1495     # types by including the name of the columns that originated the error.
   1496     raise TableMergeError(
   1497         f"The '{out_name}' columns have incompatible types: "
   1498         f"{err._incompat_types}"
   1499     ) from err

File ~/miniforge3/envs/nptest2/lib/python3.13/site-packages/astropy/table/column.py:482, in ColumnInfo.new_like(self, cols, length, metadata_conflicts, name)
    454 """
    455 Return a new Column instance which is consistent with the
    456 input ``cols`` and has ``length`` rows.
   (...)    476 
    477 """
    478 attrs = self.merge_cols_attributes(
    479     cols, metadata_conflicts, name, ("meta", "unit", "format", "description")
    480 )
--> 482 return self._parent_cls(length=length, **attrs)

File ~/miniforge3/envs/nptest2/lib/python3.13/site-packages/astropy/table/column.py:1246, in Column.__new__(cls, data, name, dtype, shape, length, description, unit, format, meta, copy, copy_indices)
   1241 if isinstance(data, MaskedColumn) and np.any(data.mask):
   1242     raise TypeError(
   1243         "Cannot convert a MaskedColumn with masked value to a Column"
   1244     )
-> 1246 self = super().__new__(
   1247     cls,
   1248     data=data,
   1249     name=name,
   1250     dtype=dtype,
   1251     shape=shape,
   1252     length=length,
   1253     description=description,
   1254     unit=unit,
   1255     format=format,
   1256     meta=meta,
   1257     copy=copy,
   1258     copy_indices=copy_indices,
   1259 )
   1260 return self

File ~/miniforge3/envs/nptest2/lib/python3.13/site-packages/astropy/table/column.py:516, in BaseColumn.__new__(cls, data, name, dtype, shape, length, description, unit, format, meta, copy, copy_indices)
    501 def __new__(
    502     cls,
    503     data=None,
   (...)    513     copy_indices=True,
    514 ):
    515     if data is None:
--> 516         self_data = np.zeros((length,) + shape, dtype=dtype)
    517     elif isinstance(data, BaseColumn) and hasattr(data, "_name"):
    518         # When unpickling a MaskedColumn, ``data`` will be a bare
    519         # BaseColumn with none of the expected attributes.  In this case
    520         # do NOT execute this block which initializes from ``data``
    521         # attributes.
    522         self_data = np.array(data.data, dtype=dtype, copy=copy)

TypeError: data type "QuadPrecDType(backend='sleef')" not understood

@mhvk

Interesting conversation! I've been meaning to try the quaddtype for astropy's Time (which currently uses two doubles), so nice to see the progress.

Directly related to here: The error message suggests it is not related to astropy since it is just trying to create a zero-filled array with np.zeros... (would check, but my office machine seems to old for numpy-quaddtype, I get an error on import, RuntimeError: Invalid slot with value 1031 passed in.)

Hey @mhvk maybe you can raise the issue on quaddtype with some more info, The only requirement should be numpy>=2.4 because there we added an extra slot, so that might be missing if using older version.

Ah its the issue from astropy, they are storing the dtype as a string and then passing that string to create zeros, basically they are doing something as following

In [3]: np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")

TypeError Traceback (most recent call last)
Cell In[3], line 1
----> 1 np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")

TypeError: data type "QuadPrecDType(backend='sleef')" not understood
They should just be storing it like a dtype object rather than string

Sounds like for now I should try a work-around, but should I make an astropy issue? Or will one of you do that?

@dlakaplan

Okay so I was checking the numpy's string parsing of dtypes and found a workaround, @dlakaplan you might need to add this line in top of your code to make that work

np._core.sctypeDict[str(QuadPrecDType())] = QuadPrecDType()

Basically this registers the string variant of numpy_quaddtype into the sctypeDict hashmap. I will check with the team if we can expose a hook and make it work implicitly within the dtype but till then this hack will work

Wow, thanks! that indeed works great.

Thanks all. I raised astropy/astropy#19197 for the astropy side.

Now excited to see how well this works with Quantity too...

Navigation Menu

Repository navigation

Investigate use with numpy-quaddtype #1944

28 remaining items

In [3]: np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Investigate use with numpy-quaddtype #1944

Description

Activity

dlakaplan commented on Aug 25, 2025

dlakaplan commented on Aug 25, 2025

SwayamInSync commented on Aug 25, 2025

dlakaplan commented on Aug 26, 2025

SwayamInSync commented on Aug 26, 2025

ngoldbaum commented on Aug 26, 2025

ngoldbaum commented on Sep 12, 2025

dlakaplan commented on Sep 12, 2025

SwayamInSync commented on Sep 12, 2025

ngoldbaum commented on Sep 12, 2025

dlakaplan commented on Oct 10, 2025

ngoldbaum commented on Oct 10, 2025

dlakaplan commented on Oct 10, 2025

dlakaplan commented on Oct 22, 2025

ngoldbaum commented on Oct 22, 2025

28 remaining items

dlakaplan commented on Jan 20, 2026

ngoldbaum commented on Jan 20, 2026

mhvk commented on Jan 20, 2026

dlakaplan commented on Jan 20, 2026

ngoldbaum commented on Jan 20, 2026

dlakaplan commented on Jan 20, 2026

ngoldbaum commented on Jan 20, 2026

SwayamInSync commented on Jan 20, 2026

SwayamInSync commented on Jan 20, 2026

dlakaplan commented on Jan 20, 2026

SwayamInSync commented on Jan 20, 2026

dlakaplan commented on Jan 20, 2026

In [3]: np.zeros((40,), dtype="QuadPrecDType(backend='sleef')")

SwayamInSync commented on Jan 20, 2026

dlakaplan commented on Jan 20, 2026

mhvk commented on Jan 20, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions