Opened 6 years ago
Last modified 6 years ago
#20438 new defect
Wrapper_cdf.call_c() drops imaginary part
Reported by: | cswiercz | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-7.2 |
Component: | cython | Keywords: | fast_callable, interpreters, wrapper_cdf |
Cc: | Merged in: | ||
Authors: | cswiercz | Reviewers: | |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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.)
Attachments (1)
Change History (9)
Changed 6 years ago by
comment:1 Changed 6 years ago by
comment:2 Changed 6 years ago by
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)
comment:3 Changed 6 years ago by
- Description modified (diff)
comment:4 Changed 6 years ago by
- Description modified (diff)
comment:5 Changed 6 years ago by
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]
comment:6 follow-up: ↓ 7 Changed 6 years ago by
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
.
comment:7 in reply to: ↑ 6 Changed 6 years ago by
Replying to cswiercz:
However, I still don't understand why
double_complex
is defined to be a double in the generate code frominterpreters.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
comment:8 Changed 6 years ago by
Apologies for the false alarm. Thank you for the fix suggestion.
Can you just put your Sage code on the ticket instead of linking to a SMC notebook?