Wrapper_cdf.call_c() drops imaginary part

[cswiercz]

A Wrapper_cdf object is created by fast_callable(f, vars=[], domain=None) when the domain is set to CDF. The standard method of evaluating the function represented by the resulting object by simple calling it behaves normally.

However, the Cython-accesible method Wrapper_cdf.call_c(double_complex* out, double_complex* in) drops the imaginary component of the input and output. A demonstration of this behavior:

cython("""
from sage.ext.interpreters.wrapper_cdf cimport Wrapper_cdf

# The following is the self-proclaimed "hack" in the auto-gened file
# ext/interpreters/wrapper_cdf.pyx
cdef extern from "complex.h":
    ctypedef double double_complex "double complex"  # cannot just change "double" to "complex"
    cdef double creal(double_complex)
    cdef double cimag(double_complex)
    cdef double_complex _Complex_I

cdef inline complex call_c(Wrapper_cdf f, complex x, complex y):
    # Quickly evaluates a fast_callable function `f` on domain `CDF`
    # (a Wrapper_cdf tupe) with complex arguments `x`, and `y`.

    cdef double_complex[2] input
    cdef double_complex* output = [0]
    cdef complex value
    input[0] = <double_complex>x
    input[1] = <double_complex>y
    f.call_c(input, output)
    value = <complex>(output[0])
    return value

from sage.all import QQ, CDF, fast_callable

R = QQ['x,y']
x,y = R.gens()
f = y**3 - 2*x**3*y + x**7
f_fast = fast_callable(f, vars=[x,y], domain=CDF)

cdef complex a = 1. + 1.j
cdef complex b = 1. - 2.j

cdef complex value_sage = f_fast(a,b)
cdef complex value_c = call_c(f_fast,a,b)

print 'sage value:  ', value_sage
print 'call_c value:', value_c
""")

sage value:   (-7-18j)
call_c value: (-7+0j)

The source code for Wrapper_cdf is generated by sage/src/sage_setup/autogen/interpreters.py. Within, there are some peculiar comments when it comes to defining the custom data type double_complex. The following is from ​Line 2533 where the Cython header for Wrapper_cdf is defined.

# We need the type double_complex to work around
#   http://trac.cython.org/ticket/869
# so this is a bit hackish.
cdef extern from "complex.h":
    ctypedef double double_complex "double complex"

But this seems to define double_complex as a double, which explains why the imaginary part is dropped. Since complex externs are not allowed in Cython(?) I tried changing this line to

ctypedef double complex double_complex "double complex"

but the following compilation error is raised

[1/3] Cythonizing sage/ext/interpreters/wrapper_cdf.pyx

Error compiling Cython file:
------------------------------------------------------------
...
cdef extern from "interp_cdf.h":
    double_complex interp_cdf(double_complex* args,
        double_complex* constants,
        PyObject** py_constants,
        double_complex* stack,
        int* code) except? -1094648119105371
                           ^
------------------------------------------------------------

sage/ext/interpreters/wrapper_cdf.pyx:64:28: Not allowed in a constant expression

This part of the wrapper code is generated by ​interpreters.py:3290. Commenting out this except? part resolves the compile issue as well as the missing imaginary part problems. Attached is a patch which makes the appropriate double_complex fixes as well as, probably erroneously, comments out the except so that the code compiles.

(Motivation: I would like to directly call the underlying C function of a Wrapper_cdf for performance purposes.)

[jdemeyer]

Can you just put your Sage code on the ticket instead of linking to a SMC notebook?

[cswiercz]

Of course! Below is the Cython input.

# input
# is there no Sage %%cython cell magic?

cython("""
from sage.ext.interpreters.wrapper_cdf cimport Wrapper_cdf

# cannot do below: it redefines double_complex from the
# auto-generated headers created by ext.interpreters.py. generates
# cython error:
#
# "Cannot assign type 'double complex [2]' to 'double_complex *'"
#
#ctypedef complex double_complex "double complex"

# Instead, use the (incorrect) definition as given in the auto-gen
# wrappers. This ctypedefs "double_complex" as a "double"! That is,
# all complex parts of calculations are lost when invoking 
# Wrapper_cdf.call_c()!
#
# The following is the self-proclaimed "hack" in the auto-gened file
# ext/interpreters/wrapper_cdf.pyx
cdef extern from "complex.h":
    ctypedef double double_complex "double complex"  # cannot just change "double" to "complex"
    cdef double creal(double_complex)
    cdef double cimag(double_complex)
    cdef double_complex _Complex_I

cdef inline complex call_c(Wrapper_cdf f, complex x, complex y):
    # Quickly evaluates a fast_callable function `f` on domain `CDF`
    # (a Wrapper_cdf tupe) with complex arguments `x`, and `y`.

    cdef double_complex[2] input
    cdef double_complex* output = [0]
    cdef complex value
    input[0] = <double_complex>x
    input[1] = <double_complex>y
    f.call_c(input, output)
    value = <complex>(output[0])
    return value

from sage.all import QQ, CDF, fast_callable

R = QQ['x,y']
x,y = R.gens()
f = y**3 - 2*x**3*y + x**7
f_fast = fast_callable(f, vars=[x,y], domain=CDF)

cdef complex a = 1. + 1.j
cdef complex b = 1. - 2.j

cdef complex value_sage = f_fast(a,b)
cdef complex value_c = call_c(f_fast,a,b)

print 'sage value:  ', value_sage
print 'call_c value:', value_c
""")

And the corresponding output of the last two print statements:

sage value:   (-7-18j)
call_c value: (-7+0j)

[jdemeyer]

A workaround:

cdef inline complex call_c(Wrapper_cdf f, complex x, complex y):
    # Quickly evaluates a fast_callable function `f` on domain `CDF`
    # (a Wrapper_cdf tupe) with complex arguments `x`, and `y`.

    cdef double complex[2] input
    cdef double complex[1] output
    input[0] = x
    input[1] = y
    f.call_c(<double_complex*>input, <double_complex*>output)
    return output[0]

[cswiercz]

Huh. That's interesting. Thank you for the fix. At least, it gives me a way to solve my current problem.

However, I still don't understand why double_complex is defined to be a double in the generate code from interpreters.py. Would this still be considered something that needed fixing? After all, calls to fast_callable are somehow unaffected by this ctypedef.

[jdemeyer]

Replying to cswiercz:

However, I still don't understand why double_complex is defined to be a double in the generate code from interpreters.py.

You answered your own question in the ticket description:

# We need the type double_complex to work around
#   http://trac.cython.org/ticket/869

[cswiercz]

Apologies for the false alarm. Thank you for the fix suggestion.