long double returns and ctypes

Question

i have a c function which returns a long double. i'd like to call this function from python using ctypes, and it mostly works. setting so.func.restype = c_longdouble does the trick -- except that python's float type is a c_double so if the returned value is larger than a double, but well within the bounds of a long double, python still gets inf as the return value. i'm on a 64 bit processor and sizeof(long double) is 16.

any ideas on getting around this (e.g. using the decimal class or numpy) without modifying the c code?

Answer 1

I'm not sure you can do it without modifying the C code. ctypes seems to have really bad support for long doubles - you can't manipulate them like numbers at all, all you can do is convert them back and forth between the native float Python type.

You can't even use a byte array as the return value instead of a c_longdouble, because of the ABI - floating-point values aren't returned in the %eax register or on the stack like normal return values, they're passed through the hardware-specific floating-point registers.

Answer 2

If you have a function return a subclass of c_longdouble, it will return the ctypes wrapped field object rather than converting to a python float. You can then extract the bytes from this (with memcpy into a c_char array, for example) or pass the object to another C function for further processing. The snprintf function can format it into a string for printing or conversion into a high-precision python numeric type.

import ctypes
libc = ctypes.cdll['libc.so.6']
libm = ctypes.cdll['libm.so.6']

class my_longdouble(ctypes.c_longdouble):
    def __str__(self):
        size = 100
        buf = (ctypes.c_char * size)()
        libc.snprintf(buf, size, '%.35Le', self)
        return buf[:].rstrip('\0')

powl = libm.powl
powl.restype = my_longdouble
powl.argtypes = [ctypes.c_longdouble, ctypes.c_longdouble]

for i in range(1020,1030):
    res = powl(2,i)
    print '2**'+str(i), '=', str(res)

Output:

2**1020 = 1.12355820928894744233081574424314046e+307
2**1021 = 2.24711641857789488466163148848628092e+307
2**1022 = 4.49423283715578976932326297697256183e+307
2**1023 = 8.98846567431157953864652595394512367e+307
2**1024 = 1.79769313486231590772930519078902473e+308
2**1025 = 3.59538626972463181545861038157804947e+308
2**1026 = 7.19077253944926363091722076315609893e+308
2**1027 = 1.43815450788985272618344415263121979e+309
2**1028 = 2.87630901577970545236688830526243957e+309
2**1029 = 5.75261803155941090473377661052487915e+309

(Note that my estimate of 35 digits of precision turned out to be excessively optimistic for long double calculations on Intel processors, which only have 64 bits of mantissa. You should use %a rather than %e/f/g if you intend to convert to a format that is not based on decimal representation.)

Answer 3

If you need high-precision floating point, have a look at GMPY.

GMPY is a C-coded Python extension module that wraps the GMP library to provide to Python code fast multiprecision arithmetic (integer, rational, and float), random number generation, advanced number-theoretical functions, and more.

GMP contains high-level floating-point arithmetic functions (mpf). This is the GMP function category to use if the C type `double' doesn't give enough precision for an application. There are about 65 functions in this category.