How to use numpy longdouble dtype?

Question

I am trying to work with np.longdouble dtype in my Python code and trying to use NumPy to control the long doubles that I get from the C module compiled with Cython.

Suppose I do this:

 import numpy as np print np.finfo(np.longdouble) Machine parameters for float128 --------------------------------------------------------------------- precision= 18 resolution= 1e-18 machep= -63 eps= 1.08420217249e-19 negep = -64 epsneg= 5.42101086243e-20 minexp=-16382 tiny= 3.36210314311e-4932 maxexp= 16384 max= 1.18973149536e+4932 nexp = 15 min= -max --------------------------------------------------------------------- a = np.longdouble(1e+346) a Out[4]: inf b = np.longdouble(1e+347) b Out[6]: inf c = a/b /usr/lib/python2.7/site-packages/spyderlib/widgets/externalshell/start_ipython_kernel.py:1: RuntimeWarning: invalid value encountered in longdouble_scalars # -*- coding: utf-8 -*- c Out[8]: nan a.dtype, b.dtype, c.dtype Out[9]: (dtype('float128'), dtype('float128'), dtype('float128')) 

In essence, this is related to the same question as in this question, and I understand that Python will first convert 1e+346 to a float, whose representation will be inf . However, can anyone suggest a workaround? Is there a way to create NumPy longdoubles that don't convert to float at first?

I have a C module that can output long doubles that I want to use in a numpy dtype np.longdouble .

Even if the solution involves re-compiling Python / NumPy, I'm ready to give it a try.

Answer 1

There are a few things you can consider.

The first is a mess. NumPy knows longdouble and float128 . Unfortunately, the names are misleading, the basic implementation is C long double, which is usually (but not always) 80-bit. (Actually, you can see it here by looking at “accuracy”, 18 digits are approximately 60 bits, and an 80-bit float has 64 bits in the mantissa. The accuracy will be about 34 digits if real 128-bit floats were used.)

There can be no direct way to pass long doubles as arguments to a C function, but if you pass pointers instead, you can avoid the problem. For example, yoy can pass the data of your array as uint8 (using myarray.view(dtype='uint8') ) and discarding the buffer pointer in long double * in your C program. At least then Python has nothing to do with conversions types. (Most likely, you do not need to accept view , because in the end you only export the pointer to the array buffer.)

Just make sure that this trick relies on a compiler that has the same type settings when compiling Python and your C program. In addition to differences in accuracy, there may be differences in byte orders (rarely, if programs run on the same computer) and differences in alignment. My Python seems to align longdouble elements with 16-byte boundaries (i.e. there are always 6 bytes of zeros per element), but the C compiler can use 10/12/16 byte alignment.

As far as I know, the details are implementation specific. Thus, this is doable, but requires some extra care, and there may be portability issues.