tuwien-cms / xprec

Full quadruple precision (double-double) data type for numpy
MIT License
12 stars 4 forks source link

Using ddouble as a constructor #9

Closed vallis closed 2 years ago

vallis commented 2 years ago

I'm interested in using xprec.ddouble as a drop-in replacement for numpy.longdouble in a Python/numpy library that needs extended precision, and therefore won't run on Apple arm64 (where numpy.longdouble is just float64).

That currently won't work because xprec.ddouble can't be used as a constructor, either from a lower type (e.g., xprec.ddouble(1.0)) or from a string (e.g., xprec.ddouble("3.141592653589793238462643383279")). Would that be difficult to add?

Could I do it myself with a derived class, or are there broader considerations that make it hard? Are there other harder aspects to achieving feature parity with numpy.longdouble?

mwallerb commented 2 years ago

Hi, @vallis!

Thanks for reaching out!

The weird thing about numpy.float64 is that it is not actually a dtype object, but a function. xprec.ddouble instead is a dtype object. A dtype object my_dtype does not support conversion directly, one instead has to do my_dtype.type(1.0) or similar. Admittedly, this is an inconsistency ... the proper way to fix this is to elevate xprec.ddouble to whatever kind of object numpy.float64 is...

As a workaround, it should be possible to do my_dtype.type(1.0) where my_dtype is either xprec.ddouble or np.dtype(numpy.float64).

vallis commented 2 years ago

Thanks Markus!

I see... in my version of numpy at least, float64 is actually a class, that's why it works as a constructor. The corresponding dtype is np.dtype('f8'), which defines the method type.

OK, for my purposes I could then monkeypatch np.longdouble so that it returns ddouble objects on some platforms.

What about conversion from string to ddouble? Where could I get that?

mwallerb commented 2 years ago

Hm, conversion from and to strings is not implemented yet but would be very much appreciated... would you like to have a go?

mwallerb commented 2 years ago

As for the workaround, you could do:

HAVE_EXT_PREC = np.finfo(np.longdouble).nmant >= 63
if HAVE_EXT_PREC:
    ext_dtype = np.dtype(np.longdouble)
else:
    import xprec
    ext_dtype = xprec.ddouble

ext_one = ext_dtype.type(1)
vallis commented 2 years ago

To convert from string, I suppose I could use qd_real::read() from the QD package, and I will need some version of PyUnicode_Decode first. But where does the conversion go in your code? Inside PyDDouble_Cast()? Somewhere else also?

vallis commented 2 years ago

In other news, is there a way to make np.finfo work on xprec.ddouble?

mwallerb commented 2 years ago

In other news, is there a way to make np.finfo work on xprec.ddouble?

Unfortunately, no. numpy provides no way to register information with np.finfo(...). As a workaround, one can use xprec.finfo(...) as drop-in replacement, but this won't work with external libraries ...

mwallerb commented 2 years ago

Closing this for lack of activity ... feel free to reopen.