CTypes test failed, complex double problem

[efimov-mikhail]

Bug report

Bug description:

Python binary has been built in the standard way:

./configure --with-pydebug && make -j

But test_ctypes fails:

-> % ./python -m unittest -v test.test_ctypes.test_libc.LibTest.test_csqrt 
test_csqrt (test.test_ctypes.test_libc.LibTest.test_csqrt) ... FAIL

======================================================================
FAIL: test_csqrt (test.test_ctypes.test_libc.LibTest.test_csqrt)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/mikhail.efimov/projects/cpython/Lib/test/test_ctypes/test_libc.py", line 30, in test_csqrt
    self.assertEqual(lib.my_csqrt(4), 2+0j)
    ~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: (5e-324+6.95324111713477e-310j) != (2+0j)

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (failures=1)

System:

-> % cat /etc/issue
Debian GNU/Linux 10 \n \l

GCC:

-> % gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/8/lto-wrapper
OFFLOAD_TARGET_NAMES=nvptx-none
OFFLOAD_TARGET_DEFAULT=1
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Debian 8.3.0-6' --with-bugurl=file:///usr/share/doc/gcc-8/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++ --prefix=/usr --with-gcc-major-version-only --program-suffix=-8 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-libmpx --enable-plugin --enable-default-pie --with-system-zlib --with-target-system-zlib --enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 8.3.0 (Debian 8.3.0-6) 

I've tried to fix it by myself but the result has not been achieved in a reasonable amount of time. There is a simple test I've provided:

-> % cat test_csqrt.c
#include <complex.h>
#include <stdio.h>

int complex my_csqrt(double complex a)
{
    double complex z1 = a;
    double complex z2 = csqrt(a);
    printf("my_csqrt (%.10f%+.10fi) = %.10f%+.10fi\n",
        creal(z1), cimag(z1), creal(z2), cimag(z2));
    return 0;
}

int main() {
    my_csqrt(4.0);
    my_csqrt(4.0+4.0j);
    my_csqrt(-1+0.01j);
    my_csqrt(-1-0.01j);
    return 0;
}

-> % gcc -lm test_csqrt.c -o test_csqrt && ./test_csqrt
my_csqrt (4.0000000000+0.0000000000i) = 2.0000000000+0.0000000000i
my_csqrt (4.0000000000+4.0000000000i) = 2.1973682269+0.9101797211i
my_csqrt (-1.0000000000+0.0100000000i) = 0.0049999375+1.0000124996i
my_csqrt (-1.0000000000-0.0100000000i) = 0.0049999375-1.0000124996i

So, it's not a problem in my libc version.

Moreover, this problem can be reproduced with standard libm.so (/lib/x86_64-linux-gnu/libm-2.28.so):

-> % cat ctypes_fix/test.py      
import ctypes
libm = ctypes.CDLL('libm.so.6')
libm.clog.restype = ctypes.c_double_complex
libm.clog.argtypes = ctypes.c_double_complex,
clog_5 = libm.clog(5.0)
clog_1000_2j = libm.clog(1000.0+2j)
print(f"{clog_5=}")
print(f"{clog_1000_2j=}")

-> % ./python ctypes_fix/test.py
clog_5=(5e-324+6.9529453382261e-310j)
clog_1000_2j=(5e-324+6.9529453382261e-310j)

IMHO, some problem lies in using ctypes.c_double_complex as an argument and return value types. FYI, with double argtype and restype clog works like classical double log:

-> % cat ctypes_fix/test2.py   
import ctypes
libm = ctypes.CDLL('libm.so.6')
libm.clog.restype = ctypes.c_double
libm.clog.argtypes = ctypes.c_double,
clog_5 = libm.clog(5.0)
clog_1000 = libm.clog(1000.0)
print(f"{clog_5=}")
print(f"{clog_1000=}")

-> % ./python ctypes_fix/test2.py
clog_5=1.6094379124341003
clog_1000=6.907755278982137

CPython versions tested on:

CPython main branch

Operating systems tested on:

Linux

Linked PRs

• gh-125322

[Eclips4]

cc @skirpichev

[skirpichev]

Looks to be libffi issue.

Here is what I see.

1. other tests don't fail (i.e. round-trip to/from c_double_complex)
2. @efimov-mikhail, you are on Debian 10 (buster), that has libffi v3.2.1; are you using this version?
3. if so, can you try to link CPython with never libffi (which, probably is available in backports)? Any v3.3+ should be fine.

It seems, that real support for complex types was added on x86_65 in v3.3-rc0. Similar failure was on Sparc during testing my pr, see this. I hoped, it was workarounded.

Edit: yes, they just did it - https://github.com/libffi/libffi/blob/v3.2.1/src/x86/ffitarget.h. Unfortunately, there is no way to detect libffi support for complex types at compilation time. There is even no version define.

[efimov-mikhail]

2. @efimov-mikhail, you are on Debian 10 (buster), that has libffi v3.2.1; are you using this version?
3. if so, can you try to link CPython with never libffi (which, probably is available in backports)? Any v3.3+ should be fine.

The answer is "yes" on both questions. All work fine with libffi-3.3-rc1.

Edit: yes, they just did it - https://github.com/libffi/libffi/blob/v3.2.1/src/x86/ffitarget.h. Unfortunately, there is no way to detect libffi support for complex types at compilation time. There is even no version define.

What do you think about adding some workaround to configure.ac? It can be in form of something like this: Write little code snippet which uses ffi_type_complex_double, calculate clog(1000.0) and exit with error code if creal(clog(1000.0)) less than 5.0. Build this snippet with linking to libffi, test the result of binary execution. And properly define FFI_ACTUALLY_HAS_COMPLEX_TYPE value.

[skirpichev]

What do you think about adding some workaround to configure.ac?

It's a bug. There should be a proper fix to exclude bad library versions. I'm thinking on just

diff --git a/configure.ac b/configure.ac
index 1864e94ace..559b513062 100644
--- a/configure.ac
+++ b/configure.ac
@@ -4043,7 +4043,7 @@ AS_VAR_IF([ac_sys_system], [Darwin], [
   ])
 ])
 AS_VAR_IF([have_libffi], [missing], [
-  PKG_CHECK_MODULES([LIBFFI], [libffi], [have_libffi=yes], [
+  PKG_CHECK_MODULES([LIBFFI], [libffi >= 3.3], [have_libffi=yes], [
     WITH_SAVE_ENV([
       CPPFLAGS="$CPPFLAGS $LIBFFI_CFLAGS"
       LDFLAGS="$LDFLAGS $LIBFFI_LIBS"

v3.3 was released in Nov 2019. Not sure if this is acceptable, however. @Eclips4 ?

Write little code snippet which uses ffi_type_complex_double, calculate clog(1000.0) and exit with error code if creal(clog(1000.0)) less than 5.0.

Or some variant along this road.

I would rather create a test function, that get expected complex input (say 1.25+0.5j) and compares that with local values and then reports result. We should call this with ffi_call() and report result to the configure script.

Let me know if you would like to work on this.

[efimov-mikhail]

It's a bug. There should be a proper fix to exclude bad library versions. I'm thinking on just <ignore libffi with versions less than 3.3>

IMO, there is no need to fully ignore libffi linking with old versions since workaround is simple and issue only relates to complex double/float/long double types. Moreover, supporting of complex numbers without C11 complex types is okay, and making harder assumptions on system libraries just for those types seems to be redundant.

expected complex input (say 1.25+0.5j)

Agree, checking

clog(1.25+0.5j) == (0.2973535538733464+0.3805063771123649j)

seems to be optimal.

Let me know if you would like to work on this.

Yes, I would like to create a PR.

[skirpichev]

Agree, checking clog(1.25+0.5j) == (0.2973535538733464+0.3805063771123649j) seems to be optimal.

Not at all. What I would suggest instead is something like

int
foo (double complex z)
{
    const double complex expected = CMPLX(1.25, -0.5);
    return z == expected;
}

Then ffi_call this with expected argument. All this is actually about passing arguments correctly.

You can try to use instead libm complex functions, but you shouldn't compare computed values exactly.

Yes, I would like to create a PR.

Ok. Feel free to ping me for review.