Current mood:
Entry tags:	gsl, lisp

Passing a complex number to a foreign function

GSL has a number of functions that take complex scalars as arguments. It defines them with either

typedef struct
  {
    double dat[2];
  }
gsl_complex;

or

typedef struct
  {
    float dat[2];
  }
gsl_complex_float;

These complex numbers are always passed by value, not as a pointer. This presents a problem to Lisp users as many foreign function interfaces and CFFI do not allow passing structs directly.

I thought maybe there was a way around this. Could it be that the real and imaginary parts are simply passed sequentially as if they were two successive arguments to the function? In fact, this does work for numbers of type (complex double-float), but when I tried that on (complex single-float), it would return without error but with the result computed as if the imaginary part was 0. I thought maybe there was a padding issue, so I packed the two single-floats into a single 64-bit word that I then interpreted as a double-float. When passing this single argument to a function that expects a (complex single-float), the correct answer is computed.

For example, the BLAS function caxpy works on single-floats:

(FOREIGN-FUNCALL "gsl_blas_caxpy" 
		 :DOUBLE (PACK-COMPLEX-AS-DOUBLE ALPHA) 
		 :POINTER (MPOINTER X) :POINTER (MPOINTER Y)
		 :INT)

where #'pack-complex-as-double does the packing on the complex scalar alpha described, and the function zaxpy works on double-floats

(FOREIGN-FUNCALL "gsl_blas_zaxpy" 
		 :DOUBLE (REALPART ALPHA) 
		 :DOUBLE (IMAGPART ALPHA) 
		 :POINTER (MPOINTER X) :POINTER (MPOINTER Y) 
		 :INT)

This has now been incorporated into GSLL. I assume it is not portable; I am on Debian amd64 but it works on x86 (32 bit) as well. It seems to be portable across CL implementations, at least it works in SBCL and CCL. So I have built in a test that tries a simple function; if it gets the wrong answer than all the functions that want a complex scalar will signal an error.