How to carry out arithmetic operations on mcomplex in LibraryFunction

Question

I want to perform some arithmetic inside the C++ code using mcomplex data type defined in WolframLibrary.h, but the compiler complains that mcomplex is an invalid binary operand for * and +. Question: What do I need to do to carry out arithmetic on mcomplex variables?

Below is the .cc file defining the LibraryFunciton called myFunc.

//mma.cc
#include "WolframLibrary.h"

//Standard LibraryLink functions
DLLEXPORT mint WolframLibrary_getVersion(){
    return WolframLibraryVersion;}

DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
    return 0;}

//User function called through LibraryFunctionLoad
EXTERN_C DLLEXPORT int myFunc(WolframLibraryData libData, mint Argc, 
   MArgument *Args, MArgument Res){

    mcomplex I0 = MArgument_getComplex(Args[0]);
    mcomplex I1 = MArgument_getComplex(Args[1]);

    mcomplex result = (I0 + I1)*(I0 + 2*I1) // <--- this is illegal.  HELP!!!

    MArgument_setComplex(Res,result);

    return LIBRARY_NO_ERROR;
}

xslittlegrass · Accepted Answer · 2017-06-27T03:26:41.027

If you look at the definition of mcomplex in WolframLibrary.h you can see that it is simply a struct with two real numbers:

typedef struct {mreal ri[2];} mcomplex;

There are no + or * been overloaded by mcomplex. So in order for them to work, we need to define them ourselves:

//mma.cc
#include "WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion(){
    return WolframLibraryVersion;}

DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
    return 0;}

mcomplex operator*(mcomplex m1, mcomplex m2){
    mcomplex m3 = {m1.ri[0]*m2.ri[0] - m1.ri[1]*m2.ri[1],
        m1.ri[0]*m2.ri[1]+m1.ri[1]*m2.ri[0]};
    return m3;
}
mcomplex operator*(mreal r1, mcomplex m2){
    mcomplex m3 = {r1*m2.ri[0], r1*m2.ri[1]};
    return m3;
}
mcomplex operator+(mcomplex m1, mcomplex m2){
    mcomplex m3 = {m1.ri[0]+m2.ri[0], m1.ri[1]+m2.ri[1]};
    return m3;
}

//User function called through LibraryFunctionLoad
EXTERN_C DLLEXPORT int myFunc(WolframLibraryData libData, mint Argc, 
   MArgument *Args, MArgument Res){

    mcomplex I0 = MArgument_getComplex(Args[0]);
    mcomplex I1 = MArgument_getComplex(Args[1]);

    mcomplex result = (I0 + I1)*(I0 + 2*I1);

    MArgument_setComplex(Res,result);

    return LIBRARY_NO_ERROR;
}

score 2 · Answer 2 · answered Oct 06 '17 at 12:00

tl;dr Just cast to a C++ std::complex<double>, and compute with it as usual.

mcomplex is defined as

typedef struct {mreal ri[2];} mcomplex;

where mreal is just a double. This means that mcomplex is always two doubles, following each other in memory with no gap.

Since C++11, std::complex<double> is guaranteed to have the same layout, which means that it is safe to reinterpret_cast to it.

std::complex has the usual arithmetic (+, -, *, /) as well as many other functions already implemented for it in the C++ standard library, which makes it very easy to work with.

Here's an example:

/* cpl.cpp */

#include <WolframLibrary.h>
#include <complex>

/* Return the version of Library Link */
EXTERN_C DLLEXPORT mint WolframLibrary_getVersion( ) {
    return WolframLibraryVersion;
}

/* Initialize Library */
EXTERN_C DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
    return LIBRARY_NO_ERROR;
}

/* Uninitialize Library */
EXTERN_C DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
    return;
}

/* Conversions between std::complex<double> and `mcomplex` */

inline const std::complex<double> &from_mcomplex(const mcomplex &z) { return reinterpret_cast<const std::complex<double> &>(z); }
inline std::complex<double> &from_mcomplex(mcomplex &z) { return reinterpret_cast<std::complex<double> &>(z); }

inline const mcomplex &to_mcomplex(const std::complex<double> &z) { return reinterpret_cast<const mcomplex &>(z); }
inline mcomplex &to_mcomplex(std::complex<double> &z) { return reinterpret_cast<mcomplex &>(z); }

EXTERN_C DLLEXPORT int mult(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
    mcomplex a = MArgument_getComplex(Args[0]);
    mcomplex b = MArgument_getComplex(Args[1]);
    MArgument_setComplex(Res, to_mcomplex( from_mcomplex(a) * from_mcomplex(b) ));
    return LIBRARY_NO_ERROR;
}

lib = CreateLibrary[{"cpl.cpp"}, "cpl"];
mult = LibraryFunctionLoad[lib, "mult", {Complex, Complex}, Complex];

mult[I, I]
(* -1. + 0. I *)

LTemplate already has this particular conversion built in. Its mma::complex_t type is just an alias for std::complex<double>.

Example:

// Cpl.h
struct Cpl {
  mma::complex_t mult(mma::complex_t a, mma::complex_t b) {
    return a*b;
  }
};

template = LClass["Cpl",
   {LFun["mult", {Complex, Complex}, Complex]}
  ];

CompileTemplate[template]

LoadTemplate[template]

obj = Make[Cpl]
(* Cpl[1] *)

obj@"mult"[I, I]
(* -1. + 0. I *)

How to carry out arithmetic operations on mcomplex in LibraryFunction

2 Answers2