Mu Parser, Dynamically adding arbitrary user functions at run time; Wanting better solutions in C or C++

Answering this question in C will be more difficult than C++. Skip to my first answer post if you want to see an answer written in "C".

Questions about muParser dynamic run time functions have been asked several times, on several sites, with good reception. The given answers are incomplete, theoretical (whiteboard design) and sometimes wrong. Therefore, the answers posted are often grossly untested. I will offer pure C or C++ examples (as you prefer) to test with.

Here's the basic question with minimal C++ code:

C++ Dynamically Define Function

The existing answers, though, are theoretical and miss several important points about muparser.

muParser and gnuplot, the background of my problem:

Programs computing math expressions (eg: muparser is one example) can use expressions compiled in source code. Practical end-user math programs often need equations entered at run time. My favorite examples are gnuplot and octave. Gnuplot is not recompiled from C source each time a function is plotted. The user functions are entered while the program runs. for example:

f(x)=x + x**2 + x**3


Gnuplot will parse the example, and compute the function as needed.

I want to do the same thing with muParser. One practical use is to improve Gnuplot for scientific purposes by adding a second higher precision calculator to it capable of __float128.

(Edit) @Ethan Merritt took exception with my original comment about the difficulty of upgrading Gnuplot directly. After a chat, he said that he might improve the documentation on Gnuplot code layout; Despite our talk, Ethan believes __float128 log() library is not portable (eg: all bets are off), and won't ever be a Gnuplot patch. So, until those docs are written, I think muParser exploration is more likely to succeed.

MuParser is coded with float size independence in mind (true portable code). MuParser does work with __float128 values and functions using my headers. Replacing standard C/C++ math library headers (#includes) with a few macro definitions is all that's required to upgrade MuParser to use Quad precision math. I give a basic, tested C++ header at the end of this post to do exactly that. A second C header is required to make the muParser C API (and ABI) handle QuadMath.

Why 128 bits + stackexchange's answers aren't enough:

I found that none of the answers shown on stackExchange actually work correctly when recoded into C in muParser-2.2.5/samples/example2/example2.c (or used straight in the C++ version).

The following run time session shows what happens when extending muparser with one partial answer from stackExchange:

x=0
f(x)=x+1
f(2)
ANS=1
x
ANS=0


This is what gnuplot's calculator does (native unmodified gcc parser).

rabbitport /etc/portage # gnuplot

G N U P L O T
Version 5.2 patchlevel 2 (Gentoo revision r0) last modified 2017-11-15
Terminal type is now 'x11'
gnuplot> x=0
gnuplot> f(x)=x+1
gnuplot> print f(2)
3
gnuplot> print x
0
gnuplot>


Gnuplot will print that f(2) is 3, because local x is 2 and 2 plus 1 is 3.

Muparser will print that f(2) is 1, because global x is zero and 0 plus 1 is 1

So, the Stack Exchange answers don't really show how to implement a practical/general user function (like I need). The answers don't even warn the user that function parameters will likely conflict with global variables.

If I recode the answer from Stack Exchange to assign a new value to x, using a single parser for both function and function call; the result is that the global variable x also becomes 2. Obviously, a bad idea.

Two parsers are needed to evaluate both an expression invoking a user function and the user function itself, without conflict.

Why I'm still looking for an improved solution

I've compiled and tinkered with dozens of programs. But, generally they all share the same problems as gnuplot and octave; coding is based on double, and extensions based on MPFR.

GnuPlot, out of the box, uses doubles; The undeflow characteristics of the Intel processor log() function ruin several of the graphs with roundoff error due to underflow and insufficient precision.

Gnu Octave can do multi-precision math. Tests with Gnu octave using MPFR library took 1 hour to complete a simple log() task. The same task hard coded in C++ with gnu quadword library took 3 seconds. In other cases, gnu octave using MPFR locked up for a week and crashed. That makes examinig plots and doing "what if" experiments very hard.

EDIT: Log expressions, similar to this example, are what I'm evaluating.
print sprintf( "%25.20Lf", lgamma(x+100.)-lgamma(x+101.) ) ) # log gamma in Gnuplot.

Log calculations are computationally expensive, time-wise, compared to + - * /
When I need to caclualte large numbers of data points, equations like this one can slow GnuOctave down significantly.
Gnuplot is fast but will calculate the result with less than double accuracy.
-4.60517018598807 226226 ( a double, space added to show 15th digit. )
-4.60517018598812 ... ( mu-parser, or C with long-double ).


That's why I am after mu parser. It's able to do calculations significantly larger than double, but without the major slowdown caused by MPFR in octave. I often need around 20 digits of an answer.

What I offer to you, and what I hope to see on Stack Overflow:

There isn't enough space in this Stack Overflow to give both the C and C++ code that will allow people to test their answers. What I would like, however, are tested answers that aren't "theoretical."

Therefore I'm only posting the C++ header required to compile muparser's C++ interface with quadfloats.

If you need to modify the muParser library, itself, to produce your answer; I don't have a problem with that. Answers that improve the parser library, itself, are welcome.

Introduction to the C++ header that I provide for Quadmath:

The header extends C++ floating point standards based on the what is done for doubles. If you are used to coding with the C++ math library and doubles, then using my header should be simple. Just replace double with a float128 type. There is a macro USE_QUADMATH that needs to be defined for the header to include Quadmath functions. Otherwise the header will default back to the standard library. The macro exists to allow compiling the same program with or without quadmath extensions, and for compiling libraries.

I think many of you will find this header very useful in other C++ projects. Since GCC/FSF doesn't distribute a C++ header for quadmath, what I give is better than nothing (or Fortran, only). The GPL2.0 license is neeeded because the header is a derivative work, translating C to C++.

File cmathe follows.

#ifndef __CMATHE
#define __CMATHE

namespace mathspace=std ;
// The namespace of all cmath functions ought to be std,
// however, if your system's broken, uncomment the following define to kludge it back to global namespace.
// or just define mathspace as a compile time macro with no value.
// #define mathspace


#ifdef USE_QUADMATH
//#warning "USING experimental quadmath routines, they have >32 digits of precision."
#include <quadmath.h>
#include <math.h>
namespace std

// 128 bit floats, have at least 32 decimal digits of precision

static inline __float128 (cos)( __float128 x ) return cosq( x );
static inline __float128 (sin)( __float128 x ) return sinq( x );
static inline __float128 (tan)( __float128 x ) return tanq( x );
static inline __float128 (acos)( __float128 x ) return acosq( x );
static inline __float128 (asin)( __float128 x ) return asinq( x );
static inline __float128 (atan)( __float128 x ) return atanq( x );
static inline __float128 (atan2)( __float128 x, __float128 y ) return atan2q( x,y );

static inline __float128 (cosh)( __float128 x ) return coshq( x );
static inline __float128 (sinh)( __float128 x ) return sinhq( x );
static inline __float128 (tanh)( __float128 x ) return tanhq( x );
static inline __float128 (acosh)( __float128 x ) return acoshq( x );
static inline __float128 (asinh)( __float128 x ) return asinhq( x );
static inline __float128 (atanh)( __float128 x ) return atanhq( x );

static inline __float128 (exp)( __float128 x ) return expq( x );
static inline __float128 (frexp)( __float128 x, int *i ) return frexpq( x,i );
static inline __float128 (ldexp)( __float128 x, int i ) return ldexpq( x,i );
static inline __float128 (log)( __float128 x ) return logq( x );
static inline __float128 (log10)( __float128 x ) return log10q( x );
static inline __float128 (modf)( __float128 x, __float128 *y ) return modfq( x,y );
static inline __float128 (exp2)( __float128 x ) return expq( x * M_LN2q ); // FIXME: best way to implement?
static inline __float128 (expm1)( __float128 x ) return expm1q( x );
static inline int (ilogb)( __float128 x ) return ilogbq( x );
static inline __float128 (log1p)( __float128 x ) return log1pq( x );
static inline __float128 (log2)( __float128 x ) return log2q( x );

static inline __float128 (logb)( __float128 x ) return (__float128)ilogbq( x );

static inline __float128 (scalbn)( __float128 x, int n ) return scalbnq( x,n );
static inline __float128 (scalbln)( __float128 x, long int ln ) return scalblnq( x,ln );

static inline __float128 (pow)( __float128 x, __float128 y ) return powq( x,y );
static inline __float128 (sqrt)( __float128 x ) return sqrtq( x );
static inline __float128 (cbrt)( __float128 x ) return cbrtq( x );
static inline __float128 (hypot)( __float128 x,__float128 y ) return hypotq( x,y );

static inline __float128 (erf)( __float128 x ) return erfq( x );
static inline __float128 (erfc)( __float128 x ) return erfcq( x );
static inline __float128 (tgamma)( __float128 x ) return tgammaq( x );
static inline __float128 (lgamma)( __float128 x ) return lgammaq( x );

static inline __float128 (ceil)( __float128 x ) return ceilq( x );
static inline __float128 (floor)( __float128 x ) return floorq( x );
static inline __float128 (fmod)( __float128 x,__float128 y ) return fmodq( x,y );
static inline __float128 (trunc)( __float128 x ) return truncq( x );
static inline __float128 (round)( __float128 x ) return roundq( x );
static inline long int (lround)( __float128 x ) return lroundq( x );
static inline long long int (llround)( __float128 x ) return llroundq( x );
static inline __float128 (rint)( __float128 x ) return rintq( x );
static inline long int (lrint)( __float128 x ) return lrintq( x );
static inline __float128 (nearbyint)( __float128 x ) return nearbyintq( x );
static inline __float128 (remainder)( __float128 x, __float128 y ) return remainderq( x,y );
static inline __float128 (remquo)( __float128 x, __float128 y, int *i ) return remquoq( x,y,i );

static inline __float128 (copysign)( __float128 x,__float128 y ) return copysignq( x,y );
static inline __float128 (nan)( char* tagp ) return nanq( tagp );
static inline __float128 (nextafter)( __float128 x,__float128 y ) return nextafterq( x,y );


static inline __float128 (fdim)( __float128 x, __float128 y ) return fdimq( x,y );
static inline __float128 (fmax)( __float128 x, __float128 y ) return fmaxq( x,y );
static inline __float128 (fmin)( __float128 x, __float128 y ) return fminq( x,y );

static inline __float128 (fabs)( __float128 x ) return fabsq( x );

// FIXME: Why is this sometimes defined as a function already !?
/*
static inline __float128 (abs)( __float128 x ) return fabsq( x );
*/
static inline __float128 (fma)( __float128 x,__float128 y, __float128 z ) return fmaq( x,y,z );

static inline int (signbit)( __float128 x ) return signbitq( x );
static inline bool (isinf)( __float128 x ) return isinfq( x );
static inline bool (isnan)( __float128 x ) return isnanq( x );

// The code below this point is an optimization and portability problem.
// I am shooting for most portable, and reasonably fast by tendency to be used.

#if __GNUC_PREREQ(2,97)
static inline bool (isunordered)( __float128 x, __float128 y ) return __builtin_isunordered(x,y);
static bool (fpclassify)( __float128 x )
__builtin_fpclassify( FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x );

static inline bool (isfinite)( __float128 x )
return (bool)__builtin_fpclassify( 0,0,1,1,1, x );

#else
#error "Compile with GCC 2.97 or above... otherwise, you must debug the boilerplate functions in the header."
static inline int (fpclassify)( __float128 x ) // FIXME not debugged functions...
if ( isnanq(x) ) return FP_NAN ;
if ( isinf(x) ) return FP_INF ;
if ( x == 0 ) return FP_ZERO;
if ( x < ?FIXME_VALUE? ) return FP_SUBNORMAL; // FIXME: This is wrong, figure out a test for subnormal.
return FP_NORMAL;

static inline bool (isunordered)( __float128 x, __float128 y )
static inline bool (isfinite)( __float128 x )
#endif

static inline bool (isnormal)( __float128 x )
return (bool)(fpclassify(x) == FP_NORMAL);

static inline bool (isgreater)( __float128 x, __float128 y )
return ! isunordered(x,y) && ( x > y );

static inline bool (isgreaterequal)( __float128 x, __float128 y )
return ! isunordered(x,y) && ( x >= y );

static inline bool (isless)( __float128 x, __float128 y )
return ! isunordered(x,y) && ( x < y );

static inline bool (islessequal)( __float128 x, __float128 y )
return ! isunordered(x,y) && ( x <= y );

static inline bool (islessgreater)( __float128 x, __float128 y )
return ! isunordered(x,y) && ( x != y );


// These have prototype issues when it comes to roundoff.
// This could be replaced by a template function, most likeley if I understood it better.

static inline __float128 (nexttoward)( __float128 x,__float128 y )
return nextafter(x,y);

static inline long double (nexttoward)( long double x,__float128 y ) (y==x)) return (long double)y;
return nexttoward( x, ((y > (__float128)x)? HUGE_VALL:- HUGE_VALL) );


/*
// FIXME: Became ambiguous over a year, not sure why.
static inline double (nexttoward)( double x,__float128 y ) (y==x)) return (double)y;
return nexttoward( x, ((y > (__float128)x)? HUGE_VAL:- HUGE_VAL) );

static inline float (nexttoward)( float x,__float128 y )
*/

// namespace

#ifdef _USE_MATH_DEFINES
#define _USE_MATH_DEFINESq
#define M_E M_Eq
#define M_LOG2E M_LOG2Eq
#define M_LOG10E M_LOG10Eq
#define M_LN2 M_LN2q
#define M_LN10 M_LN10q
#define M_PI M_PIq
#define M_PI_2 M_PI_2q
#define M_PI_4 M_PI_4q
#define M_1_PI M_1_PIq
#define M_2_PI M_2_PIq
#define M_2_SQRTPI M_2_SQRTPIq
#define M_SQRT2 M_SQRT2q
#define M_SQRT1_2 M_SQRT1_2q
#endif

#undef _USE_MATH_DEFINES

#endif // USE QUADMATH

// ------------------------------------------------------------------------
// All changes to cmath have been created, so now include cmath like normal.
// This will undef a large number of troublesome math.h macros.
#include <cmath>


#ifdef USE_QUADMATH

// Restore macro if the math defines are desired...
#ifdef _USE_MATH_DEFINESq
#define _USE_MATH_DEFINES
#else
#undef M_Eq
#undef M_LOG2Eq
#undef M_LOG10Eq
#undef M_LN2q
#undef M_LN10q
#undef M_PIq
#undef M_PI_2q
#undef M_PI_4q
#undef M_1_PIq
#undef M_2_PIq
#undef M_2_SQRTPIq
#undef M_SQRT2q
#undef M_SQRT1_2q
#endif

// Get rid of any macros ... since this is C++

#undef cos
#undef sin
#undef tan
#undef acos
#undef asin
#undef atan
#undef atan2
#undef cosh
#undef sinh
#undef tanh
#undef acosh
#undef asinh
#undef atanh
#undef exp
#undef frexp
#undef ldexp
#undef log
#undef log10
#undef modf
#undef exp2
#undef expm1
#undef ilogb
#undef log1p
#undef log2
#undef logb
#undef scalbn
#undef scalbln
#undef pow
#undef sqrt
#undef cbrt
#undef hypot
#undef erf
#undef erfc
#undef tgamma
#undef lgamma
#undef ceil
#undef floor
#undef fmod
#undef trunc
#undef round
#undef lround
#undef llround
#undef rint
#undef lrint
#undef nearbyint
#undef remainder
#undef remquo
#undef copysign
#undef nan
#undef nextafter
#undef fdim
#undef fmax
#undef fmin
#undef fabs
#undef abs
#undef fma
#undef signbit
#undef insinf
#undef isnan
#undef isunordered
#undef fpclassify
#undef isfinite
#undef isnormal
#undef igreater
#undef isgreaterequal
#undef isless
#undef islessequal
#undef islessgreater
#undef nexttoward

// -------------------------------------------------------------------------
// Also define the limits here, although it's technically the wrong header.
#include <limits>

namespace std _GLIBCXX_VISIBILITY(default)

_GLIBCXX_BEGIN_NAMESPACE_VERSION
template<>
struct numeric_limits<__float128>

static _GLIBCXX_USE_CONSTEXPR bool is_specialized = true;

static _GLIBCXX_CONSTEXPR __float128
min() _GLIBCXX_USE_NOEXCEPT return FLT128_MIN;

static _GLIBCXX_CONSTEXPR long double
max() _GLIBCXX_USE_NOEXCEPT return FLT128_MAX;

#if __cplusplus >= 201103L
static constexpr __float128
lowest() noexcept return -FLT128_MAX;
#endif

static _GLIBCXX_USE_CONSTEXPR int digits = FLT128_DIG;
static _GLIBCXX_USE_CONSTEXPR int digits10 =33; // IEEE754R minimum.
#if __cplusplus >= 201103L
#ifndef __glibcxx_max_digits10
#define __glibcxx_max_digits10(T)
(2 + (T) * 643L / 2136)
#endif
static _GLIBCXX_USE_CONSTEXPR int max_digits10
= __glibcxx_max_digits10( FLT128_MANT_DIG );
#endif
static _GLIBCXX_USE_CONSTEXPR bool is_signed = true;
static _GLIBCXX_USE_CONSTEXPR bool is_integer = false;
static _GLIBCXX_USE_CONSTEXPR bool is_exact = false;
static _GLIBCXX_USE_CONSTEXPR int radix = 2; // IEEE754R -- is binary128.

static _GLIBCXX_CONSTEXPR long double
epsilon() _GLIBCXX_USE_NOEXCEPT return FLT128_EPSILON;

static _GLIBCXX_CONSTEXPR long double
round_error() _GLIBCXX_USE_NOEXCEPT return 0.5L;

static _GLIBCXX_USE_CONSTEXPR int min_exponent = FLT128_MIN_EXP;
static _GLIBCXX_USE_CONSTEXPR int min_exponent10 = FLT128_MIN_10_EXP;
static _GLIBCXX_USE_CONSTEXPR int max_exponent = FLT128_MAX_EXP;
static _GLIBCXX_USE_CONSTEXPR int max_exponent10 = FLT128_MAX_10_EXP;

static _GLIBCXX_USE_CONSTEXPR bool has_infinity = true;
static _GLIBCXX_USE_CONSTEXPR bool has_quiet_NaN = true;
static _GLIBCXX_USE_CONSTEXPR bool has_signaling_NaN = has_quiet_NaN;
static _GLIBCXX_USE_CONSTEXPR float_denorm_style has_denorm
= denorm_present;
static _GLIBCXX_USE_CONSTEXPR bool has_denorm_loss=true; // FIXME:: implementation ?

static _GLIBCXX_CONSTEXPR __float128
infinity() _GLIBCXX_USE_NOEXCEPT return (__float128)(numeric_limits<double>::infinity());

static _GLIBCXX_CONSTEXPR __float128
quiet_NaN() _GLIBCXX_USE_NOEXCEPT return (__float128)(numeric_limits<double>::quiet_NaN());

static _GLIBCXX_CONSTEXPR __float128
signaling_NaN() _GLIBCXX_USE_NOEXCEPT return (__float128)(numeric_limits<double>::signaling_NaN()) ;

static _GLIBCXX_CONSTEXPR __float128
denorm_min() _GLIBCXX_USE_NOEXCEPT return FLT128_DENORM_MIN;

static _GLIBCXX_USE_CONSTEXPR bool is_iec559
= has_infinity && has_quiet_NaN && has_denorm == denorm_present;
static _GLIBCXX_USE_CONSTEXPR bool is_bounded = true;
static _GLIBCXX_USE_CONSTEXPR bool is_modulo = false;

static _GLIBCXX_USE_CONSTEXPR bool traps = true; // FIXME:: implementation ?
static _GLIBCXX_USE_CONSTEXPR bool tinyness_before = true; // FIXME:: implementation ?
static _GLIBCXX_USE_CONSTEXPR float_round_style round_style =
round_to_nearest;
;

_GLIBCXX_END_NAMESPACE_VERSION
// namespace

// ----------------------------------------------------------------------------
// STREAM templates for doing stream based conversion of __float128
// Note: these really don't belong here, but this is a kludge to get
// the project working quickly.
#include <iostream>

template <class T> T& cmathe_insertion_tp( T& os, __float128 x )
int precision=os.precision();
typename T::ios_base::fmtflags ff = os.flags();
char fmt[10] = "%#.*Qg";
char buf[54];

if ( ff & os.fixed ) fmt[5]= ((ff & os.uppercase)?'F':'f');
else
if ( ff & os.scientific ) fmt[5]= ((ff & os.uppercase)?'E':'e');
else if (ff & os.uppercase) fmt[5]='G';

quadmath_snprintf(buf,54,fmt,precision,x);
os << buf;
return os;

// FIXME: The following routine doesn't have any error handling for
// when the input digits overflow any realistic expectation.
// It just truncates, and returns the number parsed so far.
// There's no elegant way to do error handling, because there could
// be 1000 digits entered before the exponent; and it's likely that
// truncation destroys the exponent portion of the parse.
// It's ugly to work out all the possible cases of digit lengths
// with and without decimal point, and tracking the exponent corrections.
// This parsing routine doesn't handle locales; either.
// Parsing problems include finge cases like '+.' and '-.e314'

#include <cctype> // for isdigit.
template <class T> T& cmathe_extraction_tp( T& is, __float128 &x ) ( p == buf ) )
is.setstate( T::failbit );
return is;


x=strtoflt128( buf, NULL );
return is;



// ----------------------------------------------------------------------------
#endif // ifdef USE_QUADMATH (2nd group)
#endif


Example use of header.
File test.cpp

// Using GCC, this compiles with:
// g++ test.cpp -o test.x -lquadword
#define USE_QUADMATH // can be a command line compiler define.
#include <quadmath.h>
#include <iostream>
#include "cmathe" // cmath, extended to __float128

using namespace std;

// These are best put into a user header for convenience, depending
// on what your system's extended float size is.
typedef __float128 floate_t; // extended float type, definied as generic float.
ostream& operator << (ostream& cout,__float128 x)
return cmathe_insertion_tp<ostream>(cout,x);


int main(void)
floate_t a,aa;
long double b;
double c;

a = lgamma( (floate_t)100. ) - lgamma( (floate_t)101. );
aa= lgamma(100.) - lgamma( 101.);
b = lgamma( (long double)100. ) - lgamma( (long double)101. );
c = lgamma( (double)100. ) - lgamma( (double)101. );

cout.setf(ios::fixed); cout.precision(20);
cout << c << " double" << endl
<< b << " long double" << endl
<< a << " QUADWORD" << endl
<< aa << " promotion test " << endl;
return 0;

// Run on Intel, 4960X processor.
-4.60517018598807226226 double
-4.60517018598809138585 long double
-4.60517018598809136804 QUADWORD
-4.60517018598807226226 promotion test


@halfer

2 Answers
2

I have nothing to say about MuParser, but as a gnuplot developer I can say with some confidence that you have misjudged how difficult/easy it would be to change gnuplot evaluation to use quad precision rather than double precision. All reference to floating point values by the evaluation stack is based on

struct cmplx
double real, imag;
;


That is to say, the only change required to the gnuplot evaluation stack would be redefining struct cmplx to explicitly use 128bit precision. The problem is that to my knowledge there is no uniform standard for doing this. It is compiler and architecture dependent.

struct cmplx

It may be instructive to look at the small number of changes needed for the recent switch from 32bit to 64bit integer arithmetic in the gnuplot evaluation stack. I don't see why increasing the floating point precision would be any harder, except for the lack of a uniform standard to code against.

Here's my first answer, and its written purely in C. There are deficiencies that can be improved, such as error checking. A C++ implementation will be upcoming, later this week, on gitlab.

My approach uses the muparser add variable factory function to selectively link user function variables to the command line parser's global variable space. Each new user function is stored in a new muparser instance with all parameters being kept local, except the function name which is global. Separation of parser variable spaces prevents global and parameter names clobbering each other. Thanks goes to @n.m. who's chat comments helped me see the possibility of this answer. ;)

Deficiencies of my answer:

1) My answer relies on a sequential allocation order of parameters (local variables) so that the "C" function call "stack" parameters can be block copies into local user function variables. This technique would not work if local parameter memory was allocated in random locations with malloc, calloc, etc. The same kind of problem will arise in C++.

2) The Parameters, being local, can be corrupted by a recursive function call overwriting the parameters. In order to undo corruption, a second copy of call parameters needs to be done immediately after any user function is evaluated; since this takes extra time, I have not implemented a second copy. Copying the stack once, already slows my answer down. Therefore, recursive function calls are not supported.

3) In order to keep user function operation stable, error checking of parameter size and function trampoline number is done every time a user function is called. No optimizations are attempted. A better answer would include optimization that removed error checking code when safe. eg: implement the user call function callback with different levels of security and link to the appropriate callback depending on whether the function has been executed before or not.

4) MuParser variables are capable of storing muParser handles. Therefore, all function calls can be managed by the parser's variables, directly. However, I wasn't sure how to do this without allowing users to tamper with the function addresses, and risk "execution" arbitrary memory locations (security risk). I chose to implement an un-necessary static global function "trampoline" that can't be tampered with by the user. This has the undesirable side effect of introducing thread safety unknowns, because of static variables. I don't know if my answer is thread safe, or if it will work with MP (multi-processor lib). I don't know how to work around or test this issue, yet.

5) Function parameter allocation necessarily has memory leaks when a function is re-defined. eg: because of issue #1, blocks of memory are allocated with random length and can't be reused, safely. I don't know how to work around this yet.

6) Traditional function definitions, like gnuplot's fn(a,b)=a*b , were not easy for me to implement. To simplify the definition code, I implemented function definition of 'foo' with the syntax def(foo , ...)=function__body . A better solution would match Gnuplot's style of definition. Unfortunately, f()= style requires catching the error when an open '(' follows a variable 'f' that is undefined. Restarting the parser after an error, mid parse, is not documented. defining "(" as a unary operator is also not documented; and my attempt failed.
Therefore, to prevent writing hairy code, I used an easier to parse syntax.

7) MuParser documentation does not show many examples of manipulating multi-element variables (lists). Functions can return tuples, and (maybe?) accept them, but I've made no attempt to handle multi-value lists other than to print them. Lists are an asset in Gnuplot that should be considered/tested out. (A TODO).

Bonuses of this technique:

1) All user functions can be called by all other user functions.

2) User defined functions are not limited to linear functions ( as in the stackexchange example ), but can be truly arbitrary. Built in functions like sin(), powers, etc. can all be called from a user function as well.

3) Since each parser instance (AKA user function) compiles to bytecode separately, each user function will save it's optimized byte code after being called once. The speed enhancements of reusing bytecode is maintained. Using a single parser instance can not keep more than one expression in byte-code, but multiple parser instances can. Only redefinition of a function will slow the calculation down during a re-parse.

The following code is most of the differences made to file: muparser-2.2.5/samples/example2/example2.c . This code snippet shows how I implement the linking of user function variables to global function space, and creation of "local" function parameters so that arbitrary user functions can be made.

//---------------------------------------------------------------------------
// User function implement ....

static muParserHandle_t hParser; // make command line parser, global.
static muParserHandle_t hUser[MAX_UFUNS + 1]; // number of user functions
static muInt_t szUser[MAX_UFUNS+1]; // number of parameters in function
static muFloat_t* pFUser[MAX_UFUNS+1]; // Starting parameter* of funciton

static muInt_t userNew=MAX_UFUNS; // next traampoline to allocate
static unsigned int ufi=0; // user function index number passback. (messy!)

/* Call user functions, with limits checking and error handling for parameters.
The extra error checking, naturally, slows down the function call.
Copying the ?stack? pushed parameters also slows it down.
Optimizing these two issues out, however, is not trivial.
Note:
muparser forces multi-functions, to have const on the muFloat_t params...
This const suggests that the parameters are not truly on a stack.
Swapping of stack parameters with local varaibles, is one way to
implement recursive safe user functions..
That's not possible with "const" applied to the parameter list.
The reason muparser has this restriction, requires more research.
*/
muFloat_t UFuns( const muFloat_t *pS0, muInt_t n )

muFloat_t* LinkGlobal(const muChar_t* a_szName, void *_hFunction )

muParserHandle_t *hFunction = _hFunction;
const muChar_t *pName;
muFloat_t *pVal;

mupSetExpr( hParser, a_szName ); (void)mupEval( hParser );
if (mupError(hParser)) return NULL; // Something went wrong.
mupGetExprVar( hParser, 0, &pName, &pVal );
if (hFunction)
return pVal;



The following snippet of code shows how the user string input is parsed to detect the definition of a function, and how the parser is exploited to link a user function by redirecting auto-variable (AKA parameter) creation to a linker function.

if (!mystrncmp(szLine,_T("def("),4 )) // define a function.
muParserHandle_t hFunction = mupCreate(muBASETYPE_FLOAT);
muChar_t *p=szLine+4;
int nPar=0;
while( *p && *p != _T(')') ) ++p; if (*p) *p++=0;
if (*p && *p++ != _T('='))
myprintf("syntax error, character is not '=' after def() n");
mupRelease(hFunction);
continue;

mupSetVarFactory( hFunction, LinkGlobal, &hFunction );
mupSetErrorHandler(hFunction,OnError);
mupSetExpr(hFunction,szLine+4); (void)mupEvalMulti(hFunction,&nPar);
if (mupError(hFunction)) continue;
szUser[ufi]=nPar; // non-threadsafe, global solution.
myprintf(_T("number of parameters is %dn"),nPar);

// Function is set up, along with locals. Parse f() body.
if (!*p) p=_T("NAN"); myprintf(_T("Empty function=NAN!n"));
myprintf(_T("function body=%sn"),p);
mupSetExpr(hFunction,p);
mupSetVarFactory(hFunction,LinkGlobal,NULL);

continue;



I am preparing a gitlab repository, that includes the QUADLIB header file for C and C++. This will allow you to build, test, and modify this solution with the complete QUADWORD upgrade I need for gnuplot. I am not posting the diffs needed to compile muparser library as QUADWORD for C in this post. They will exist, however, in the gitlab repository so that muparser lib compiles 128 bit by default.

I will try to have gitlab up by Wed Aug22,2018; and will edit this post to keep you updated.

A test run of the parser, with all debugging output included, follows:

Executing "/tmp/muparser/muparser-2.2.5/samples/example2/example2.x" (argc=1)
__________
_____ __ ________ _____ _______ ______ ____ _______
/ | | | ___/__ _ __ / ___/_/ __ \_ __
| Y Y | | /| | / __ _| | /___ ___/ | | /
|__|_| /|____/ |____| (____ /|__| /____ > ___ >|__|
/ / / /
Version 2.2.5 (20150427; GC; 64BIT; RELEASE; ASCII) (DLL)
Sample build with ASCII support
(C) 2015 Ingo Berg
---------------------------------------
Commands:
list var - list parser variables
list exprvar - list expression variables
list const - list all numeric parser constants
locale de - switch to german locale
locale en - switch to english locale
locale reset - reset locale
test bulk - test bulk mode
quit - exits the parser

---------------------------------------
Constants:
"_e" 2.718281828459045235360287
"_pi" 3.141592653589793238462643
---------------------------------------
Please enter an expression:
def(fn,a,b)=a*b
Generating new variable "fn" (slots left: 100; context pointer: 0x0)
Generating new variable "(null)" (slots left: 99; context pointer: 0x0)
Generating new variable "a" (slots left: 98; context pointer: 0x0)
Generating new variable "b" (slots left: 97; context pointer: 0x0)
number of parameters is 3
function body=a*b
a,b
ANS=1.000000,2.000000,
U(fn,7,8)
Parameters=7.000000,8.000000, --> a=7.000000,b=8.000000,fn=100.000000,
ANS=56.000000,
a,b
ANS=1.000000,2.000000,
def(fn,a)=sin(a)
Generating new variable "(null)" (slots left: 96; context pointer: 0x0)
Generating new variable "a" (slots left: 95; context pointer: 0x0)
number of parameters is 2
function body=sin(a)
U(fn,1.57)
Parameters=1.570000, --> a=1.570000,fn=100.000000,
ANS=1.000000,
U(fn,0.5)
Parameters=0.500000, --> a=0.500000,fn=100.000000,
ANS=0.479426,


As you can see, it works. :)

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.