Math Parser for Complex Numbers in C (ExprTk)

Question

I have been using the ExprTk library quite frequently in the past in order to further process large output files generated with Mathematica (containing mathematical expressions) in C. Until now, I exclusively used this library to process expressions that yield values of the type <double>, for which the library works flawlessly by defining the types

typedef exprtk::symbol_table<double> symbol_table_t;
typedef exprtk::expression<double>     expression_t;
typedef exprtk::parser<double>             parser_t;

and storing "everything" in a struct

struct readInExpression 
{
  double a, b;

  symbol_table_t symbol_table;

  expression_t expression;
};

Reading in a text file that contains the variables a and b as well as e.g. the user-defined function

double my_function(double a, double b) {
    return a+b;
} 

can be achieved by means of

void readInFromFile(readInExpression* f, parser_t* p) {
      std::string file = "xyz.txt";

      std::ifstream ifs(file);
      std::string content( (std::istreambuf_iterator<char>(ifs) ),
                           (std::istreambuf_iterator<char>()    ) );

      f->symbol_table.add_variable("a",f->a);
      f->symbol_table.add_variable("b",f->b);

      f->symbol_table.add_function("my_function",my_function);

      f->expression.register_symbol_table(f->symbol_table);

      p->compile(content,f->expression);
}

One may then evaluate the read-in expression for arbitrary values of a and b by using

double evaluateFile(readInExpression* f, double a, double b) {
    f->a = a;
    f->b = b;

    return f->expression.value();
}

Recently, I ran into problems when trying to process text files that contain complex numbers and functions that return complex values of the type std::complex<double>. More specifically, I have a .txt file that contains expressions of the form

2*m*A0(m*m) + Complex(1.0,2.0)*B0(M*M,0.0,m*m)

where A0(a) and B0(a,b,c) are the scalar loop integrals that arise from the Passarino-Veltman reduction of (tensor) loop integrals in high-energy physics. These can be evaluated numerically in C using LoopTools, where it is to be noted that they take complex values for certain values of a, b, and c. Simply replacing <double> by std::complex<double> in the typedefs above throws tons of errors when compiling. I am not sure whether the ExprTk library is able to handle complex numbers at all -- I know that it cannot deal with custom classes, but from what I understand, it should be able to handle native datatypes (as I found here, ExprTk is able to at least deal with vectors, but given the complexity of the expressions I need to process, I do not think it will be possible to somehow rewrite everything in form of vectors, in particular due to the difference in doing algebra with complex numbers and vectors). Note that neither can I split the expressions into real and imaginary part because I have to evaluate the expressions for many different values of the variables.

Although I dealt with complex numbers and the mentioned functions A0(a) and B0(a,b,c) in text files before, I solved this by simply including the .txt files in C using #include "xyz.txt", implemented in a corresponding function, which, however, seems impossible given the size of the text files at hand (the compiler throws an error if I try to do so).

Does anybody know if and how ExprTk can deal with complex numbers? (A MWE would be highly appreciated.) If that is not the case, can anyone here suggest a different math parser that is user friendly and can deal with complex numbers of the type std::complex<double>, at the same time allowing to define custom functions that themselves return such complex values?

A MWE:

/************/
/* Includes */
/************/

#include <iostream> // Input and Output on the console
#include <fstream> // Input and Output to files

#include <string> // In order to work with strings -- needed to loop through the Mathematica output files

#include "exprtk.hpp" // Parser to evaluate string expressions as mathematical/arithmetic input

#include <math.h> // Simple Math Stuff
#include <gsl/gsl_math.h> // GSL math stuff
#include <complex> // Complex Numbers

/**********/
/* Parser */
/**********/

// Type definitions for the parser
typedef exprtk::symbol_table<double> symbol_table_t; // (%)
typedef exprtk::expression<double>     expression_t; // (%)
typedef exprtk::parser<double>             parser_t; // (%)

/* This struct is used to store certain information of the Mathematica files in
order to later evaluate them for different variables with the parser library. */
struct readInExpression
{
  double a,b; // (%)

  symbol_table_t symbol_table;

  // Instantiate expression
  expression_t expression;
};

/* Global variable where the read-in file/parser is stored. */
readInExpression file;
parser_t parser;

/*******************/
/* Custom function */
/*******************/

double my_function(double a, double b) {
  return a+b;
}

/***********************************/
/* Converting Mathematica Notation */
/***********************************/

/* Mathematica prints complex numbers as Complex(x,y), so we need a function to convert to C++ standard. */
std::complex<double> Complex(double a, double b) { // (%)
  std::complex<double> c(a,b);
  return c;
}

/************************************/
/* Processing the Mathematica Files */
/************************************/

double evaluateFileDoubleValuedInclude(double a, double b) {
    return
    #include "xyz.txt"
    ;
}

std::complex<double> evaluateFileComplexValuedInclude(double a, double b) {
    return
    #include "xyzC.txt"
    ;
}

void readInFromFile(readInExpression* f, parser_t* p) {
  std::string file = "xyz.txt"; // (%)

  std::ifstream ifs(file);
  std::string content( (std::istreambuf_iterator<char>(ifs) ),
                       (std::istreambuf_iterator<char>()    ) );

  // Register variables with the symbol_table
  f->symbol_table.add_variable("a",f->a);
  f->symbol_table.add_variable("b",f->b);

  // Add custom functions to the evaluation list (see definition above)
  f->symbol_table.add_function("my_function",my_function); // (%)
  // f->symbol_table.add_function("Complex",Complex); // (%)

  // Register symbol_table to instantiated expression
  f->expression.register_symbol_table(f->symbol_table);

  // Compile the expression with the instantiate parser
  p->compile(content,f->expression);
}

std::complex<double> evaluateFile(readInExpression* f, double a, double b) { // (%)
  // Set the values of the struct to the input values
  f->a = a;
  f->b = b;

  // Evaluate the result for the upper values
  return f->expression.value();
}

int main() {

    exprtk::symbol_table<std::complex<double> > st1; // Works
    exprtk::expression<std::complex<double> >     e1; // Works
    // exprtk::parser<std::complex<double> >             p1; // Throws an error

    double a = 2.0;
    double b = 3.0;

    std::cout << "Evaluating the text file containing only double-valued functions via the #include method: \n" << evaluateFileDoubleValuedInclude(a,b) << "\n \n";
    std::cout << "Evaluating the text file containing complex-valued functions via the #include method: \n" << evaluateFileComplexValuedInclude(a,b) << "\n \n";

    readInFromFile(&file,&parser);
    std::cout<< "Evaluating either the double-valued or the complex-valued file [see the necessary changes tagged with (%)]:\n" << evaluateFile(&file,a,b) << "\n";

    return 0;
}

xyz.txt

a + b * my_function(a,b)

xyzC.txt

2.0*Complex(a,b) + 3.0*a

To get the MWE to work, put the exprtk.hpp file in the same folder where you compile.

Note that the return type of the evaluateFile(...) function can be/is std::complex<double>, even though only double-valued types are returned. Lines tagged with // (%) are subject to change when trying out the complex-valued file xyzC.txt.

Instantiating exprtk::parser<std::complex<double> > throws (among others)

./exprtk.hpp:1587:10: error: no matching function for call to 'abs_impl'
                                                              exprtk_define_unary_function(abs  )

while all other needed types seem to not complain about the type std::complex<double>.

Answer 1

I actually know next to nothing about ExprTk (just what I just read in its documentation and a bit of its code -- EDIT: now somewhat more of its code), but it seems to me unlikely that you'll be able to accomplish what you want without doing some major surgery to the package.

The basis of its API, as you demonstrate in your question, are template objects specialised on a single datatype. The documentation says that that type "…can be any floating point type. This includes… any custom type conforming to an interface comptaible (sic) with the standard floating point type." Unfortunately, it doesn't clarify what they consider the interface of the standard floating point type to be; if it includes every standard library function which could take a floating point argument, it's a very big interface indeed. However, the distribution includes the adaptor used to create a compatible interface for the MPFR package which gives some kind of idea what is necessary.

But the issue here is that I suspect you don't want an evaluator which can only handle complex numbers. It seems to me that you want to be able to work with both real and complex numbers. For example, there are expressions which are unambiguous for real numbers and somewhat arbitrary for complex numbers; these include a < b and max(a, b), neither of which are implemented for complex types in C++. However, they are quite commonly used in real-valued expressions, so just eliminating them from the evaluation language would seem a bit arbitrary. [Note 1] ExprTK does assume that the numeric type is ordered for some purposes, including its incorrect (imho) delta equality operator, and those comparisons are responsible for a lot of the error messages which you are receiving.

The ExprTK header does correctly figure out that std::complex<double> is a complex type, and tags it as such. However, no implementation is provided for any standard function called on complex numbers, even though C++ includes implementations for many of them. This absence is basically what triggers the rest of the errors, including the one you mention for abs, which is an example of a math function which C++ does implement for complex numbers. [Note 2]

That doesn't mean that there is no solution; just that the solution is going to involve a certain amount of work. So a first step might be to fill in the implementations for complex types, as per the MPFR adaptor linked to above (although not everything goes through the _impl methods, as noted above with respect to ordered comparison operators).

One option would be to write your own "real or complex" datatype; in a simple implementation, you could use a discriminated union like:

template<typename T>
class RealOrComplex {
  public:
    RealOrComplex(T a = 0)
        :  is_complex_(false), value_.real(a) {}
    RealOrComplex(std::complex<T> a)
        : is_complex_(true), value_.cplx(a) {}

    // Operator implementations, omitted

  private:
    bool is_complex_;
    union {
      T real;
      std::complex<T> cmplx;
    } value_;
};

A possible simpler but more problematic approach would be to let a real number simply be a complex number whose imaginary part is 0 and then write shims for any missing standard library math functions. But you might well need to shim a large part of the Loop library as well.

So while all that is doable, actually doing it is, I'm afraid, too much work for an SO answer.

Since there is some indication in the ExprTK source that the author is aware of the existence of complex numbers, you might want to contact them directly to ask about the possibility of a future implementation.

---

Notes

1. It seems that Mathematica throws an error, if these operations are attempted on complex arguments. On the other hand, MatLab makes an arbitrary choice: ordered comparison looks only at the real part, but maximum is (inconsistently) handled by converting to polar coordinates and then comparing component-wise.

2. Forcing standard interfaces to be explicitly configured seems to me to be a curious implementation choice. Surely it would have been better to allow standard functions to be the default implementation, which would also avoid unnecessary shims like the explicit reimplementation of log1p and other standard math functions. (Perhaps these correspond to known inadequacies in certain math libraries, though.)

Answer 2

... can anyone here suggest a different math parser that is user friendly and can deal with complex numbers of the type std::complex(double), at the same time allowing to define custom functions that themselves return such complex values?

I came across only one math parser that handled complex numbers, Foreval: https://sourceforge.net/projects/foreval/

Implementation as dll library. You can connect real and complex variables of the "double" and "extended" type to Foreval in any form, passing the addresses of the variables. Has built-in standard functions with complex variables. Also, you can connect external functions with complex variables. But only the type of complex variables, passed and returned in functions is own, internal and will be different from std::complex(double). There are examples for GCC in the source.

There are two disadvantages: There is only a 32-bit version of Foreval.dll (connection is possible only for 32-bit program). And only for OS Windows.

But there are also advantages: Foreval.dll is a compiler, that generates machine code (fast calculations). There is a real type with floating point - "extended" ("double" is too), also for complex numbers.

Answer 3

The ExprTk library provides a simple type adaptor that wraps std::complex, which can be found here:

https://www.partow.net/programming/exprtk/index.html#downloads

It can be used as the basis of a more complete complex number type implementation that can be used with ExprTk. The following is a simple usage example:

#include "exprtk_complex_adaptor.hpp"
#include "exprtk.hpp"
 .
 .
 .

typedef exprtk::symbol_table<cmplx::complex_t> symbol_table_t;
typedef exprtk::expression<cmplx::complex_t>   expression_t;
typedef exprtk::parser<cmplx::complex_t>       parser_t;

const std::string expression = "(x + i) / sin(y + 2i)";

symbol_table_t symbol_table;

T i = T(0.0,1.0);
T x = T(1.1,0.0);
T y = T(2.2,0.0);

symbol_table.add_variable ("i" , i );
symbol_table.add_variable ("x" , x );
symbol_table.add_variable ("y" , y );

expression_t expression;
expression.register_symbol_table(symbol_table);

parser_t parser;
parser.compile(expression, expression);

const auto result = expression.value();