cfun.h

Go to the documentation of this file.

//                  CVM Class Library
//                  http://cvmlib.com
//
//          Copyright Sergei Nikolaev 1992-2014
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#ifndef _CFUN_H
#define _CFUN_H

#include "cvm.h"        // http://cvmlib.com

#include <ctype.h>
#include <vector>

#define CFUN_K_FUN                  1
#define CFUN_K_VAR                  2
#define CFUN_K_CONST                3
#define CFUN_K_PLUS                 10
#define CFUN_K_MINUS                11
#define CFUN_K_MULT                 12
#define CFUN_K_DIV                  13
#define CFUN_K_POWER                14
#define CFUN_K_EXP                  16
#define CFUN_K_SQRT                 17
#define CFUN_K_LOG                  18
#define CFUN_K_LOG10                19
#define CFUN_K_SIN                  20
#define CFUN_K_COS                  21
#define CFUN_K_TAN                  22
#define CFUN_K_ARCSIN               23
#define CFUN_K_ARCCOS               24
#define CFUN_K_ARCTAN               25
#define CFUN_K_SINH                 26
#define CFUN_K_COSH                 27
#define CFUN_K_TANH                 28
#define CFUN_K_UMINUS               29
#define CFUN_K_SIGN                 30
#define CFUN_K_ABS                  31
#define CFUN_K_DELTA                32
#define CFUN_K_INFINITY             33
#define CFUN_K_MNINFINITY           34
#define CFUN_K_INTSIN               35
#define CFUN_K_INTCOS               36
#define CFUN_K_SAT                  37
#define CFUN_K_DSOLVE               38
#define CFUN_K_IIF                  39
#define CFUN_O_BRACE                '{'
#define CFUN_C_BRACE                '}'
#define CFUN_O_BRACKET              '['
#define CFUN_C_BRACKET              ']'
#define CFUN_O_PARENTH              '('
#define CFUN_C_PARENTH              ')'
#define CFUN_O_SPARENTH             "("
#define CFUN_C_SPARENTH             ")"
#define CFUN_COMMA                  ','
#define CFUN_PLUS                   '+'
#define CFUN_MINUS                  '-'
#define CFUN_POWER                  '^'
#define CFUN_MULT                   '*'
#define CFUN_DIV                    '/'
#define CFUN_POINT                  '.'
#define CFUN_SIGNS                  "()+-^*/"
#define CFUN_INF                    "INF"

#define CFUN_EXP                    "exp"
#define CFUN_SQRT                   "sqrt"
#define CFUN_LOG                    "log"
#define CFUN_LOG10                  "log10"
#define CFUN_SIN                    "sin"
#define CFUN_COS                    "cos"
#define CFUN_TAN                    "tan"
#define CFUN_ASIN                   "asin"
#define CFUN_ACOS                   "acos"
#define CFUN_ATAN                   "atan"
#define CFUN_SINH                   "sinh"
#define CFUN_COSH                   "cosh"
#define CFUN_TANH                   "tanh"
#define CFUN_SI                     "sinint"
#define CFUN_CI                     "cosint"
#define CFUN_SMINUS                 "-"
#define CFUN_SPLUS                  "+"
#define CFUN_SMULT                  "*"
#define CFUN_SDIV                   "/"
#define CFUN_SPOWER                 "^"
#define CFUN_SCOMMA                 ","
#define CFUN_SSPACE                 " "
#define CFUN_SIGN                   "sign"
#define CFUN_ABS                    "abs"
#define CFUN_IIF                    "iif"
#define CFUN_DELTA                  "delta"
#define CFUN_POWERS                 "power"
#define CFUN_SAT                    "sat"
#define CFUN_I                      "i"
#define CFUN_SPACES                 " \t\r\n"
#define CFUN_SINF                   CFUN_O_SPARENTH CFUN_INF CFUN_C_SPARENTH
#define CFUN_SMINF                  CFUN_O_SPARENTH CFUN_SMINUS CFUN_INF CFUN_C_SPARENTH

// find* functions failure code
#define CFUN_NOT_FOUND              size_t(~0)

#if defined (CVM_FLOAT)
#   define             CFUN_M_ZERO                         0.F
#   define             CFUN_M_ONE                          1.F
#   define             CFUN_M_TWO                          2.F
#   define             CFUN_M_MONE                         -1.F
#   define             CFUN_M_HALF                         0.5F
#   define             CFUN_M_E                            2.7182818284590452353602874713527F
#   define             CFUN_M_LN2                          0.69314718055994530941723212145818F
#   define             CFUN_M_PI                           3.1415926535897932384626433832795F
#   define             CFUN_M_PI_2                         1.5707963267948966192313216916398F
#   define             CFUN_M_LN_10                        2.3025850929940456840179914546844F
#   define             CFUN_M_GAMMA                        0.57721566490153286060651209008240243F
#else
#   define             CFUN_M_ZERO                         0.L
#   define             CFUN_M_ONE                          1.L
#   define             CFUN_M_TWO                          2.L
#   define             CFUN_M_MONE                         -1.L
#   define             CFUN_M_HALF                         0.5L
#   define             CFUN_M_E                            2.7182818284590452353602874713527L
#   define             CFUN_M_LN2                          0.69314718055994530941723212145818L
#   define             CFUN_M_PI                           3.1415926535897932384626433832795L
#   define             CFUN_M_PI_2                         1.5707963267948966192313216916398L
#   define             CFUN_M_LN_10                        2.3025850929940456840179914546844L
#   define             CFUN_M_GAMMA                        0.57721566490153286060651209008240243L
#endif

#define                CFUN_OK                            CVM_OK
#define                CFUN_PARSEERROR                    CVM_THE_LAST_ERROR_CODE + 1  
#define                CFUN_DOMAINERROR                   CVM_THE_LAST_ERROR_CODE + 2  
#define                CFUN_DOMAINERROR_C                 CVM_THE_LAST_ERROR_CODE + 3  
#define                CFUN_CONVERGENCEERROR              CVM_THE_LAST_ERROR_CODE + 4  
#define                CFUN_CONVERGENCEERROR_C            CVM_THE_LAST_ERROR_CODE + 5  
#define                CFUN_SUBSTPARAMETERERROR           CVM_THE_LAST_ERROR_CODE + 6  
#define                CFUN_VARSDONTMATCH                 CVM_THE_LAST_ERROR_CODE + 7  
#define                CFUN_NULLPOINTERERROR              CVM_THE_LAST_ERROR_CODE + 8  
#define                CFUN_PARAMETER_RECURSION           CVM_THE_LAST_ERROR_CODE + 9  

#define                CFUN_SIMPS_STACK_DEPTH             32
#define                CFUN_MAX_EI_ITERATIONS             1000


CVM_NAMESPACE_BEG



// specializations for some elementary functions
template <typename T>
class ElementaryFunctions 
{
    typedef std::complex<T> TC;

public:
    static T asin(const T& v) throw(cvmexception)
    {
        if (v < CFUN_M_MONE || v > CFUN_M_ONE)
            throw cvmexception(CFUN_DOMAINERROR, "asin", v);
        return ::asin(v);
    }
    static T acos(const T& v) throw(cvmexception)
    {
        if (v < CFUN_M_MONE || v > CFUN_M_ONE)
            throw cvmexception(CFUN_DOMAINERROR, "acos", v);
        return ::acos(v);
    }
    static T tan(const T& v)
    {
        return ::tan(v);
    }
    static T atan(const T& v)
    {
        return ::atan(v);
    }
    static T sinh(const T& v)
    {
        return ::sinh(v);
    }
    static T cosh(const T& v)
    {
        return ::cosh(v);
    }
    static T tanh(const T& v)
    {
        return ::tanh(v);
    }
    static T sqrt(const T& v) throw(cvmexception)
    {
        if (v < CFUN_M_ZERO)
            throw cvmexception(CFUN_DOMAINERROR, "sqrt", v);
        return ::sqrt(v);
    }
    static T log(const T& v) throw(cvmexception)
    {
        if (v < cvm::basic_cvmMachMin<T>())
            throw cvmexception(CFUN_DOMAINERROR, "log", v);
        return ::log(v);
    }
    static T log10(const T& v) throw(cvmexception)
    {
        if (v < cvm::basic_cvmMachMin<T>())
            throw cvmexception(CFUN_DOMAINERROR, "log", v);
        return ::log10(v);
    }
    static T pow(const T& vBase, const T& vPower)
    {
        return ::pow(vBase, vPower);
    }

    // Sine integral
    static T sinint(const T& x, const T& eps) throw(cvmexception)
    {
        T si = CFUN_M_ZERO;
        T xa = cvm::_abs(x);

        if (xa > CFUN_M_TWO) {
            si = CFUN_M_PI_2 + e1(xa, eps).imag();
        }
        else {
            si = xa;
            if (xa > eps) {
                T t;
                T f = xa;
                int sign = -1;
                for (int i = 3; i < CFUN_MAX_EI_ITERATIONS; i += 2) {
                    f *= xa * xa / T((i - 1) * i);
                    t = f / T(i);
                    si += T(sign) * t;
                    sign = -sign;
                    if (t / cvm::_abs(si) < eps)
                        break;
                    if (i >= CFUN_MAX_EI_ITERATIONS - 2) {
                        throw cvmexception(CFUN_CONVERGENCEERROR, "sinint", x);
                    }
                }
            }
        }
        if (x < CFUN_M_ZERO) {
            si = -si;
        }
        return si;
    }

    // Cosine integral
    static T cosint(const T& x, const T& eps) throw(cvmexception)
    {
        if (x <= CFUN_M_ZERO) {
            throw cvmexception(CFUN_DOMAINERROR, "cosint", x);
        }

        T ci = CFUN_M_ZERO;
        if (x > CFUN_M_TWO) {
            ci = -e1(x, eps).real();
        }
        else {
            T t;
            T f = CFUN_M_ONE;
            int sign = -1;
            for (int i = 2; i < CFUN_MAX_EI_ITERATIONS; i += 2) {
                f *= x * x / T((i - 1) * i);
                t = f / T(i);
                ci += T(sign) * t;
                sign = -sign;
                if (t / cvm::_abs(ci) < eps)
                    break;
                if (i >= CFUN_MAX_EI_ITERATIONS - 2) {
                    throw cvmexception(CFUN_CONVERGENCEERROR, "cosint", x);
                }
            }
            ci += log(x) + CFUN_M_GAMMA;
        }
        return ci;
    }

    // Exponential integral routine
    static TC e1(const T& x, const T& eps) throw(cvmexception)
    {
        static const TC c1(CFUN_M_ONE, CFUN_M_ZERO);
        static const TC c2(CFUN_M_TWO, CFUN_M_ZERO);
        TC a;
        TC b(CFUN_M_ONE, x);
        TC c(CFUN_M_ONE / (eps * eps), CFUN_M_ZERO);
        TC d(c1 / b);
        TC e(d);

        for (int i = 1; i < CFUN_MAX_EI_ITERATIONS; ++i) {
            a = TC(T(-i * i), CFUN_M_ZERO);
            b += c2;
            d = c1 / (a * d + b);
            c = b + a / c;
            a = c * d;
            e *= a;
            if (cvm::_abs(CFUN_M_ONE - a.real()) + cvm::_abs(a.imag()) < eps)
                break;
            if (i >= CFUN_MAX_EI_ITERATIONS - 1) {
                throw cvmexception(CFUN_CONVERGENCEERROR, "ElementaryFunctions<T>::e1", x);
            }
        }
        return TC(cos(x), -sin(x)) * e;
    }
};

template <typename T>
class ElementaryFunctions<std::complex<T> >
{
    typedef std::complex<T> TC;
public:
    static TC asin(const TC& v)
    {
        static const TC c1(CFUN_M_ONE, CFUN_M_ZERO);
        static const TC cmi(CFUN_M_ZERO, CFUN_M_MONE);
        return cmi * std::log(std::sqrt(c1 - v * v) - cmi * v);
    }
    static TC acos(const TC& v)
    {
        static const TC c1(CFUN_M_ONE, CFUN_M_ZERO);
        static const TC cmi(CFUN_M_ZERO, CFUN_M_MONE);
        return cmi * std::log(v - cmi * std::sqrt(c1 - v * v));
    }
    static TC tan(const TC& v)
    {
        const T rm = v.real() * CFUN_M_TWO;
        const T im = v.imag() * CFUN_M_TWO;
        if (cvm::_abs(im) >= basic_cvmLogMachMax<T>()) {
            return TC(CFUN_M_ZERO, im > CFUN_M_ZERO ? CFUN_M_ONE : CFUN_M_MONE);
        }
        else {
            const T d = ::cos(rm) + ::cosh(im);
            return TC(::sin(rm) / d, ::sinh(im) / d);
        }
    }
    static TC atan(const TC& v)
    {
        static const TC ci(CFUN_M_ZERO, CFUN_M_ONE);
        static const TC ci2(CFUN_M_ZERO, CFUN_M_HALF);
        return ci2 * std::log((ci + v) / (ci - v));
    }
    static TC sinh(const TC& v)
    {
        return std::sinh(v);
    }
    static TC cosh(const TC& v)
    {
        return std::cosh(v);
    }
    static TC tanh(const TC& v)
    {
        return std::tanh(v);
    }
    static TC sqrt(const TC& v)
    {
        return std::sqrt(v);
    }
    static TC log(const TC& v)
    {
        return std::log(v);
    }
    static TC log10(const TC& v)
    {
        return std::log10(v);
    }
    static TC pow(const TC& vBase, const TC& vPower)
    {
        return std::pow(vBase, vPower);
    }

    // Sine integral
    static TC sinint(const TC& x, const T& eps) throw(cvmexception)
    {
        TC si(x);
        if (cvm::_abs(x) > eps) {
            TC t;
            TC f(x);
            int sign = -1;
            for (int i = 3; i < CFUN_MAX_EI_ITERATIONS; i += 2) {
                f *= x * x / T((i - 1) * i);
                t = f / T(i);
                si += t * T(sign);
                sign = -sign;
                if (cvm::_abs(t) / cvm::_abs(si) < eps)
                    break;
                if (i >= CFUN_MAX_EI_ITERATIONS - 2) {
                    throw cvmexception(CFUN_CONVERGENCEERROR_C, "sinint", x.real(), x.imag());
                }
            }
        }
        return si;
    }

    // Cosine integral
    static TC cosint(const TC& x, const T& eps) throw(cvmexception)
    {
        static const TC gamma(CFUN_M_GAMMA, CFUN_M_ZERO);
        if (cvm::_abs(x) < eps) {
            throw cvmexception(CFUN_DOMAINERROR_C, "cosint", x.real(), x.imag());
        }

        TC ci(CFUN_M_ZERO);
        TC t;
        TC f(CFUN_M_ONE);
        int sign = -1;
        for (int i = 2; i < CFUN_MAX_EI_ITERATIONS; i += 2) {
            f *= x * x / T((i - 1) * i);
            t = f / T(i);
            ci += t * T(sign);
            sign = -sign;
            if (cvm::_abs(t) / cvm::_abs(ci) < eps)
                break;
            if (i >= CFUN_MAX_EI_ITERATIONS - 2) {
                throw cvmexception(CFUN_CONVERGENCEERROR_C, "cosint", x.real(), x.imag());
            }
        }
        return ci + log(x) + gamma;
    }
};


// cross-type stuff
template <typename T1, typename T2>
class Comparator
{
public:
    static bool lt(const T1& l, const T2& r);
    static bool le(const T1& l, const T2& r);
    static bool gt(const T1& l, const T2& r);
    static bool ge(const T1& l, const T2& r);
    static bool eq(const T1& l, const T2& r);
};


template <typename T>
class shared_ptr_pair
{
    std::shared_ptr<T> mp1;  
    std::shared_ptr<T> mp2;  
public:

    shared_ptr_pair(T* p1, T* p2)
      : mp1(p1),
        mp2(p2)
    {
    }

    shared_ptr_pair(const std::shared_ptr<T>& p1, const std::shared_ptr<T>& p2)
      : mp1(p1),
        mp2(p2)
    {
    }

    shared_ptr_pair(std::shared_ptr<T>&& p1, std::shared_ptr<T>&& p2)
      : mp1(std::move(p1)),
        mp2(std::move(p2))
    {
    }

    std::shared_ptr<T> operator[](size_t i) throw(cvmexception)
    {
        if (i == 0) {
            return mp1;
        }
        else if (i == 1) {
            return mp2;
        }
        throw cvmexception(CVM_INDEX_GT, i, 1);
    }
};


typedef std::vector<std::string> string_array; 


inline size_t __not(size_t i)
{
    return i ? 0 : 1;
}

inline void __clean_array(string_array& a)
{
    if (a.size() > 0) {
        a.erase(a.begin(), a.end());
    }
}

template <typename T>
inline bool __parse_num(const std::string& s, T& result)
{
    bool ret = true;
    char* stop;
    result = T(::strtod(s.c_str(), &stop));
    if ((result == CFUN_M_ZERO && s[0] != '0') || result == HUGE_VAL || result == -HUGE_VAL || *stop != '\0') {
        ret = false;
    }
    return ret;
}

// 7.0 FIX: these big things must be inlined too!
// MS doesn't like to do memory allocation and deallocation in different versions of RTL.
// Therefore, to make it possible to mix MSVS 2012 and 2010 generated code we have to inline this stuff.

inline std::string::const_iterator __trim_beg(const std::string& s, size_t nEdge = 0)
{
    static const char szSpaces[] = CFUN_SPACES;
    std::string::const_iterator it = s.end();
    size_t nLen = s.length();

    if (nLen > nEdge * 2) {       // both sides
        nLen = s.find_first_not_of(szSpaces, nEdge);

        if (nLen != CFUN_NOT_FOUND) {
            it = s.begin() + nLen;
        }
    }
    return it;
}

inline std::string::const_iterator __trim_end(const std::string& s, size_t nEdge = 0)
{
    static const char szSpaces[] = CFUN_SPACES;
    std::string::const_iterator it = s.end();
    size_t nLen = s.length();

    if (nLen > nEdge * 2) {       // both sides
        nLen = s.find_last_not_of(szSpaces, nLen - nEdge - 1);

        if (nLen != CFUN_NOT_FOUND) {
            it = s.begin() + (nLen + 1);
        }
    }
    return it;
}

// returns false if empty
inline bool __strip_parenth(const std::string& s, size_t& nFirst, size_t& nLast)
{
    static const char szSpaces[] = CFUN_SPACES;
    static const std::string sSigns(CFUN_SIGNS);
    bool bRet = true;
    int nLevel;
    char c;
    const size_t nLastOld = nLast;
    if (s.length() <= nLast || (s[nLast] && sSigns.find(s[nLast]) != CFUN_NOT_FOUND)) {
        --nLast;
    }

    for (;;) {
        nFirst = s.find_first_not_of(szSpaces, nFirst);
        nLast = s.find_last_not_of(szSpaces, nLast);
        if (nFirst == CFUN_NOT_FOUND || nLast == CFUN_NOT_FOUND) {
            bRet = false;
            break;
        }
        if (s[nFirst] == CFUN_O_PARENTH && s[nLast] == CFUN_C_PARENTH) {    // if (xxxxx)
            nLevel = 1;
            for (size_t i = nFirst + 1; i < nLast; i++) {                   // but like (x) + (y) or (x*(2+y)-3)
                c = s[i];
                if (c == CFUN_C_PARENTH) {
                    nLevel--;
                    if (!nLevel)
                        break;
                }
                if (c == CFUN_O_PARENTH) {
                    nLevel++;
                }
            }
            if (nLevel) {
                nFirst++;
                nLast--;
                continue;
            }
            break;
        }
        else {
            break;
        }
    }
    if (nLastOld > nLast) {
        ++nLast;             // -> '\0'
    }
    return bRet;
}

inline size_t __tassign(std::string& sDest, const std::string& sSrc, size_t nFirst, size_t nLast)
{
    size_t nRet = 0;
    if (__strip_parenth(sSrc, nFirst, nLast)) {
        if (&sDest == &sSrc) {
            if (nLast < sDest.length()) {
                sDest.erase(sDest.begin() + nLast, sDest.end());
            }
            if (nFirst) {
                sDest.erase(sDest.begin(), sDest.begin() + nFirst);
            }
        }
        else {
            sDest.assign(sSrc.begin() + nFirst, sSrc.begin() + nLast);
        }
        nRet = nLast - nFirst;
    }
    else {
        sDest.erase();
    }
    return nRet;
}

// "(x, 2)" -> "x" and "2"
// 0 - empty string
// 1 - there was a comma
// 2 - there wasn't a comma
inline int __separate(const std::string& sBody, size_t nFirst, size_t nLast,
                      std::string& sLeft, std::string& sRight)
{
    if (__strip_parenth(sBody, nFirst, nLast)) {
        size_t nOPos = sBody.find(CFUN_O_PARENTH, nFirst);
        size_t nCPos = sBody.find(CFUN_C_PARENTH, nFirst);

        size_t i = nFirst;
        for (;;) {
            i = sBody.find(CFUN_COMMA, i);
            if (i == CFUN_NOT_FOUND) {
                __tassign(sLeft, sBody, nFirst, nLast);
                sRight.erase();
                return 2;
            }
            if (nOPos == CFUN_NOT_FOUND ||
                nCPos == CFUN_NOT_FOUND ||
                (nCPos > nOPos && (i < nOPos || i > nCPos))) {
                __tassign(sLeft, sBody, nFirst, i);
                __tassign(sRight, sBody, i + 1, nLast);
                return 1;
            }
            i++;
        }
    }
    sLeft.erase();
    sRight.erase();
    return 0;
}

inline void __push_back_unique(string_array& a, const std::string& v)
{
    string_array::const_iterator it = a.begin();
    string_array::const_iterator end = a.end();
    while (it != end) {
        if ((*it) == v)
            return;
        ++it;
    }
    a.push_back(v);
}

inline bool __arrays_equal(const string_array& a1, const string_array& a2)
{
    bool ret = a1.size() == a2.size();
    if (ret) {
        string_array::const_iterator i1 = a1.begin();
        string_array::const_iterator e1 = a1.end();
        string_array::const_iterator i2 = a2.begin();
        while (i1 != e1) {
            if ((*i1) != (*i2)) {
                ret = false;
                break;
            }
            ++i1;
            ++i2;
        }
    }
    return ret;
}

inline void __add_to_array(string_array& a, const std::string& sSrc,
                           std::string::const_iterator iBeg, std::string::const_iterator iEnd)
{
    static const char szSpaces[] = CFUN_SPACES;
    if (iEnd > iBeg) {
        size_t nFirst = sSrc.find_first_not_of(szSpaces, iBeg - sSrc.begin());

        if (nFirst != CFUN_NOT_FOUND) {
            size_t nLast = sSrc.find_last_not_of(szSpaces, iEnd - sSrc.begin() - 1);

            if (nLast != CFUN_NOT_FOUND &&
                nLast >= nFirst) {
                a.push_back(std::string(sSrc.begin() + nFirst, sSrc.begin() + (nLast + 1)));
            }
        }
    }
}

inline void __parse_vars(const std::string& sVars, size_t nFirst, size_t nLast, string_array& saResult)
{
    size_t nL;
    __clean_array(saResult);
    for (;;) {
        nL = sVars.find(CFUN_COMMA, nFirst);
        if (nL == CFUN_NOT_FOUND || nL >= nLast) {
            __add_to_array(saResult, sVars, sVars.begin() + nFirst, sVars.begin() + nLast);
            break;
        }
        __add_to_array(saResult, sVars, sVars.begin() + nFirst, sVars.begin() + nL);
        nFirst = nL + 1;
    }
}

inline std::string __format_vars(const string_array& a, bool emptyBraces = true)
{
    std::stringstream stream;
    bool showBraces = emptyBraces || a.size() > 0;
    if (showBraces) {
        stream << CFUN_O_BRACE;
    }
    string_array::const_iterator it = a.begin();
    string_array::const_iterator end = a.end();
    while (it != end) {
        stream << *it;
        ++it;
        if (it != end)
            stream << CFUN_SCOMMA;
    }
    if (showBraces) {
        stream << CFUN_C_BRACE;
    }
    return stream.str();
}

// MAC is missing isalpha while linking...
inline int is_alpha(int c)
{
#if defined (MAC)
    if ((c >= 0141 && c <= 0172) || (c >= 0101 && c <= 0132)) {
        return 1;
    }
    return 0;
#else
    return isalpha(c);
#endif
}

// k*x & k=a+b ---> (a+b)*x
inline bool __subst_parameter(std::string& sBody, const string_array& saParameters,
                              const string_array& saMeanings) throw(cvmexception)
{
    bool ret = false;
    const size_t nPsz = saParameters.size();

    for (size_t i = 0; i < nPsz; i++) {
        if (sBody == saParameters[i]) {
            ret = true;
            if (i < saMeanings.size()) {
                const std::string& sMeaning = saMeanings[i];
                size_t nBeg = sMeaning.find_first_of(sBody);
                if (nBeg != CFUN_NOT_FOUND) {
                    size_t nEnd = nBeg + sBody.length();
                    if ((nBeg <= 0 || !is_alpha(sMeaning.at(nBeg - 1))) &&  // p -> {p+p}  ||  p -> {2*p-x}  ||  p -> {22-p}
                        (nEnd >= sMeaning.length() || !is_alpha(sMeaning.at(nEnd)))) {
//                  if ((nBeg <= 0 || !isalpha(sMeaning.at(nBeg - 1))) &&  // p -> {p+p}  ||  p -> {2*p-x}  ||  p -> {22-p}
//                      (nEnd >= sMeaning.length() || !isalpha(sMeaning.at(nEnd)))) {
                        throw cvmexception(CFUN_PARAMETER_RECURSION, sBody.c_str(), sMeaning.c_str());
                    }
                }
                sBody.assign(__trim_beg(sMeaning), __trim_end(sMeaning));
            }
            else {
                throw cvmexception(CFUN_SUBSTPARAMETERERROR, sBody.c_str());
            }
            break;
        }
    }
    return ret;
}

// returns true if body looks like 1.e+12 or 1.e-12
inline bool __ifbodynum(const std::string& sBody, size_t nBodyLength, size_t nSignPos)
{
    static const char ccExp = 'e', ccExpb = 'E', ccPoint = '.';
    return nSignPos < nBodyLength - 1 &&
           isdigit(sBody[nSignPos + 1]) &&
           nSignPos > 1 &&
           (sBody[nSignPos - 1] == ccExp || sBody[nSignPos - 1] == ccExpb) &&
           (sBody[nSignPos - 2] == ccPoint || isdigit(sBody[nSignPos - 2]));
}

inline size_t __tassign(std::string& sDest, const std::string& sSrc, size_t nFirst)
{
    return __tassign(sDest, sSrc, nFirst, sSrc.length());
}

inline size_t __tassign(std::string& sDest, const std::string& sSrc)
{
    return __tassign(sDest, sSrc, 0, sSrc.length());
}

template <typename T> class BaseFunction;
template <typename T> class UnaryFunction;
template <typename T> class BinaryFunction;
template <typename T> class Fconst;
template <typename T> class Finfinity;
template <typename T> class Fmninfinity;
template <typename T> class Fplus;
template <typename T> class Fminus;
template <typename T> class Fmult;
template <typename T> class Fdiv;
template <typename T> class Fpower;
template <typename T> class Fsat;
template <typename T> class Fexp;
template <typename T> class Fsqrt;
template <typename T> class Flog;
template <typename T> class Flog10;
template <typename T> class Fsin;
template <typename T> class Fcos;
template <typename T> class Ftan;
template <typename T> class Fasin;
template <typename T> class Facos;
template <typename T> class Fatan;
template <typename T> class Fsinh;
template <typename T> class Fcosh;
template <typename T> class Ftanh;
template <typename T> class Fsinint;
template <typename T> class Fcosint;
template <typename T> class Fuminus;
template <typename T> class Fsign;
template <typename T> class Fabs;
template <typename T> class Fdelta;
template <typename T> class Fvar;
template <typename T> class Fiif;

#ifdef CVM_USE_VARIADIC_TEMPLATES
#   define CFUN_NEW_FUNC(cls, ...) std::make_shared<cls<T>>(__VA_ARGS__)
#else
#   define CFUN_NEW_FUNC(cls, ...) BasePointer(new cls<T>(__VA_ARGS__))
#endif



template <typename T>
class FunctionFactory
{
public:
    typedef typename std::shared_ptr<BaseFunction<T>> BasePointer; 

    static BasePointer compile(const std::string& sPar, size_t nFirst, size_t nLast,
                               const string_array& saVars,
                               const string_array& saParameters,
                               const string_array& saMeanings,
                               bool subst = true) throw(cvmexception)
    {
        size_t i;
        static const char pccSigns[] = CFUN_SIGNS;
        static const char pccInfinity[] = CFUN_INF;
        std::string sBody, sLeft;
        size_t nBodyLength = __tassign(sBody, sPar, nFirst, nLast);

        if (!nBodyLength) {
            throw cvmexception(CFUN_PARSEERROR, sPar.c_str(), __format_vars(saVars).c_str());
        }

        // trying to handle real or complex number first
        BasePointer pfConst = CFUN_NEW_FUNC(Fconst, sBody, false); // std::make_shared<Fconst<T>>(sBody, false);
        if (pfConst->is_valid()) {
            return pfConst;
        }

        if (subst && __subst_parameter(sBody, saParameters, saMeanings)) { // k*x & k=a+b ---> (a+b)*x
            return FunctionFactory<T>::compile(sBody, 0, sBody.length(), saVars, saParameters, saMeanings, false);
        }

        int nLevel = 0;
        for (i = nBodyLength - 1; i < CFUN_NOT_FOUND; --i) {
            switch (sBody[i]) {
            case CFUN_O_PARENTH:
                ++nLevel;
                break;
            case CFUN_C_PARENTH:
                --nLevel;
                break;
            case CFUN_PLUS:
                if (!nLevel) {
                    if (!i) {                                              // +x -> x
                        return FunctionFactory<T>::compile(sBody, 1, nBodyLength, saVars, saParameters, saMeanings, subst);
                    }
                    if (__ifbodynum(sBody, nBodyLength, i)) {
                        break;
                    }
                    return CFUN_NEW_FUNC(Fplus, sBody, 0, i, i + 1, nBodyLength, saVars, saParameters, saMeanings);
                }
                break;
            case CFUN_MINUS:
                if (!nLevel) {
                    if (__ifbodynum(sBody, nBodyLength, i)) {
                        break;
                    }
                    if (!i) {
                        return CFUN_NEW_FUNC(Fuminus, sBody, 1, nBodyLength, saVars, saParameters, saMeanings);
                    }
                    return CFUN_NEW_FUNC(Fminus, sBody, 0, i, i + 1, nBodyLength, saVars, saParameters, saMeanings);
                }
                break;
            default:
                break;
            }
        }
        if (nLevel) {                                                       // ((...) or (...))
            throw cvmexception(CFUN_PARSEERROR, sBody.c_str(), __format_vars(saVars).c_str());
        }

        for (i = nBodyLength - 1; i < CFUN_NOT_FOUND; --i) {
            switch (sBody[i]) {
            case CFUN_O_PARENTH:
                ++nLevel;
                break;
            case CFUN_C_PARENTH:
                --nLevel;
                break;
            case CFUN_MULT:                                                 // x * y
                if (!i) {
                    throw cvmexception(CFUN_PARSEERROR, sBody.c_str(), __format_vars(saVars).c_str());
                }
                if (!nLevel) {
                    return CFUN_NEW_FUNC(Fmult, sBody, 0, i, i + 1, nBodyLength, saVars, saParameters, saMeanings);
                }
                break;
            case CFUN_DIV:                                                  // x / y
                if (!i) {
                    throw cvmexception(CFUN_PARSEERROR, sBody.c_str(), __format_vars(saVars).c_str());
                }
                if (!nLevel) {
                    return CFUN_NEW_FUNC(Fdiv, sBody, 0, i, i + 1, nBodyLength, saVars, saParameters, saMeanings);
                }
                break;
            default:
                break;
            }
        }

        for (i = 0; i < nBodyLength; ++i) {
            switch (sBody[i]) {
            case CFUN_O_PARENTH:
                ++nLevel;
                break;
            case CFUN_C_PARENTH:
                --nLevel;
                break;
            case CFUN_POWER:                                                // x ^ y
                if (!i) {
                    throw cvmexception(CFUN_PARSEERROR, sBody.c_str(), __format_vars(saVars).c_str());
                }
                if (!nLevel) {
                    return CFUN_NEW_FUNC(Fpower, sBody, 0, i, i + 1, nBodyLength, saVars, saParameters, saMeanings);
                }
                break;
            default:
                break;
            }
        }

        if (isdigit(sBody[0]) || sBody[0] == CFUN_POINT) {                  // 1.e-3 or .2e+12 or 99
            return CFUN_NEW_FUNC(Fconst, sBody);
        }
        if (sBody == pccInfinity) {
            return CFUN_NEW_FUNC(Finfinity);
        }
        for (i = 0; i < saVars.size(); ++i) {
            if (saVars[i] == sBody) {
                return std::make_shared<Fvar<T>>(sBody, i);
            }
        }

        size_t nLeft = sBody.find_first_of(pccSigns, 0);                    // sin(... cos(... etc. or x1 or y
        __tassign(sLeft, sBody, 0, nLeft);

        if (sLeft.length() == 0) {
            throw cvmexception(CFUN_PARSEERROR, sBody.c_str(), __format_vars(saVars).c_str());
        }

        if (sLeft == CFUN_EXP)
            return CFUN_NEW_FUNC(Fexp, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SQRT)
            return CFUN_NEW_FUNC(Fsqrt, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_LOG)
            return CFUN_NEW_FUNC(Flog, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_LOG10)
            return CFUN_NEW_FUNC(Flog10, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SIN)
            return CFUN_NEW_FUNC(Fsin, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_COS)
            return CFUN_NEW_FUNC(Fcos, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_TAN)
            return CFUN_NEW_FUNC(Ftan, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_ASIN)
            return CFUN_NEW_FUNC(Fasin, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_ACOS)
            return CFUN_NEW_FUNC(Facos, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_ATAN)
            return CFUN_NEW_FUNC(Fatan, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SINH)
            return CFUN_NEW_FUNC(Fsinh, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_COSH)
            return CFUN_NEW_FUNC(Fcosh, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_TANH)
            return CFUN_NEW_FUNC(Ftanh, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SI)
            return CFUN_NEW_FUNC(Fsinint, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_CI)
            return CFUN_NEW_FUNC(Fcosint, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SIGN)
            return CFUN_NEW_FUNC(Fsign, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_ABS)
            return CFUN_NEW_FUNC(Fabs, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_IIF)
            return CFUN_NEW_FUNC(Fiif, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_DELTA)
            return CFUN_NEW_FUNC(Fdelta, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_POWERS)
            return CFUN_NEW_FUNC(Fpower, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_SAT)
            return CFUN_NEW_FUNC(Fsat, sBody, nLeft, nBodyLength, saVars, saParameters, saMeanings);
        else if (sLeft == CFUN_I)
            return CFUN_NEW_FUNC(Fconst, T(1));
        else
            throw cvmexception(CFUN_PARSEERROR, sLeft.c_str(), __format_vars(saVars).c_str());
        return CFUN_NEW_FUNC(Fconst, T(0));
    }
};


template <typename T>
class BaseFunction
{
protected:
    typedef typename FunctionFactory<T>::BasePointer BasePointer; 

    virtual BasePointer _simpl() const = 0;                     // returns simplified function as a new object
    virtual int _depth(bool) const = 0;                         // call stack acquire and release

public:
    BaseFunction()
    {
    }

    virtual ~BaseFunction()
    {
    }

    virtual int _kind() const = 0;                              // a kind of function type id
    virtual T _value(const T* pdVars) const = 0;                // value of scalar function of few variables
    virtual BasePointer _clone() const = 0;                     // a new clone of an object
    virtual BasePointer _deriv(size_t nVarNum) const = 0;       // _derivative by given variable as a new object
    virtual bool _equals(const BasePointer& pfSrc) const = 0;   // checks whether they are identical
    virtual std::string _format(int nPrecision) const = 0;      // returns formatted body

    virtual bool is_valid() const {                             // to override in FConst
        return true;
    }

    virtual BasePointer _getArg(size_t) const {                 // nArgNum-th variable, 0-based, const version
        return nullptr;
    }

    BasePointer _simp() const
    {
        BasePointer pRet;
        if (_depth(true) > CFUN_SIMPS_STACK_DEPTH) {
            pRet = this->_clone();
#ifdef CFUN_DEBUG
            std::cout << "WARNING! Stack depth " << CFUN_SIMPS_STACK_DEPTH << " reached" << std::endl;
#endif
        }
        else {
            pRet = this->_simpl();
        }
        _depth(false);
        return pRet;
    }
};


template <typename T>
class UnaryFunction : public BaseFunction<T>
{
protected:
    typedef typename BaseFunction<T>::BasePointer BasePointer; 
    BasePointer mf0; 

    virtual BasePointer _simpl() const = 0;                     // returns simplified function as a new object
    virtual int _depth(bool) const = 0;                         // call stack acquire and release

public:
    UnaryFunction()
      : mf0()
    {
    }

    UnaryFunction(const UnaryFunction<T>& rfSrc)
      : mf0(rfSrc.mf0)
    {
    }

    UnaryFunction(UnaryFunction<T>&& rfSrc)
      : mf0(std::move(rfSrc.mf0))
    {
    }

    UnaryFunction(const BasePointer& pfSrc) 
      : mf0(pfSrc)
    {
    }

    UnaryFunction(BasePointer&& pfSrc)
      : mf0(std::move(pfSrc))
    {
    }

    UnaryFunction(const std::string& sArg, const string_array& saVars,
                  const string_array& saParameters, const string_array& saMeanings)
      : mf0(FunctionFactory<T>::compile(sArg, 0, sArg.length(), saVars, saParameters, saMeanings))
    {
    }

    UnaryFunction(const std::string& sArg, size_t nF, size_t nL, const string_array& saVars,
                  const string_array& saParameters, const string_array& saMeanings)
      : mf0(FunctionFactory<T>::compile(sArg, nF, nL, saVars, saParameters, saMeanings))
    {
    }

    UnaryFunction& operator = (const UnaryFunction<T>& rfSrc)
    {
        mf0 = rfSrc.mf0;
        return *this;
    }

    UnaryFunction& operator = (UnaryFunction<T>&& rfSrc)
    {
        mf0 = std::move(rfSrc.mf0);
        return *this;
    }

    virtual int _kind() const = 0;
    virtual T _value(const T* pdVars) const = 0;
    virtual BasePointer _clone() const = 0;
    virtual BasePointer _deriv(size_t nVarNum) const = 0;
    virtual const char* _name() const = 0;

    bool _equals(const BasePointer& pfSrc) const override
    {
        return pfSrc &&
               mf0 &&
               this->_kind() == pfSrc->_kind() &&
               mf0->_equals(pfSrc->_getArg(0));
    }

    BasePointer _getArg(size_t nArgNum) const override
    {
        return nArgNum ? nullptr : mf0;
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << this->_name() << CFUN_O_PARENTH << mf0->_format(nPrecision) << CFUN_C_PARENTH;
        return osBuf.str();
    }
};


template <typename T>
class BinaryFunction : public UnaryFunction<T>
{
protected:
    typedef typename UnaryFunction<T>::BasePointer BasePointer; 
    BasePointer mf1; 

    void BinCtr(const std::string& sBodyCommaBody,              // for bodies like pow(a,b) and sat(a,b)
                size_t nFirst, size_t nLast,
                const string_array& saVars,
                const string_array& saParameters, const string_array& saMeanings) throw(cvmexception)
    {
        std::string sLeft, sRight;
        if (__separate(sBodyCommaBody, nFirst, nLast, sLeft, sRight) == 1) {
            this->mf0 = FunctionFactory<T>::compile(sLeft, 0, sLeft.length(), saVars, saParameters, saMeanings);
            this->mf1 = FunctionFactory<T>::compile(sRight, 0, sRight.length(), saVars, saParameters, saMeanings);
        }
        else {
            throw cvmexception(CFUN_PARSEERROR, sBodyCommaBody.c_str(), __format_vars(saVars).c_str());
        }
    }

    virtual BasePointer _simpl() const = 0;                     // returns simplified function as a new object
    virtual int _depth(bool) const = 0;                         // call stack acquire and release

public:
    BinaryFunction()
      : mf1()
    {
    }

    BinaryFunction(const BinaryFunction<T>& rfSrc)
      : UnaryFunction<T>(rfSrc.mf0),
        mf1(rfSrc.mf1)
    {
    }

    BinaryFunction(BinaryFunction<T>&& rfSrc)
      : UnaryFunction<T>(std::move(rfSrc.mf0)),
        mf1(std::move(rfSrc.mf1))
    {
    }

    BinaryFunction(const BaseFunction<T>& rfSrc0, const BaseFunction<T>& rfSrc1)
      : UnaryFunction<T>(rfSrc0),
        mf1(rfSrc1.mf1)
    {
    }

    BinaryFunction(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : UnaryFunction<T>(pfSrc0),
        mf1(pfSrc1)
    {
    }

    BinaryFunction(const std::string& sArg1, const std::string& sArg2,
                   const string_array& saVars,
                   const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg1, 0, sArg1.length(), saVars, saParameters, saMeanings),
        mf1(FunctionFactory<T>::compile(sArg2, 0, sArg2.length(), saVars, saParameters, saMeanings))
    {
    }

    BinaryFunction(const std::string& sArgs, 
                   size_t nF1, size_t nL1, size_t nF2, size_t nL2,
                   const string_array& saVars,
                   const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArgs, nF1, nL1, saVars, saParameters, saMeanings),
        mf1(FunctionFactory<T>::compile(sArgs, nF2, nL2, saVars, saParameters, saMeanings))
    {
    }

    BinaryFunction(const std::string& sArg, size_t nFirst, size_t nLast,
                   const string_array& saVars,
                   const string_array& saParameters, const string_array& saMeanings)
        : mf1()
    {
        std::string sArg1, sArg2;
        if (__separate(sArg, nFirst, nLast, sArg1, sArg2) == 1) {
            this->mf0 = FunctionFactory<T>::compile(sArg1, 0, sArg1.length(), saVars, saParameters, saMeanings);
            this->mf1 = FunctionFactory<T>::compile(sArg2, 0, sArg2.length(), saVars, saParameters, saMeanings);
        }
    }

    BinaryFunction& operator = (const BinaryFunction<T>& rfSrc)
    {
        this->mf0 = rfSrc.mf0;
        this->mf1 = rfSrc.mf1;
        return *this;
    }

    BinaryFunction& operator = (BinaryFunction<T>&& rfSrc)
    {
        this->mf0 = std::move(rfSrc.mf0);
        this->mf1 = std::move(rfSrc.mf1);
        return *this;
    }

    virtual int _kind() const = 0;
    virtual T _value(const T* pdVars) const = 0;
    virtual BasePointer _clone() const = 0;
    virtual BasePointer _deriv(size_t nVarNum) const = 0;
    virtual std::string _format(int nPrecision) const = 0;

    bool _equals(const BasePointer& pfSrc) const override
    {
        return pfSrc && this->mf0 && this->mf1 &&
               this->_kind() == pfSrc->_kind() &&
               this->mf0->_equals(pfSrc->_getArg(0)) &&
               this->mf1->_equals(pfSrc->_getArg(1));
    }

    BasePointer _getArg(size_t nArgNum) const override
    {
        return nArgNum ? this->mf1 : this->mf0;
    }

    const char* _name() const override
    {
#ifdef _DEBUG
        assert(false);
#endif
        static const char sz_exp[] = "";         // no name for binaries
        return sz_exp;
    }
};


template <typename T>
class Fconst : public BaseFunction<T>
{
private:
    T mdConst;
    bool mbValid;

    int _kind() const override {
        return CFUN_K_CONST;
    }

protected:
    typedef typename BaseFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override {
        return this->_clone();
    }

    int _depth(bool) const override {
        return 1;
    }

public:
    Fconst()
      : mdConst(CFUN_M_ZERO),
        mbValid(true)
    {
    }

    Fconst(const Fconst& rfSrc)
      : BaseFunction<T>(),
        mdConst(rfSrc.mdConst),
        mbValid(rfSrc.mbValid)
    {
    }

    Fconst(const T& rdSrc)
      : mdConst(rdSrc),
        mbValid(true)
    {
    }

    Fconst(const std::string& sConst, bool bThrow = true)
      : mdConst(CFUN_M_ZERO), 
        mbValid(Helper<T>::convert(sConst, mdConst, bThrow))
    {
    }

    bool is_valid() const override {
        return mbValid;
    }

    T _value(const T*) const override {
        return mdConst;
    }

    BasePointer _clone() const override {
        return std::make_shared<Fconst<T>>(*this);
    }

    BasePointer _deriv(size_t) const override {
        return std::make_shared<Fconst<T>>();
    }

    std::string _format(int nPrecision) const override {
        return Helper<T>::_format(mdConst, nPrecision);
    }

    bool _equals(const BasePointer& pfSrc) const override {
        return this->_kind() == pfSrc->_kind() && Comparator<T,T>::eq(mdConst, pfSrc->_value(nullptr));
    }

protected:
    template<typename U>
    class Helper {
    public:
        static std::string _format(const U& v, int nPrecision)
        {
            std::ostringstream osBuf;
            if (nPrecision > 0) {
                osBuf.precision(nPrecision);
                osBuf.setf(std::ios::scientific | std::ios::showpoint | std::ios::left);
            }
            if (v < CFUN_M_ZERO) {
                osBuf << CFUN_O_PARENTH << v << CFUN_C_PARENTH;
            }
            else {
                osBuf << v;
            }
            return osBuf.str();
        }
        static bool convert(const std::string& s, U& u, bool bThrow) throw(cvmexception)
        {
            bool ret = __parse_num<U>(s, u);
            if (bThrow && !ret) {
                throw cvmexception(CFUN_PARSEERROR, s.c_str(), "n/a");
            }
            return ret;
        }
    };

    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static std::string _format(const UC& v, int nPrecision)
        {
            std::ostringstream osBuf;
            if (nPrecision > 0) {
                osBuf.precision(nPrecision);
                osBuf.setf(std::ios::scientific | std::ios::showpoint | std::ios::left);
            }
            osBuf << v;
            return osBuf.str();
        }
        static bool convert(const std::string& s, UC& uc, bool bThrow) throw(cvmexception)
        {
            bool ret1 = true, ret2 = true;
            std::string sRe, sIm;
            if (__separate(s, 0, s.length(), sRe, sIm) == 1) {
                U re, im;
                ret1 = __parse_num<U>(sRe, re);
                if (bThrow && !ret1) {
                    throw cvmexception(CFUN_PARSEERROR, sRe.c_str(), "n/a");
                }
                ret2 = __parse_num<U>(sIm, im);
                if (bThrow && !ret2) {
                    throw cvmexception(CFUN_PARSEERROR, sIm.c_str(), "n/a");
                }
                uc = UC(re, im);
            }
            else {
                U re;
                ret1 = __parse_num<U>(s, re);
                if (bThrow && !ret1) {
                    throw cvmexception(CFUN_PARSEERROR, s.c_str(), "n/a");
                }
                uc = UC(re, CFUN_M_ZERO);
            }
            return ret1 && ret2;
        }
    };
};

template <typename T>
class Finfinity : public BaseFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_INFINITY;
    }

protected:
    typedef typename BaseFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override {
        return this->_clone();
    }

    int _depth(bool) const override {
        return 1;
    }

public:
    Finfinity()
      : BaseFunction<T>()
    {
    }

    T _value(const T*) const override {
        return basic_cvmMachMax<T>();
    }

    BasePointer _clone() const {
        return std::make_shared<Finfinity<T>>();
    }

    BasePointer _deriv(size_t) const override {
        return this->_clone();
    }

    std::string _format(int) const override {
        return CFUN_SINF;
    }

    bool _equals(const BasePointer& pfSrc) const override {
        return this->_kind() == pfSrc->_kind();
    }
};

template <typename T>
class Fmninfinity : public BaseFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_MNINFINITY;
    }

protected:
    typedef typename BaseFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override {
        return this->_clone();
    }

    int _depth(bool) const override {
        return 1;
    }

public:
    Fmninfinity()
      : BaseFunction<T>()
    {
    }

    T _value(const T*) const override {
        return -basic_cvmMachMax<T>();
    }

    BasePointer _clone() const override {
        return std::make_shared<Fmninfinity<T>>();
    }

    BasePointer _deriv(size_t) const override {
        return this->_clone();
    }

    std::string _format(int) const override {
        return CFUN_SMINF;
    }

    bool _equals(const BasePointer& pfSrc) const override {
        return this->_kind() == pfSrc->_kind();
    }
};

template <typename T>
class Fplus : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_PLUS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        int i, j, k;
        BasePointer fZero = std::make_shared<Fconst<T>>();
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());

        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_equals(fZero))
            {                                                               // 0 + f(x) = f(x) | f(x) + 0 = f(x)
                return fCopy[__not(i)];
            }
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_INFINITY)
            {                                                               // inf + f(x) = inf | f(x) + inf = inf
                return std::make_shared<Finfinity<T>>();
            }
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_MNINFINITY)
            {                                                               // -inf + f(x) = -inf | f(x) + -inf = -inf
                return std::make_shared<Fmninfinity<T>>();
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_CONST &&
            fCopy[1]->_kind() == CFUN_K_CONST)
        {                                                                   // 1 + 2 = 3
            return std::make_shared<Fconst<T>>(fCopy[0]->_value(nullptr) + fCopy[1]->_value(nullptr));
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_CONST)
            {                                                               // x + (-1) = x-1
                T val = fCopy[i]->_value(nullptr);
                if (Comparator<T,T>::lt(val, CFUN_M_ZERO))
                {
                    return std::make_shared<Fminus<T>>(fCopy[__not(i)], std::make_shared<Fconst<T>>(-val))->_simp();
                }
            }
        }
        if (fCopy[0]->_equals(fCopy[1]))
        {                                                                   // x+x = 2*x
            return std::make_shared<Fmult<T>>(std::make_shared<Fconst<T>>(CFUN_M_TWO), fCopy[0])->_simp();
        }

        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_UMINUS)
            {                                                               // -a + b = b - a     a + -b = a - b
                return std::make_shared<Fminus<T>>(fCopy[__not(i)], fCopy[i]->_getArg(0))->_simp();
            }
        }

        if (fCopy[0]->_kind() == CFUN_K_MULT &&
            fCopy[1]->_kind() == CFUN_K_MULT)
        {
            for (i = 0; i <= 1; ++i)
            {                                                               // C*A+C*B = C*(A+B)
                for (j = 0; j <= 1; ++j)
                {                                                           // C*A+B*C = C*(A+B)
                    if (fCopy[0]->_getArg(i)->_equals(fCopy[1]->_getArg(j)))
                    {                                                       // A*C+C*B = C*(A+B)
                        return
                        std::make_shared<Fmult<T>>(fCopy[0]->_getArg(i),    // A*C+B*C = C*(A+B)
                                                   std::make_shared<Fplus<T>>(fCopy[0]->_getArg(__not(i)),
                                                                              fCopy[1]->_getArg(__not(j))))->_simp();
                    }
                }
            }
        }
        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_MULT)
            {                                                               // a*x + x = (a+1)*x
                for (i = 0; i <= 1; ++i)
                {
                    if (fCopy[j]->_getArg(i)->_equals(fCopy[__not(j)]))
                    {
                        return
                        std::make_shared<Fmult<T>>(fCopy[__not(j)],
                                                   std::make_shared<Fplus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                              std::make_shared<Fconst<T>>(CFUN_M_ONE)))->_simp();
                    }
                }
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_DIV &&
            fCopy[1]->_kind() == CFUN_K_DIV &&                              // x/a+y/a = (x+y)/a
            fCopy[0]->_getArg(1)->_equals(fCopy[1]->_getArg(1)))
        {
            return
            std::make_shared<Fdiv<T>>(std::make_shared<Fplus<T>>(fCopy[0]->_getArg(0),
                                                                 fCopy[1]->_getArg(0)),
                                      fCopy[0]->_getArg(1))->_simp();
        }
        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_PLUS)
            {                                                               // x+a + a
                for (i = 0; i <= 1; ++i)
                {
                    BasePointer fTmp =
                    std::make_shared<Fplus<T>>(fCopy[j]->_getArg(i),
                                               fCopy[__not(j)]);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fplus<T>>(fTmps,
                                                   fCopy[j]->_getArg(__not(i)))->_simp();
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_PLUS)
                {                                                           // x+a + x+b
                    for (i = 0; i <= 1; ++i)
                    {
                        for (k = 0; k <= 1; ++k)
                        {
                            BasePointer fTmp =
                            std::make_shared<Fplus<T>>(fCopy[j]->_getArg(i),
                                                       fCopy[__not(j)]->_getArg(k));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fplus<T>>(fTmps,
                                                           std::make_shared<Fplus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                      fCopy[__not(j)]->_getArg(__not(k))))->_simp();
                            }
                        }
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_MINUS)
                {                                                           // x+a + x-b | x-a + x+b
                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp =
                        std::make_shared<Fminus<T>>(fCopy[j]->_getArg(i),
                                                    fCopy[__not(j)]->_getArg(1));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fplus<T>>(fTmps,
                                                       std::make_shared<Fplus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                  fCopy[__not(j)]->_getArg(0)))->_simp();
                        }
                    }
                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp =
                        std::make_shared<Fplus<T>>(fCopy[j]->_getArg(i),
                                                   fCopy[__not(j)]->_getArg(0));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fplus<T>>(fTmps,
                                                       std::make_shared<Fminus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                   fCopy[__not(j)]->_getArg(1)))->_simp();
                        }
                    }
                }
            }
            if (fCopy[j]->_kind() == CFUN_K_MINUS)
            {
                {                                                           // a-x + a
                    BasePointer fTmp =
                    std::make_shared<Fplus<T>>(fCopy[j]->_getArg(0),
                                               fCopy[__not(j)]);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fminus<T>>(fTmps,
                                                    fCopy[j]->_getArg(1))->_simp();
                    }
                }
                {
                    BasePointer fTmp =
                    std::make_shared<Fminus<T>>(fCopy[__not(j)],            // x-a + a
                                                fCopy[j]->_getArg(1));
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fplus<T>>(fTmps,
                                                   fCopy[j]->_getArg(0))->_simp();
                    }
                }

                if (fCopy[__not(j)]->_kind() == CFUN_K_MINUS)
                {                                                           // x-a + x-b
                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp =
                        std::make_shared<Fplus<T>>(fCopy[j]->_getArg(i),
                                                   fCopy[__not(j)]->_getArg(i));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            BasePointer fTmp =
                            i ?
                            std::make_shared<Fminus<T>>(std::make_shared<Fplus<T>>(fCopy[j]->_getArg(0),
                                                                                   fCopy[__not(j)]->_getArg(0)),
                                                        fTmps)
                            :
                            std::make_shared<Fminus<T>>(fTmps,
                                                        std::make_shared<Fplus<T>>(fCopy[j]->_getArg(1),
                                                                                   fCopy[__not(j)]->_getArg(1)));
                            return fTmp->_simp();
                        }
                    }
                    BasePointer fTmp =
                    std::make_shared<Fminus<T>>(fCopy[j]->_getArg(0),       // x-a + b-x
                                                fCopy[__not(j)]->_getArg(1));
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fplus<T>>(std::make_shared<Fminus<T>>(fCopy[__not(j)]->_getArg(0),
                                                                               fCopy[j]->_getArg(1)),
                                                   fTmps)->_simp();
                    }
                }
            }
        }
        return std::make_shared<Fplus<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fplus(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fplus(const std::string& sAddend1, const std::string& sAddend2,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sAddend1, sAddend2, saVars, saParameters, saMeanings)
    {
    }

    Fplus(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return this->mf0->_value(pdVars) + this->mf1->_value(pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fplus<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fplus<T>>(this->mf0->_deriv(nVarNum),
                                   this->mf1->_deriv(nVarNum))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << this->mf0->_format(nPrecision) << CFUN_SPLUS;
        if (this->mf1->_kind() == CFUN_K_UMINUS) {
            osBuf << CFUN_O_PARENTH << this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
        }
        else {
            osBuf << this->mf1->_format(nPrecision);
        }
        return osBuf.str();
    }
};

template <typename T>
class Fminus : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_MINUS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        int i, j, k;
        BasePointer fZero = std::make_shared<Fconst<T>>();
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        if (fCopy[0]->_equals(fZero))
        {                                                                   // 0 - x = -x
            return std::make_shared<Fuminus<T>>(fCopy[1])->_simp();
        }
        if (fCopy[1]->_equals(fZero))                                       // x - 0 = x
        {
            return fCopy[0];
        }
        if (fCopy[0]->_kind() == CFUN_K_INFINITY)                           // inf - x = inf
        {
            return std::make_shared<Finfinity<T>>();
        }
        if (fCopy[1]->_kind() == CFUN_K_MNINFINITY)                         // x - -inf = inf
        {
            return std::make_shared<Finfinity<T>>();
        }
        if (fCopy[0]->_kind() == CFUN_K_MNINFINITY)                         // -inf - x = -inf
        {
            return std::make_shared<Fmninfinity<T>>();
        }
        if (fCopy[1]->_kind() == CFUN_K_INFINITY)                           // x - inf = -inf
        {
            return std::make_shared<Fmninfinity<T>>();
        }
        if (fCopy[0]->_kind() == CFUN_K_CONST &&
            fCopy[1]->_kind() == CFUN_K_CONST)                              // 3 - 2 = 1
        {
            return std::make_shared<Fconst<T>>(fCopy[0]->_value(nullptr)-
                                               fCopy[1]->_value(nullptr));
        }
        if (fCopy[0]->_equals(fCopy[1]))                                    // x - x = 0
        {
            return std::make_shared<Fconst<T>>();
        }
        if (fCopy[0]->_kind() == CFUN_K_UMINUS)                             // -a - b = -(a + b)
        {
            return
            std::make_shared<Fuminus<T>>(std::make_shared<Fplus<T>>(fCopy[0]->_getArg(0),
                                                                    fCopy[1])->_simp())->_simp();
        }
        if (fCopy[1]->_kind() == CFUN_K_UMINUS)                              // a - -b = a + b
        {
            return
            std::make_shared<Fplus<T>>(fCopy[0],
                                       fCopy[1]->_getArg(0))->_simp();
        }
        if (fCopy[0]->_kind() == CFUN_K_MULT &&
            fCopy[1]->_kind() == CFUN_K_MULT)
        {
            for (i = 0; i <= 1; ++i)                                        // C*A-C*B = C*(A-B)
            {
                for (j = 0; j <= 1; ++j)                                    // C*A-B*C = C*(A-B)
                {
                    if (fCopy[0]->_getArg(i)->_equals(fCopy[1]->_getArg(j))) // A*C-C*B = C*(A-B)
                    {
                        return
                        std::make_shared<Fmult<T>>(fCopy[0]->_getArg(i),    // A*C-B*C = C*(A-B)
                                                  std::make_shared<Fminus<T>>(fCopy[0]->_getArg(__not(i)),
                                                                              fCopy[1]->_getArg(__not(j))))->_simp();
                    }
                }
            }
        }
        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_MULT)
            {                                                               // a*x - x = (a-1)*x
                for (i = 0; i <= 1; ++i)
                {
                    if (fCopy[j]->_getArg(i)->_equals(fCopy[__not(j)]))
                    {
                        BasePointer fTmp =
                        std::make_shared<Fmult<T>>(fCopy[__not(j)],
                                                   j ?
                                                   std::make_shared<Fminus<T>>(std::make_shared<Fconst<T>>(CFUN_M_ONE),
                                                                               fCopy[j]->_getArg(__not(i)))
                                                   :
                                                   std::make_shared<Fminus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                               std::make_shared<Fconst<T>>(CFUN_M_ONE)));
                        return fTmp->_simp();
                    }
                }
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_DIV &&
            fCopy[1]->_kind() == CFUN_K_DIV &&                              // x/a-y/a = (x+y)/a
            fCopy[0]->_getArg(1)->_equals(fCopy[1]->_getArg(1)))
        {
            return
            std::make_shared<Fdiv<T>>(std::make_shared<Fminus<T>>(fCopy[0]->_getArg(0),
                                                                  fCopy[1]->_getArg(0)),
                                      fCopy[0]->_getArg(1))->_simp();
        }
        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_PLUS)                           // [x+a] - a  |  a - [x+a]
            {
                for (i = 0; i <= 1; ++i)
                {
                    BasePointer fTmp = std::make_shared<Fminus<T>>(fCopy[j]->_getArg(i),
                                                                   fCopy[__not(j)]);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        if (j)
                        {
                            return
                            std::make_shared<Fminus<T>>(fTmps,
                                                        fCopy[j]->_getArg(__not(i)))->_simp();
                        }
                        else {
                            return
                            std::make_shared<Fplus<T>>(fCopy[j]->_getArg(__not(i)),
                                                       fTmps)->_simp();
                        }
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_PLUS)                // x+a - x+b
                {
                    for (i = 0; i <= 1; ++i)
                    {
                        for (k = 0; k <= 1; ++k)
                        {
                            BasePointer fTmp = std::make_shared<Fminus<T>>(fCopy[j]->_getArg(i),
                                                                           fCopy[__not(j)]->_getArg(k));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fminus<T>>(fTmps,
                                                            std::make_shared<Fplus<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                       fCopy[__not(j)]->_getArg(__not(k))))->_simp();
                            }
                        }
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_MINUS)               // x+a - x-b | x-a - x+b
                {
                    if (!j)
                    {
                        return
                        std::make_shared<Fplus<T>>(fCopy[0],
                                                   std::make_shared<Fminus<T>>(fCopy[1]->_getArg(1),
                                                                               fCopy[1]->_getArg(0)))->_simp();
                    }
                    else                                                    // j==1, i.e.    x-a - x+b
                    {
                        for (i = 0; i <= 1; ++i)
                        {
                            BasePointer fTmp = std::make_shared<Fminus<T>>(fCopy[0]->_getArg(0),
                                                                           fCopy[1]->_getArg(i));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fminus<T>>(fTmps,
                                                            std::make_shared<Fplus<T>>(fCopy[0]->_getArg(1),
                                                                                       fCopy[1]->_getArg(__not(i))))->_simp();
                            }
                        }
                        for (i = 0; i <= 1; ++i)                            // a-x - x+b
                        {
                            BasePointer fTmp = std::make_shared<Fplus<T>>(fCopy[0]->_getArg(1),
                                                                          fCopy[1]->_getArg(i));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fminus<T>>(std::make_shared<Fminus<T>>(fCopy[0]->_getArg(0),
                                                                                        fCopy[1]->_getArg(__not(i))),
                                                            fTmps)->_simp();
                            }
                        }
                    }
                }
            }
            if (fCopy[j]->_kind() == CFUN_K_MINUS)                          // a-x - a  |  a - a-x
            {
                if (j)                                                      // a - a-x == a + x-a
                {
                    return
                    std::make_shared<Fplus<T>>(fCopy[0],
                                               std::make_shared<Fminus<T>>(fCopy[1]->_getArg(1),
                                                                           fCopy[1]->_getArg(0)))->_simp();
                }
                else
                {
                    {                                                       // j==0 i.e. a-x - a
                        BasePointer fTmp = std::make_shared<Fminus<T>>(fCopy[j]->_getArg(0),
                                                                       fCopy[__not(j)]);
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fminus<T>>(fTmps,
                                                        fCopy[j]->_getArg(1))->_simp();
                        }
                    }

                    {                                                       // x-a - a
                        BasePointer fTmp = std::make_shared<Fplus<T>>(fCopy[j]->_getArg(1),
                                                                      fCopy[__not(j)]);
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fminus<T>>(fCopy[j]->_getArg(0),
                                                        fTmps)->_simp();
                        }
                    }
                }
            }
        }
        return std::make_shared<Fminus<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fminus(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fminus(const std::string& sMinuend, const std::string& sSubtrahend,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sMinuend, sSubtrahend, saVars, saParameters, saMeanings)
    {
    }

    Fminus(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return this->mf0->_value(pdVars) - this->mf1->_value(pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fminus<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fminus<T>>(this->mf0->_deriv(nVarNum),
                                    this->mf1->_deriv(nVarNum))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << this->mf0->_format(nPrecision) << CFUN_SMINUS;
        if (this->mf1->_kind() == CFUN_K_UMINUS) {
            osBuf << CFUN_O_PARENTH << this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
        }
        else {
            osBuf << this->mf1->_format(nPrecision);
        }
        return osBuf.str();
    }
};

template <typename T>
class Fmult : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_MULT;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        int i, j, k;
        BasePointer fZero = std::make_shared<Fconst<T>>();
        BasePointer fOne = std::make_shared<Fconst<T>>(CFUN_M_ONE);
        BasePointer fMOne = std::make_shared<Fconst<T>>(CFUN_M_MONE);
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_INFINITY)
            {
                return std::make_shared<Finfinity<T>>();
            }
            if (fCopy[i]->_kind() == CFUN_K_MNINFINITY)
            {
                return std::make_shared<Fmninfinity<T>>();
            }
        }
        if (fCopy[0]->_equals(fZero) ||
            fCopy[1]->_equals(fZero))                                       // 0 * x = 0 | x * 0 = 0
        {
            return fZero;
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_equals(fOne)) {                                  // 1 * x = x
                return fCopy[__not(i)];
            }
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_equals(fMOne))                                   // -1 * x = -x
            {
                return
                std::make_shared<Fuminus<T>>(fCopy[__not(i)])->_simp();
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_CONST &&
            fCopy[1]->_kind() == CFUN_K_CONST)                              // 3 * 2 = 6
        {
            return std::make_shared<Fconst<T>>(fCopy[0]->_value(nullptr) * fCopy[1]->_value(nullptr));
        }
        if (fCopy[0]->_equals(fCopy[1]))                                    // x * x = x^2
        {
            return
            std::make_shared<Fpower<T>>(fCopy[0],
                                        std::make_shared<Fconst<T>>(CFUN_M_TWO))->_simp();
        }
        if (fCopy[0]->_kind() == CFUN_K_UMINUS &&
            fCopy[1]->_kind() == CFUN_K_UMINUS)                             // (-a)*(-b) = a*b)
        {
            return
            std::make_shared<Fmult<T>>(fCopy[0]->_getArg(0),
                                       fCopy[1]->_getArg(0))->_simp();
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_UMINUS)                         // (-a)*b = -a*b
            {
                return
                std::make_shared<Fuminus<T>>(std::make_shared<Fmult<T>>(fCopy[i]->_getArg(0),
                                                                        fCopy[__not(i)]))->_simp();
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_POWER &&                            // x^a*x^b=x^(a+b)
            fCopy[1]->_kind() == CFUN_K_POWER &&
            fCopy[0]->_getArg(0)->_equals(fCopy[1]->_getArg(0)))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[0]->_getArg(0),
                                        std::make_shared<Fplus<T>>(fCopy[0]->_getArg(1),
                                                                   fCopy[1]->_getArg(1)))->_simp();
        }
        if (fCopy[0]->_kind() == CFUN_K_EXP &&                              // exp(a)*exp(b)=exp(a+b)
            fCopy[1]->_kind() == CFUN_K_EXP)
        {
            return
            std::make_shared<Fexp<T>>(std::make_shared<Fplus<T>>(fCopy[0]->_getArg(0),
                                                                 fCopy[1]->_getArg(0)))->_simp();
        }
        if (fCopy[0]->_kind() == CFUN_K_POWER &&                            // x^a*x=x^(a+1)
            fCopy[0]->_getArg(0)->_equals(fCopy[1]))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[1],
                                        std::make_shared<Fplus<T>>(fCopy[0]->_getArg(1),
                                                                   std::make_shared<Fconst<T>>(CFUN_M_ONE)))->_simp();
        }
        if (fCopy[1]->_kind() == CFUN_K_POWER &&                            // x*x^a=x^(a+1)
            fCopy[1]->_getArg(0)->_equals(fCopy[0]))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[0],
                                        std::make_shared<Fplus<T>>(fCopy[1]->_getArg(1),
                                                                   fOne))->_simp();
        }

        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_MULT)                           // x*a * a
            {
                for (i = 0; i <= 1; ++i)
                {
                    BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(i),
                                                                  fCopy[__not(j)]);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fmult<T>>(fTmps,
                                                   fCopy[j]->_getArg(__not(i)))->_simp();
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_MULT)                // x*a * x*b
                {
                    for (i = 0; i <= 1; ++i)
                    {
                        for (k = 0; k <= 1; ++k)
                        {
                            BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(i),
                                                                          fCopy[__not(j)]->_getArg(k));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fmult<T>>(fTmps,
                                                           std::make_shared<Fmult<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                      fCopy[__not(j)]->_getArg(__not(k))))->_simp();
                            }
                        }
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_DIV)                 // x*a * x/b | x/a * x*b
                {

                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(i),
                                                                     fCopy[__not(j)]->_getArg(1));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fmult<T>>(fTmps,
                                                       std::make_shared<Fmult<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                  fCopy[__not(j)]->_getArg(0)))->_simp();
                        }
                    }
                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(i),
                                                                      fCopy[__not(j)]->_getArg(0));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fmult<T>>(fTmps,
                                                       std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                 fCopy[__not(j)]->_getArg(1)))->_simp();
                        }
                    }
                }
            }
            if (fCopy[j]->_kind() == CFUN_K_DIV)                            // a/x * a
            {
                {
                    BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(0),
                                                                  fCopy[__not(j)]);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fdiv<T>>(fTmps,
                                                  fCopy[j]->_getArg(1))->_simp();
                    }
                }
                {                                                           // x/a * a
                    BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[__not(j)],
                                                                 fCopy[j]->_getArg(1));
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fmult<T>>(fTmps,
                                                   fCopy[j]->_getArg(0))->_simp();
                    }
                }

                if (fCopy[__not(j)]->_kind() == CFUN_K_DIV)                 //   x/a * x/b
                {
                    for (i = 0; i <= 1; ++i)
                    {
                        BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(i),
                                                                      fCopy[__not(j)]->_getArg(i));
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            BasePointer fTmp =
                            i ?
                            std::make_shared<Fdiv<T>>(std::make_shared<Fmult<T>>(fCopy[j]->_getArg(0),
                                                                                 fCopy[__not(j)]->_getArg(0)),
                                                      fTmps)
                            :
                            std::make_shared<Fdiv<T>>(fTmps,
                                                      std::make_shared<Fmult<T>>(fCopy[j]->_getArg(1),
                                                                                 fCopy[__not(j)]->_getArg(1)));
                            return fTmp->_simp();
                        }
                    }                                                       // x/a * b/x
                    BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(0),
                                                                 fCopy[__not(j)]->_getArg(1));
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        return
                        std::make_shared<Fmult<T>>(std::make_shared<Fdiv<T>>(fCopy[__not(j)]->_getArg(0),
                                                                             fCopy[j]->_getArg(1)),
                                                   fTmps)->_simp();
                    }
                }
            }
        }
        return std::make_shared<Fmult<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fmult(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fmult(const std::string& sMultiplicand, const std::string& sMultiplier,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sMultiplicand, sMultiplier, saVars, saParameters, saMeanings)
    {
    }

    Fmult(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return this->mf0->_value(pdVars) * this->mf1->_value(pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fmult<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fplus<T>>(std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                                              this->mf1),
                                   std::make_shared<Fmult<T>>(this->mf0,
                                                              this->mf1->_deriv(nVarNum)))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        switch (this->mf0->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
            osBuf << CFUN_O_SPARENTH << this->mf0->_format(nPrecision) << CFUN_C_PARENTH << CFUN_SMULT;
            break;
        default:
            osBuf << this->mf0->_format(nPrecision) << CFUN_SMULT;
            break;
        }
        switch (this->mf1->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_UMINUS:
            osBuf << CFUN_O_SPARENTH << this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
            break;
        default:
            osBuf << this->mf1->_format(nPrecision);
            break;
        }
        return osBuf.str();
    }
};

template <typename T>
class Fdiv : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_DIV;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        int i, j, k;
        BasePointer fZero = std::make_shared<Fconst<T>>();
        BasePointer fOne = std::make_shared<Fconst<T>>(CFUN_M_ONE);
        BasePointer fMOne = std::make_shared<Fconst<T>>(CFUN_M_MONE);
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        if (fCopy[0]->_equals(fZero))                                       // 0 / x = 0
        {
            return fZero;
        }
        if (fCopy[1]->_equals(fZero))
        {                                                                   // x / 0 = inf
            return std::make_shared<Finfinity<T>>();
        }
        if (fCopy[1]->_kind() == CFUN_K_INFINITY ||
            fCopy[1]->_kind() == CFUN_K_MNINFINITY)                         // x / inf = 0 | x / -inf = 0
        {
            return fZero;
        }
        if (fCopy[1]->_equals(fOne))
        {                                                                   // x / 1 = x
            return fCopy[0];
        }
        if (fCopy[1]->_equals(fMOne))
        {                                                                   // x / -1 = -x
            return std::make_shared<Fuminus<T>>(fCopy[0]);
        }
        if (fCopy[0]->_kind() == CFUN_K_CONST &&
            fCopy[1]->_kind() == CFUN_K_CONST)
        {                                                                   // 6 / 2 = 3
            T dNomin = fCopy[0]->_value(nullptr);
            T dDenom = fCopy[1]->_value(nullptr);
            if (cvm::_abs(dDenom) <= cvm::cvmMachMin())
            {
                if (Comparator<T,T>::lt(dNomin, CFUN_M_ZERO))
                {
                    return std::make_shared<Fmninfinity<T>>();
                }
                else
                {
                    return std::make_shared<Finfinity<T>>();
                }
            }
            else
            {
                return std::make_shared<Fconst<T>>(dNomin / dDenom);
            }
        }
        if (fCopy[0]->_equals(fCopy[1]))
        {                                                                   // x / x = 1
            return std::make_shared<Fconst<T>>(CFUN_M_ONE);
        }
        if (fCopy[0]->_kind() == CFUN_K_UMINUS &&
            fCopy[1]->_kind() == CFUN_K_UMINUS)
        {                                                                   // (-a)/(-b) = a/b)
            return
            std::make_shared<Fdiv<T>>(fCopy[0]->_getArg(0),
                                      fCopy[1]->_getArg(0))->_simp();
        }
        for (i = 0; i <= 1; ++i)
        {
            if (fCopy[i]->_kind() == CFUN_K_UMINUS)
            {                                                               // (-a)/b = -a/b
                BasePointer fTmp =
                i ?
                std::make_shared<Fuminus<T>>(std::make_shared<Fdiv<T>>(fCopy[__not(i)],
                                                                       fCopy[i]->_getArg(0)))
                :
                std::make_shared<Fuminus<T>>(std::make_shared<Fdiv<T>>(fCopy[i]->_getArg(0),
                                                                       fCopy[__not(i)]));
                return fTmp->_simp();
            }
        }
        if (fCopy[0]->_kind() == CFUN_K_POWER &&                            // x^a/x^b=x^(a-b)
            fCopy[1]->_kind() == CFUN_K_POWER &&
            fCopy[0]->_getArg(0)->_equals(fCopy[1]->_getArg(0)))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[0]->_getArg(0),
                                        std::make_shared<Fminus<T>>(fCopy[0]->_getArg(1),
                                                                    fCopy[1]->_getArg(1)))->_simp();
        }
        if (fCopy[0]->_kind() == CFUN_K_POWER &&                            // x^a/x=x^(a-1)
            fCopy[0]->_getArg(0)->_equals(fCopy[1]))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[1],
                                        std::make_shared<Fminus<T>>(fCopy[0]->_getArg(1),
                                                                    std::make_shared<Fconst<T>>(CFUN_M_ONE)))->_simp();
        }
        if (fCopy[1]->_kind() == CFUN_K_POWER &&                            // x/x^a=x^(1-a)
            fCopy[1]->_getArg(0)->_equals(fCopy[0]))
        {
            return
            std::make_shared<Fpower<T>>(fCopy[0],
                                        std::make_shared<Fminus<T>>(std::make_shared<Fconst<T>>(CFUN_M_ONE),
                                                                    fCopy[1]->_getArg(1)))->_simp();
        }

        for (j = 0; j <= 1; ++j)
        {
            if (fCopy[j]->_kind() == CFUN_K_MULT)
            {                                                               // [x*a] / a  |  a / [x*a]
                for (i = 0; i <= 1; ++i)
                {
                    BasePointer pNum = fCopy[j]->_getArg(i);
                    BasePointer pDen = fCopy[__not(j)];
                    BasePointer fTmp = std::make_shared<Fdiv<T>>(j ? pDen : pNum,
                                                                 j ? pNum : pDen);
                    BasePointer fTmps(fTmp->_simp());
                    if (!fTmp->_equals(fTmps))
                    {
                        if (j)
                        {
                            return
                            std::make_shared<Fdiv<T>>(fTmps,
                                                      fCopy[j]->_getArg(__not(i)))->_simp();
                        }
                        else
                        {
                            return
                            std::make_shared<Fmult<T>>(fCopy[j]->_getArg(__not(i)),
                                                       fTmps)->_simp();
                        }
                    }
                }
                if (fCopy[__not(j)]->_kind() == CFUN_K_MULT)
                {                                                           //   x*a / x*b
                    for (i = 0; i <= 1; ++i)
                    {
                        for (k = 0; k <= 1; ++k)
                        {
                            BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(i),
                                                                         fCopy[__not(j)]->_getArg(k));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fdiv<T>>(fTmps,
                                                          std::make_shared<Fmult<T>>(fCopy[j]->_getArg(__not(i)),
                                                                                     fCopy[__not(j)]->_getArg(__not(k))))->_simp();
                            }
                        }
                    }
                }

                if (fCopy[__not(j)]->_kind() == CFUN_K_DIV)
                {                                                           // x*a / x/b | x/a / x*b
                    if (!j)
                    {
                        return
                        std::make_shared<Fmult<T>>(fCopy[0],
                                                   std::make_shared<Fdiv<T>>(fCopy[1]->_getArg(1),
                                                                             fCopy[1]->_getArg(0)))->_simp();
                    }
                    else
                    {                                                       // j==1, i.e. x/a / x*b
                        for (i = 0; i <= 1; ++i)
                        {
                            BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[0]->_getArg(0),
                                                                         fCopy[1]->_getArg(i));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fdiv<T>>(fTmps,
                                                          std::make_shared<Fmult<T>>(fCopy[0]->_getArg(1),
                                                                                     fCopy[1]->_getArg(__not(i))))->_simp();
                            }
                        }
                        for (i = 0; i <= 1; ++i)
                        {                                                   //  a/x / x*b
                            BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[0]->_getArg(1),
                                                                          fCopy[1]->_getArg(i));
                            BasePointer fTmps(fTmp->_simp());
                            if (!fTmp->_equals(fTmps))
                            {
                                return
                                std::make_shared<Fdiv<T>>(std::make_shared<Fdiv<T>>(fCopy[0]->_getArg(0),
                                                                                    fCopy[1]->_getArg(__not(i))),
                                                          fTmps)->_simp();
                            }
                        }
                    }
                }
            }
            if (fCopy[j]->_kind() == CFUN_K_DIV)
            {                                                               // a/x / a  |  a / a/x
                if (j)
                {                                                           // a / a/x == a * x/a
                    return
                    std::make_shared<Fmult<T>>(fCopy[0],
                                               std::make_shared<Fdiv<T>>(fCopy[1]->_getArg(1),
                                                                         fCopy[1]->_getArg(0)))->_simp();
                }
                else
                {
                    {                                                       // j==0 i.e. a/x / a
                        BasePointer fTmp = std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(0),
                                                                     fCopy[__not(j)]);
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fdiv<T>>(fTmps,
                                                      fCopy[j]->_getArg(1))->_simp();
                        }
                    }
                    {
                        BasePointer fTmp = std::make_shared<Fmult<T>>(fCopy[j]->_getArg(1),    // x/a / a
                                                                      fCopy[__not(j)]);
                        BasePointer fTmps(fTmp->_simp());
                        if (!fTmp->_equals(fTmps))
                        {
                            return
                            std::make_shared<Fdiv<T>>(fCopy[j]->_getArg(0),
                                                      fTmps)->_simp();
                        }
                    }
                }
            }
        }
        return std::make_shared<Fdiv<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fdiv(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fdiv(const std::string& sDividend, const std::string& sDivisor,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sDividend, sDivisor, saVars, saParameters, saMeanings)
    {
    }

    Fdiv(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override
    {
        T dDenom = this->mf1->_value(pdVars);
        if (Comparator<T,T>::eq(dDenom, CFUN_M_ZERO)) {
            throw cvmexception(CVM_DIVISIONBYZERO);
        }
        return this->mf0->_value(pdVars) / dDenom;
    }

    BasePointer _clone() const override {
        return std::make_shared<Fdiv<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(std::make_shared<Fminus<T>>(std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                                                                         this->mf1),
                                                              std::make_shared<Fmult<T>>(this->mf0,
                                                                                         this->mf1->_deriv(nVarNum))),
                                  std::make_shared<Fpower<T>>(this->mf1,
                                                              std::make_shared<Fconst<T>>(CFUN_M_TWO)))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        switch (this->mf0->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_DIV:
            osBuf << CFUN_O_SPARENTH << this->mf0->_format(nPrecision) << CFUN_C_PARENTH << CFUN_SDIV;
            break;
        default:
            osBuf << this->mf0->_format(nPrecision) << CFUN_SDIV;
            break;
        }
        switch (this->mf1->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_UMINUS:
        case CFUN_K_MULT:
        case CFUN_K_DIV:
            osBuf << CFUN_O_SPARENTH << this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
            break;
        default:
            osBuf << this->mf1->_format(nPrecision);
            break;
        }
        return osBuf.str();
    }
};

template <typename T>
class Fpower : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_POWER;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        BasePointer fZero = std::make_shared<Fconst<T>>();
        BasePointer fOne = std::make_shared<Fconst<T>>(CFUN_M_ONE);
        BasePointer fMOne = std::make_shared<Fconst<T>>(CFUN_M_MONE);
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        if (fCopy[1]->_equals(fZero))
        {                                                                   // x^0 = 1
            return fOne;
        }
        else if (fCopy[0]->_equals(fZero))
        {                                                                   // 0^x = 0
            return fZero;
        }
        else if (fCopy[1]->_equals(fOne))
        {                                                                   // x^1 = x
            return fCopy[0];
        }
        else if (fCopy[0]->_kind() == CFUN_K_CONST &&
                 fCopy[1]->_kind() == CFUN_K_CONST)
        {                                                                   // 6^2 = 36
            return
            std::make_shared<Fconst<T>>(std::pow(fCopy[0]->_value(nullptr), fCopy[1]->_value(nullptr)));
        }
        else if (fCopy[1]->_equals(fMOne))
        {                                                                   // x^-1 = 1/x
            return
            std::make_shared<Fdiv<T>>(fOne,
                                      fCopy[0]);
        }
        else if (fCopy[0]->_kind() == CFUN_K_POWER)
        {                                                                   // [x^a]^b = x^[a*b]
            return
            std::make_shared<Fpower<T>>(fCopy[0]->_getArg(0),
                                        std::make_shared<Fmult<T>>(fCopy[0]->_getArg(1),
                                                                   fCopy[1]))->_simp();
        }
        return std::make_shared<Fpower<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fpower(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fpower(const std::string& sBasis, const std::string& sExponent,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sBasis, sExponent, saVars, saParameters, saMeanings)
    {
    }

    Fpower(const std::string& sBasisCommaExponent, size_t nFirst, size_t nLast,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
    {
        this->BinCtr(sBasisCommaExponent, nFirst, nLast, saVars, saParameters, saMeanings);
    }

    Fpower(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return std::pow(this->mf0->_value(pdVars), this->mf1->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fpower<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fplus<T>>(std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                                              std::make_shared<Fmult<T>>(this->mf1,
                                                                                         std::make_shared<Fpower<T>>(this->mf0,
                                                                                                                     std::make_shared<Fplus<T>>(this->mf1,
                                                                                                                                                std::make_shared<Fconst<T>>(CFUN_M_MONE))))),
                                   std::make_shared<Fmult<T>>(this->mf1->_deriv(nVarNum),
                                                              std::make_shared<Fmult<T>>(std::make_shared<Fpower<T>>(this->mf0,
                                                                                                                     this->mf1),
                                                                                         std::make_shared<Flog<T>>(this->mf0))))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        switch (this->mf0->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_MULT:
        case CFUN_K_DIV:
        case CFUN_K_POWER:
        case CFUN_K_UMINUS:
            osBuf << CFUN_O_SPARENTH << this->mf0->_format(nPrecision) << CFUN_C_SPARENTH << CFUN_SPOWER;
            break;
        default:
            osBuf << this->mf0->_format(nPrecision) << CFUN_SPOWER;
            break;
        }
        switch (this->mf1->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_MULT:
        case CFUN_K_DIV:
        case CFUN_K_POWER:
        case CFUN_K_UMINUS:
            osBuf << CFUN_O_SPARENTH << this->mf1->_format(nPrecision) << CFUN_C_SPARENTH;
            break;
        default:
            osBuf << this->mf1->_format(nPrecision);
            break;
        }
        return osBuf.str();
    }
};

template <typename T>
class Fsat : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_SAT;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        BasePointer fZero = std::make_shared<Fconst<T>>();
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        if (fCopy[0]->_equals(fZero)) {
            return fZero;
        }
        return std::make_shared<Fsat<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fsat(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fsat(const std::string& sFun, const std::string& sDelta,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sFun, sDelta, saVars, saParameters, saMeanings)
    {
    }

    Fsat(const std::string& sFunCommaDelta, size_t nFirst, size_t nLast,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
    {
        this->BinCtr(sFunCommaDelta, nFirst, nLast, saVars, saParameters, saMeanings);
    }

    Fsat(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::_value(this->mf0.get(), this->mf1.get(), pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsat<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t) const override
    {
        return
        std::make_shared<Fplus<T>>(std::make_shared<Fdelta<T>>(this->mf0,
                                                               this->mf1),
                                   std::make_shared<Fdelta<T>>(this->mf0,
                                                               std::make_shared<Fuminus<T>>(this->mf1)))->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << CFUN_SAT << CFUN_O_PARENTH << this->mf0->_format(nPrecision) << CFUN_SCOMMA <<
            this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
        return osBuf.str();
    }

protected:
    template<typename U>
    class Helper {
    public:
        static U _value(const BaseFunction<U>* pfArg0, BaseFunction<U>* pfArg1, const U* pdVars)
        {
            U dRes = CFUN_M_ZERO;
            const U dMeans = pfArg0->_value(pdVars);
            const U dDelta = cvm::_abs(pfArg1->_value(pdVars));
            if (dMeans > dDelta) {
                dRes = CFUN_M_ONE;
            }
            else if (dMeans < -dDelta) {
                dRes = CFUN_M_MONE;
            }
            return dRes;
        }
    };

    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static UC _value(const BaseFunction<UC>* pfArg0, BaseFunction<UC>* pfArg1, const UC* pdVars)
        {
            UC dRes = CFUN_M_ZERO;
            const U dMeans = pfArg0->_value(pdVars).real();
            const U dDelta = cvm::_abs(pfArg1->_value(pdVars).real());
            if (dMeans > dDelta) {
                dRes = CFUN_M_ONE;
            }
            else if (dMeans < -dDelta) {
                dRes = CFUN_M_MONE;
            }
            return dRes;
        }
    };
};

template <typename T>
class Fexp : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_EXP;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(std::exp(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fexp<T>>(fCopy);
    }

public:
    Fexp(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fexp(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fexp(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return std::exp(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fexp<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                   std::make_shared<Fexp<T>>(this->mf0))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_exp[] = CFUN_EXP;
        return sz_exp;
    }
};

template <typename T>
class Fsqrt : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_SQRT;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::sqrt(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fsqrt<T>>(fCopy);
    }

public:
    Fsqrt(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fsqrt(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fsqrt(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::sqrt(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsqrt<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fmult<T>>(std::make_shared<Fconst<T>>(CFUN_M_TWO),
                                                             std::make_shared<Fsqrt<T>>(this->mf0)))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_sqrt[] = CFUN_SQRT;
        return sz_sqrt;
    }
};

template <typename T>
class Flog : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_LOG;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fZero = std::make_shared<Fconst<T>>();
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_equals(fZero))
        {                                                                   // log 0 = - inf
            return std::make_shared<Fmninfinity<T>>();
        }
        else if (fCopy->_kind() == CFUN_K_CONST)
        {                                                                   // log 1 = 0
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::log(fCopy->_value(nullptr)));
        }
        else if (fCopy->_kind() == CFUN_K_SQRT)
        {                                                                   // log sqrt (x) = 0.5 * log (x)
            return
            std::make_shared<Fmult<T>>(std::make_shared<Flog<T>>(fCopy->_getArg(0)),
                                       std::make_shared<Fconst<T>>(CFUN_M_HALF))->_simp();
        }
        else if (fCopy->_kind() == CFUN_K_POWER)
        {                                                                   // log (x^y) = y * log (x)
            return
            std::make_shared<Fmult<T>>(std::make_shared<Flog<T>>(fCopy->_getArg(0)),
                                       fCopy->_getArg(1))->_simp();
        }
        return std::make_shared<Flog<T>>(fCopy);
    }

public:
    Flog(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Flog(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Flog(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::log(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Flog<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  this->mf0)->_simp();
    }

    const char* _name() const override
    {
        static const char sz_log[] = CFUN_LOG;
        return sz_log;
    }
};

template <typename T>
class Flog10 : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_LOG10;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fZero = std::make_shared<Fconst<T>>();
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_equals(fZero))
        {
            return std::make_shared<Fmninfinity<T>>();
        }
        else if (fCopy->_kind() == CFUN_K_CONST)
        {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::log10(fCopy->_value(nullptr)));
        }
        else if (fCopy->_kind() == CFUN_K_SQRT)
        {                                                                   // log sqrt (x) = 0.5 * log (x)
            return
            std::make_shared<Fmult<T>>(std::make_shared<Flog10<T>>(fCopy->_getArg(0)),
                                       std::make_shared<Fconst<T>>(CFUN_M_HALF))->_simp();
        }
        else if (fCopy->_kind() == CFUN_K_POWER)
        {                                                                   // log (x^y) = y * log (x)
            return
            std::make_shared<Fmult<T>>(std::make_shared<Flog10<T>>(fCopy->_getArg(0)),
                                       fCopy->_getArg(1))->_simp();
        }
        return std::make_shared<Flog10<T>>(fCopy);
    }

public:
    Flog10(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Flog10(const std::string& sArg,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Flog10(const std::string& sArg, size_t nF, size_t nL,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::log10(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Flog10<T>>(*this);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fmult<T>>(std::make_shared<Fconst<T>>(CFUN_M_LN_10),
                                                             this->mf0))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_log10[] = CFUN_LOG10;
        return sz_log10;
    }
};

template <typename T>
class Fsin : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_SIN;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(std::sin(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fsin<T>>(fCopy);
    }

public:
    Fsin(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fsin(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fsin(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return std::sin(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsin<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                   std::make_shared<Fcos<T>>(this->mf0))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_sin[] = CFUN_SIN;
        return sz_sin;
    }
};

template <typename T>
class Fcos : public UnaryFunction<T> {
private:
    int _kind() const override {
        return CFUN_K_COS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(std::cos(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fcos<T>>(fCopy);
    }

public:
    Fcos(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fcos(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fcos(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return std::cos(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fcos<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                   std::make_shared<Fuminus<T>>(std::make_shared<Fsin<T>>(this->mf0)))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_cos[] = CFUN_COS;
        return sz_cos;
    }
};

template <typename T>
class Ftan : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_TAN;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::tan(fCopy->_value(nullptr)));
        }
        return std::make_shared<Ftan<T>>(fCopy);
    }

public:
    Ftan(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Ftan(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Ftan(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::tan(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Ftan<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fpower<T>>(std::make_shared<Fcos<T>>(this->mf0),
                                                              std::make_shared<Fconst<T>>(CFUN_M_TWO)))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_tan[] = CFUN_TAN;
        return sz_tan;
    }
};

template <typename T>
class Fasin : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_ARCSIN;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::asin(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fasin<T>>(fCopy);
    }

public:
    Fasin(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fasin(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fasin(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::asin(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fasin<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fsqrt<T>>(std::make_shared<Fminus<T>>(std::make_shared<Fconst<T>>(CFUN_M_ONE),
                                                                                         std::make_shared<Fpower<T>>(this->mf0,
                                                                                                                     std::make_shared<Fconst<T>>(CFUN_M_TWO)))))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_asin[] = CFUN_ASIN;
        return sz_asin;
    }
};

template <typename T>
class Facos : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_ARCCOS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::acos(fCopy->_value(nullptr)));
        }
        return std::make_shared<Facos<T>>(fCopy);
    }

public:
    Facos(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Facos(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Facos(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::acos(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Facos<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fuminus<T>>(std::make_shared<Fasin<T>>(this->mf0)->_deriv(nVarNum))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_acos[] = CFUN_ACOS;
        return sz_acos;
    }
};

template <typename T>
class Fatan : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_ARCTAN;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::atan(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fatan<T>>(fCopy);
    }

public:
    Fatan(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fatan(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fatan(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::atan(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fatan<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fplus<T>>(std::make_shared<Fconst<T>>(CFUN_M_ONE),
                                                             std::make_shared<Fpower<T>>(this->mf0,
                                                                                         std::make_shared<Fconst<T>>(CFUN_M_TWO))))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_atan[] = CFUN_ATAN;
        return sz_atan;
    }
};


template <typename T>
class Fsinh : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_SINH;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::sinh(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fsinh<T>>(fCopy);
    }

public:
    Fsinh(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fsinh(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fsinh(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::sinh(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsinh<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                   std::make_shared<Fcosh<T>>(this->mf0))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_sinh[] = CFUN_SINH;
        return sz_sinh;
    }
};


template <typename T>
class Fcosh : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_COSH;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::cosh(fCopy->_value(nullptr)));
        }
        return std::make_shared<Fcosh<T>>(fCopy);
    }

public:
    Fcosh(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fcosh(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fcosh(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::cosh(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fcosh<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(this->mf0->_deriv(nVarNum),
                                   std::make_shared<Fuminus<T>>(std::make_shared<Fsinh<T>>(this->mf0)))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_cosh[] = CFUN_COSH;
        return sz_cosh;
    }
};


template <typename T>
class Ftanh : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_TANH;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(ElementaryFunctions<T>::tanh(fCopy->_value(nullptr)));
        }
        return std::make_shared<Ftanh<T>>(fCopy);
    }

public:
    Ftanh(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Ftanh(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Ftanh(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return ElementaryFunctions<T>::tanh(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Ftanh<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fdiv<T>>(this->mf0->_deriv(nVarNum),
                                  std::make_shared<Fpower<T>>(std::make_shared<Fcosh<T>>(this->mf0),
                                                              std::make_shared<Fconst<T>>(CFUN_M_TWO)))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_tanh[] = CFUN_TANH;
        return sz_tanh;
    }
};

template <typename T>
class Fsinint : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_INTSIN;
    }

protected:
    int _depth(bool bAcquire) const override {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(Helper<T>::_value(fCopy.get(), nullptr));
        }
        return std::make_shared<Fsinint<T>>(fCopy);
    }

public:
    Fsinint(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fsinint(const std::string& sArg,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fsinint(const std::string& sArg, size_t nF, size_t nL,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::_value(this->mf0.get(), pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsinint<T>>(this->mf0);
    }

    BasePointer _deriv(size_t) const override
    {
        return
        std::make_shared<Fdiv<T>>(std::make_shared<Fsin<T>>(this->mf0),
                                  this->mf0)->_simp();
    }

    const char* _name() const override
    {
        static const char sz_si[] = CFUN_SI;
        return sz_si;
    }

protected:
    template<typename U>
    class Helper {
    public:
        static U _value(const BaseFunction<U>* pfArg, const U* pdVars)
        {
            return ElementaryFunctions<U>::sinint(pfArg->_value(pdVars), basic_cvmMachSp<U>());
        }
    };
    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static UC _value(const BaseFunction<UC>* pfArg, const UC* pdVars)
        {
            return ElementaryFunctions<UC>::sinint(pfArg->_value(pdVars), basic_cvmMachSp<U>());
        }
    };
};

template <typename T>
class Fcosint : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_INTCOS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(Helper<T>::_value(fCopy.get(), nullptr));
        }
        return std::make_shared<Fcosint<T>>(fCopy);
    }

public:
    Fcosint(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fcosint(const std::string& sArg,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fcosint(const std::string& sArg, size_t nF, size_t nL,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::_value(this->mf0.get(), pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fcosint<T>>(this->mf0);
    }

    BasePointer _deriv(size_t) const override
    {
        return
        std::make_shared<Fdiv<T>>(std::make_shared<Fcos<T>>(this->mf0),
                                  this->mf0)->_simp();
    }

    const char* _name() const override
    {
        static const char sz_ci[] = CFUN_CI;
        return sz_ci;
    }

protected:
    template<typename U>
    class Helper {
    public:
        static U _value(const BaseFunction<U>* pfArg, const U* pdVars)
        {
            return ElementaryFunctions<U>::cosint(pfArg->_value(pdVars), basic_cvmMachSp<U>());
        }
    };
    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static UC _value(const BaseFunction<UC>* pfArg, const UC* pdVars)
        {
            return ElementaryFunctions<UC>::cosint(pfArg->_value(pdVars), basic_cvmMachSp<U>());
        }
    };
};

template <typename T>
class Fuminus : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_UMINUS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_UMINUS)
        {                                                                   //  - - b
            return fCopy->_getArg(0);
        }
        if (fCopy->_kind() == CFUN_K_INFINITY)
        {
            return std::make_shared<Fmninfinity<T>>();
        }
        if (fCopy->_kind() == CFUN_K_MNINFINITY)
        {
            return std::make_shared<Finfinity<T>>();
        }
        if (fCopy->_kind() == CFUN_K_CONST)
        {
            return std::make_shared<Fconst<T>>(- fCopy->_value(nullptr));
        }
        if (fCopy->_kind() == CFUN_K_MINUS)
        {                                                                   // - (a - b) = b - a
            return
            std::make_shared<Fminus<T>>(fCopy->_getArg(1),
                                        fCopy->_getArg(0));
        }
        if (fCopy->_kind() == CFUN_K_MULT)
        {                                                                   // - x*y
            for (size_t i = 0; i <= 1; ++i)
            {
                BasePointer fTmp = std::make_shared<Fuminus<T>>(fCopy->_getArg(i));
                BasePointer fTmps(fTmp->_simp());
                if (!fTmp->_equals(fTmps))
                {
                    return
                    std::make_shared<Fmult<T>>(fTmps,
                                               fCopy->_getArg(__not(i)))->_simp();
                }
            }
        }
        if (fCopy->_kind() == CFUN_K_DIV)
        {                                                                   // - x/y
            for (size_t i = 0; i <= 1; ++i)
            {
                BasePointer fTmp = std::make_shared<Fuminus<T>>(fCopy->_getArg(i));
                BasePointer fTmps(fTmp->_simp());
                if (!fTmp->_equals(fTmps))
                {
                    BasePointer fTmp(
                    i ?
                    std::make_shared<Fdiv<T>>(fCopy->_getArg(__not(i)),
                                              fTmps)
                    :
                    std::make_shared<Fdiv<T>>(fTmps,
                                              fCopy->_getArg(__not(i))));
                    return fTmp->_simp();
                }
            }
        }
        return std::make_shared<Fuminus<T>>(fCopy);
    }

public:
    Fuminus(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fuminus(const std::string& sArg,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fuminus(const std::string& sArg, size_t nF, size_t nL,
            const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return - this->mf0->_value(pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fuminus<T>>(this->mf0);
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << this->_name();
        switch (this->mf0->_kind()) {
        case CFUN_K_PLUS:
        case CFUN_K_MINUS:
        case CFUN_K_UMINUS:
            osBuf << CFUN_O_SPARENTH << this->mf0->_format(nPrecision) << CFUN_C_PARENTH;
            break;
        default:
            osBuf << this->mf0->_format(nPrecision);
            break;
        }
        return osBuf.str();
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fuminus<T>>(this->mf0->_deriv(nVarNum))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_uminus[] = CFUN_SMINUS;
        return sz_uminus;
    }
};

template <typename T>
class Fsign : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_SIGN;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST)
        {
            return std::make_shared<Fconst<T>>(Helper<T>::_value(fCopy.get(), nullptr));
        }
        return std::make_shared<Fsign<T>>(fCopy);
    }

public:
    Fsign(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fsign(const std::string& sArg,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fsign(const std::string& sArg, size_t nF, size_t nL,
          const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::_value(this->mf0.get(), pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fsign<T>>(this->mf0);
    }

    BasePointer _deriv(size_t) const override
    {
        return
        std::make_shared<Fdelta<T>>(this->mf0,
                                    std::make_shared<Fconst<T>>())->_simp();
    }

    const char* _name() const override
    {
        static const char sz_sign[] = CFUN_SIGN;
        return sz_sign;
    }

protected:
    template<typename U>
    class Helper {
    public:
        static U _value(const BaseFunction<U>* pfArg, const U* pdVars)
        {
            const U dV = pfArg->_value(pdVars);
            return dV > CFUN_M_ZERO ? CFUN_M_ONE :
                                      (dV < CFUN_M_ZERO ? CFUN_M_MONE : CFUN_M_ZERO);
        }
    };

    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static UC _value(const BaseFunction<UC>* pfArg, const UC* pdVars)
        {
            const U dV = pfArg->_value(pdVars).real();
            return dV > CFUN_M_ZERO ? CFUN_M_ONE :
                                      (dV < CFUN_M_ZERO ? CFUN_M_MONE : CFUN_M_ZERO);
        }
    };
};

template <typename T>
class Fabs : public UnaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_ABS;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    BasePointer _simpl() const override
    {
        BasePointer fCopy(this->mf0->_simp());
        if (fCopy->_kind() == CFUN_K_CONST) {
            return std::make_shared<Fconst<T>>(cvm::_abs(fCopy->_value(nullptr)));
        }
        else if (fCopy->_kind() == CFUN_K_ABS) {
            return fCopy;
        }
        return std::make_shared<Fabs<T>>(fCopy);
    }

public:
    Fabs(const BasePointer& pfSrc)
      : UnaryFunction<T>(pfSrc)
    {
    }

    Fabs(const std::string& sArg,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, saVars, saParameters, saMeanings)
    {
    }

    Fabs(const std::string& sArg, size_t nF, size_t nL,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : UnaryFunction<T>(sArg, nF, nL, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return cvm::_abs(this->mf0->_value(pdVars));
    }

    BasePointer _clone() const override {
        return std::make_shared<Fabs<T>>(this->mf0);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fmult<T>>(std::make_shared<Fsign<T>>(this->mf0),
                                   this->mf0->_deriv(nVarNum))->_simp();
    }

    const char* _name() const override
    {
        static const char sz_abs[] = CFUN_ABS;
        return sz_abs;
    }
};

template <typename T>
class Fdelta : public BinaryFunction<T>
{
private:
    int _kind() const override {
        return CFUN_K_DELTA;
    }

protected:
    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    typedef typename BinaryFunction<T>::BasePointer BasePointer; 
    typedef shared_ptr_pair<BaseFunction<T>> PointerPair;

    BasePointer _simpl() const override
    {
        PointerPair fCopy(this->mf0->_simp(), this->mf1->_simp());
        if (fCopy[0]->_kind() == CFUN_K_CONST && fCopy[1]->_kind() == CFUN_K_CONST) {
            T dM1 = fCopy[0]->_value(nullptr);
            T dM2 = fCopy[1]->_value(nullptr);
            if (cvm::_abs(dM1 - dM2) <= cvm::cvmMachMin()) {
                return std::make_shared<Finfinity<T>>();
            }
            else {
                return std::make_shared<Fconst<T>>();
            }
        }
        return std::make_shared<Fdelta<T>>(fCopy[0], fCopy[1]);
    }

public:
    Fdelta(const BasePointer& pfSrc0, const BasePointer& pfSrc1)
      : BinaryFunction<T>(pfSrc0, pfSrc1)
    {
    }

    Fdelta(const std::string& sArg1, const std::string& sArg2,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArg1, sArg2, saVars, saParameters, saMeanings)
    {
    }

    Fdelta(const std::string& sArgs, size_t nF1, size_t nL1, size_t nF2, size_t nL2,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArgs, nF1, nL1, nF2, nL2, saVars, saParameters, saMeanings)
    {
    }

    Fdelta(const std::string& sArg, size_t nFirst, size_t nLast,
           const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : BinaryFunction<T>(sArg, nFirst, nLast, saVars, saParameters, saMeanings)
    {
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::_value(this->mf0.get(), this->mf1.get(), pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fdelta<T>>(this->mf0, this->mf1);
    }

    BasePointer _deriv(size_t) const override
    {
        return
        std::make_shared<Fdelta<T>>(this->mf0,
                                    this->mf1)->_simp();
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << CFUN_DELTA << CFUN_O_PARENTH << this->mf0->_format(nPrecision) << CFUN_SCOMMA <<
            this->mf1->_format(nPrecision) << CFUN_C_PARENTH;
        return osBuf.str();
    }

protected:
    template<typename U>
    class Helper {
    public:
        static U _value(const BaseFunction<U>* pfArg0, const BaseFunction<U>* pfArg1, const U* pdVars)
        {
            return cvm::_abs(pfArg0->_value(pdVars) - pfArg1->_value(pdVars)) <= basic_cvmMachMin<U>() ?
                   basic_cvmMachMax<U>() : CFUN_M_ZERO;
        }
    };

    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static UC _value(const BaseFunction<UC>* pfArg0, const BaseFunction<UC>* pfArg1, const UC* pdVars)
        {
            return cvm::_abs((pfArg0->_value(pdVars) - pfArg1->_value(pdVars)).real()) <= basic_cvmMachMin<U>() ?
                   basic_cvmMachMax<U>() : CFUN_M_ZERO;
        }
    };
};

template <typename T>
class Fvar : public BaseFunction<T>
{
private:
    size_t mnVarPos;
    std::string msVarName;

    int _kind() const override {
        return CFUN_K_VAR;
    }

protected:
    int _depth(bool) const override {
        return 1;
    }

    typedef typename BaseFunction<T>::BasePointer BasePointer; 


    BasePointer _simpl() const override
    {
        return this->_clone();
    }

public:
    Fvar()
      : mnVarPos(0),
        msVarName()
    {
    }

    Fvar(const Fvar<T>& rfSrc)
      : mnVarPos(rfSrc.mnVarPos),
        msVarName(rfSrc.msVarName)
    {
    }

    Fvar(std::string& sVN, size_t nVP)
      : mnVarPos(nVP),
        msVarName(sVN)
    {
    }

    T _value(const T* pdVars) const override
    {
        if (pdVars == nullptr) {
            throw cvmexception(CFUN_NULLPOINTERERROR, "Fvar<T>::_value");
        }
        return pdVars[mnVarPos];
    }

    BasePointer _clone() const override {
        return std::make_shared<Fvar<T>>(*this);
    }

    std::string _format(int) const override {
        return msVarName;
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        if (mnVarPos == nVarNum) {
            return std::make_shared<Fconst<T>>(CFUN_M_ONE);
        }
        return std::make_shared<Fconst<T>>();
    }

    bool _equals(const BasePointer& pfSrc) const override
    {
        return this->_kind() == pfSrc->_kind() &&
               // since kind() is the same we can do this downcasting
               mnVarPos == reinterpret_cast<const Fvar<T>*>(pfSrc.get())->mnVarPos &&
               msVarName == reinterpret_cast<const Fvar<T>*>(pfSrc.get())->msVarName;
    }

    BasePointer _getArg(size_t) const override {
        return nullptr;
    }
};

template <typename T>
class Fiif : public UnaryFunction<T>
{
private:
    typedef typename UnaryFunction<T>::BasePointer BasePointer; 

    int _kind() const override {
        return CFUN_K_IIF;
    }

    BasePointer mpflt;  // if mf0 < 0
    BasePointer mpfge;  // if mf0 >= 0

protected:
    void TriCtr(const std::string& sBodyCommaBodyCommaBody, size_t nFirst, size_t nLast,
                const string_array& saVars,
                const string_array& saParameters, const string_array& saMeanings) throw(cvmexception)
    {
        std::string sF, sLeft, sRight;
        if (__separate(sBodyCommaBodyCommaBody, nFirst, nLast, sF, sRight) == 1) {
            this->mf0 = FunctionFactory<T>::compile(sF, 0, sF.length(), saVars, saParameters, saMeanings);
            if (__separate(sRight, 0, sRight.length(), sLeft, sRight) == 1) {
                mpflt = FunctionFactory<T>::compile(sLeft, 0, sLeft.length(), saVars, saParameters, saMeanings);
                mpfge = FunctionFactory<T>::compile(sRight, 0, sRight.length(), saVars, saParameters, saMeanings);
                return;
            }
        }
        throw cvmexception(CFUN_PARSEERROR, sBodyCommaBodyCommaBody.c_str(), __format_vars(saVars).c_str());
    }

    int _depth(bool bAcquire) const override
    {
        static int nDepth;
        return bAcquire ? ++nDepth : --nDepth;
    }

    BasePointer _simpl() const override
    {
        return std::make_shared<Fiif<T>>(this->mf0->_simp(), this->mpflt->_simp(), this->mpfge->_simp());
    }

public:
    Fiif(const BasePointer& rfSrc, const BasePointer& rfLT, const BasePointer& rfGE)
      : UnaryFunction<T>(rfSrc),
        mpflt(rfLT),
        mpfge(rfGE)
    {
    }

    Fiif(const std::string& sBody,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : mpflt(nullptr),
        mpfge(nullptr)
    {
        TriCtr(sBody, 0, sBody.length(), saVars, saParameters, saMeanings);
    }

    Fiif(const std::string& sBody, size_t nFirst, size_t nLast,
         const string_array& saVars, const string_array& saParameters, const string_array& saMeanings)
      : mpflt(nullptr),
        mpfge(nullptr)
    {
        TriCtr(sBody, nFirst, nLast, saVars, saParameters, saMeanings);
    }

    T _value(const T* pdVars) const override {
        return Helper<T>::left(this->mf0.get(), pdVars) ? this->mpflt->_value(pdVars) :
                                                          this->mpfge->_value(pdVars);
    }

    BasePointer _clone() const override {
        return std::make_shared<Fiif<T>>(*this);
    }

    BasePointer _deriv(size_t nVarNum) const override
    {
        return
        std::make_shared<Fiif<T>>(this->mf0, this->mpflt->_deriv(nVarNum), this->mpfge->_deriv(nVarNum));
    }

    std::string _format(int nPrecision) const override
    {
        std::ostringstream osBuf;
        osBuf << this->_name() << CFUN_O_PARENTH << this->mf0->_format(nPrecision) << CFUN_SCOMMA <<
            this->mpflt->_format(nPrecision) << CFUN_SCOMMA <<
            this->mpfge->_format(nPrecision) << CFUN_C_PARENTH;
        return osBuf.str();
    }

    const char* _name() const override
    {
        static const char sz_iif[] = CFUN_IIF;
        return sz_iif;
    }

protected:
    template<typename U>
    class Helper {
    public:
        static bool left(const BaseFunction<U>* pfArg0, const U* pdVars)
        {
            return pfArg0->_value(pdVars) < CFUN_M_ZERO;
        }
    };

    template<typename U>
    class Helper<std::complex<U> > {
        typedef std::complex<U> UC;
    public:
        static bool left(const BaseFunction<UC>* pfArg0, const UC* pdVars)
        {
            return pfArg0->_value(pdVars).real() < CFUN_M_ZERO;
        }
    };
};



template <typename T>
class basic_function
{
public:
    typedef typename FunctionFactory<T>::BasePointer BasePointer; 

protected:
    mutable string_array mvars; 
    BasePointer mp; 

    const string_array& _check_vars(const string_array& saVars) const throw(cvmexception)
    {
        if (this->vars_num() > 0 && saVars.size() > 0 && !__arrays_equal(mvars, saVars)) {
            throw cvmexception(CFUN_VARSDONTMATCH, __format_vars(mvars).c_str(), __format_vars(saVars).c_str());
        }
        return this->vars_num() > 0 ? mvars : saVars;      // here we cover cases like '5+x'
    }

    // protected utility constructor
    basic_function(const string_array& saVars, const BasePointer& pf)
      : mvars(saVars),
        mp(pf)
    {
    }

    basic_function(string_array&& saVars, BasePointer&& pf)
      : mvars(std::move(saVars)),
        mp(std::move(pf))
    {
    }

public:
    basic_function()
      : mvars(),
        mp(std::make_shared<Fconst<T>>())
    {
    }

    explicit basic_function(const string_array& saVars)
      : mvars(saVars),
        mp(std::make_shared<Fconst<T>>())
    {
    }

    explicit basic_function(const T& dConst)
      : mvars(),
        mp(std::make_shared<Fconst<T>>(dConst))
    {
    }

    basic_function(const basic_function& rf)
      : mvars(rf.mvars),
        mp(rf.mp)
    {
    }

    basic_function(basic_function&& rf)
      : mvars(std::move(rf.mvars)),
        mp(std::move(rf.mp))
    {
    }

    explicit basic_function(const std::string& sInput)
      : mvars(),
        mp()
    {
        string_array saParameters, saMeanings;
        size_t nAF = sInput.find(CFUN_O_BRACE, 0);
        size_t nAL = sInput.find(CFUN_C_BRACE, 0);
        if (nAL != CFUN_NOT_FOUND && nAL > nAF) {
            __parse_vars(sInput, nAF + 1, nAL, mvars);
        }
        mp = FunctionFactory<T>::compile(sInput, nAL + 1, sInput.length(), mvars, saParameters, saMeanings);
    }

    basic_function(const string_array& saVars, const std::string& sBody)
      : mvars(saVars),
        mp()
    {
        string_array saParameters, saMeanings;
        mp = FunctionFactory<T>::compile(sBody, 0, sBody.length(), mvars, saParameters, saMeanings);
    }

    basic_function(const string_array& saVars, const std::string& sBody,
                   const string_array& saParameters, const string_array& saMeanings)
        : mvars(saVars),
          mp(FunctionFactory<T>::compile(sBody, 0, sBody.length(), mvars, saParameters, saMeanings))
    {
    }

    virtual ~basic_function()
    {
    }

    const string_array& vars() const {
        return mvars;
    }

    size_t vars_num() const {
        return mvars.size();
    }

    T value(const T* pd) const {
        return mp->_value(pd);
    }


    T operator () (const T* pd) const {
        return mp->_value(pd);
    }

    T operator () (const T& d) const {
        return mp->_value(&d);
    }

    T operator () (const T& d1, const T& d2) const
    {
        T dVals[2];
        dVals[0] = d1;
        dVals[1] = d2;
        return mp->_value(dVals);
    }

    T operator () (const T& d1, const T& d2, const T& d3) const
    {
        T dVals[3];
        dVals[0] = d1;
        dVals[1] = d2;
        dVals[2] = d3;
        return mp->_value(dVals);
    }

    T operator () () const
    {
        static const T dZero = T();
        return mp->_value(&dZero);
    }

    std::string format(int nPrecision = 0) const {
        return mp->_format(nPrecision);
    }

    std::string vformat(int nPrecision = 0) const
    {
        std::stringstream stream;
        std::string vars = __format_vars(this->mvars, false);
        if (vars.size() > 0) {
            stream << vars << CFUN_SSPACE;
        }
        stream << mp->_format(nPrecision);
        return stream.str();
    }

    bool operator == (const basic_function& rf) const {
        return vars_num() == rf.vars_num() && mp->_equals(rf.mp);
    }

    bool operator != (const basic_function& rf) const {
        return !operator == (rf);
    }

    basic_function drv(size_t nVarNum = 0) const {
        return basic_function(mvars, mp->_deriv(nVarNum));
    }

    basic_function& simp()
    {
        mp = mp->_simp();
        return *this;
    }

    basic_function& operator = (const basic_function& rfSrc)
    {
        if (mp != rfSrc.mp) {
            // we don't check vars here, we just override...
            mvars = rfSrc.mvars;
            mp = rfSrc.mp;
        }
        return *this;
    }

    basic_function& operator = (basic_function&& rfSrc)
    {
        if (mp != rfSrc.mp) {
            // we don't check vars here, we just override...
            mvars = std::move(rfSrc.mvars);
            mp = std::move(rfSrc.mp);
        }
        return *this;
    }

    basic_function& operator = (const T& rdSrc)
    {
        mp = std::make_shared<Fconst<T>>(rdSrc);
        return *this;
    }

    basic_function operator - () const {
        return basic_function(mvars, std::make_shared<Fuminus<T>>(mp));
    }

    basic_function operator + (const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fplus<T>>(mp, rf.mp));
    }

    basic_function operator - (const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fminus<T>>(mp, rf.mp));
    }

    basic_function operator *(const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fmult<T>>(mp, rf.mp));
    }

    basic_function operator / (const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fdiv<T>>(mp, rf.mp));
    }

    basic_function operator ^(const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fpower<T>>(mp, rf.mp));
    }

    basic_function& operator += (const basic_function& rf) throw(cvmexception)
    {
        mvars = this->_check_vars(rf.mvars);
        mp = std::make_shared<Fplus<T>>(mp, rf.mp);
        return *this;
    }

    basic_function& operator -= (const basic_function& rf) throw(cvmexception)
    {
        mvars = this->_check_vars(rf.mvars);
        mp = std::make_shared<Fminus<T>>(mp, rf.mp);
        return *this;
    }

    basic_function& operator *= (const basic_function& rf) throw(cvmexception)
    {
        mvars = this->_check_vars(rf.mvars);
        mp = std::make_shared<Fmult<T>>(mp, rf.mp);
        return *this;
    }

    basic_function& operator /= (const basic_function& rf) throw(cvmexception)
    {
        mvars = this->_check_vars(rf.mvars);
        mp = std::make_shared<Fdiv<T>>(mp, rf.mp);
        return *this;
    }

    basic_function operator + (const T& d) const throw(cvmexception) {
        return basic_function(mvars, std::make_shared<Fplus<T>>(mp, std::make_shared<Fconst<T>>(d)));
    }

    basic_function operator - (const T& d) const throw(cvmexception) {
        return basic_function(mvars, std::make_shared<Fminus<T>>(mp, std::make_shared<Fconst<T>>(d)));
    }

    basic_function operator *(const T& d) const throw(cvmexception) {
        return basic_function(mvars, std::make_shared<Fmult<T>>(mp, std::make_shared<Fconst<T>>(d)));
    }

    basic_function operator / (const T& d) const throw(cvmexception) {
        return basic_function(mvars, std::make_shared<Fdiv<T>>(mp, std::make_shared<Fconst<T>>(d)));
    }

    basic_function operator ^(const T& d) const throw(cvmexception) {
        return basic_function(mvars, std::make_shared<Fpower<T>>(mp, std::make_shared<Fconst<T>>(d)));
    }

    basic_function& operator += (const T& d) throw(cvmexception)
    {
        mp = std::make_shared<Fplus<T>>(mp, std::make_shared<Fconst<T>>(d));
        return *this;
    }

    basic_function& operator -= (const T& d) throw(cvmexception)
    {
        mp = std::make_shared<Fminus<T>>(mp, std::make_shared<Fconst<T>>(d));
        return *this;
    }

    basic_function& operator *= (const T& d) throw(cvmexception)
    {
        mp = std::make_shared<Fmult<T>>(mp, std::make_shared<Fconst<T>>(d));
        return *this;
    }

    basic_function& operator /= (const T& d) throw(cvmexception)
    {
        mp = std::make_shared<Fdiv<T>>(mp, std::make_shared<Fconst<T>>(d));
        return *this;
    }

    basic_function& operator ^= (const T& d) throw(cvmexception)
    {
        mp = std::make_shared<Fpower<T>>(mp, std::make_shared<Fconst<T>>(d));
        return *this;
    }

    basic_function sat(const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fsat<T>>(mp, rf.mp));
    }

    basic_function exp() const {
        return basic_function(mvars, std::make_shared<Fexp<T>>(mp));
    }

    basic_function sqrt() const {
        return basic_function(mvars, std::make_shared<Fsqrt<T>>(mp));
    }

    basic_function log() const {
        return basic_function(mvars, std::make_shared<Flog<T>>(mp));
    }

    basic_function log10() const {
        return basic_function(mvars, std::make_shared<Flog10<T>>(mp));
    }

    basic_function sin() const {
        return basic_function(mvars, std::make_shared<Fsin<T>>(mp));
    }

    basic_function cos() const {
        return basic_function(mvars, std::make_shared<Fcos<T>>(mp));
    }

    basic_function tan() const {
        return basic_function(mvars, std::make_shared<Ftan<T>>(mp));
    }

    basic_function asin() const {
        return basic_function(mvars, std::make_shared<Fasin<T>>(mp));
    }

    basic_function acos() const {
        return basic_function(mvars, std::make_shared<Facos<T>>(mp));
    }

    basic_function atan() const {
        return basic_function(mvars, std::make_shared<Fatan<T>>(mp));
    }

    basic_function sinh() const {
        return basic_function(mvars, std::make_shared<Fsinh<T>>(mp));
    }

    basic_function cosh() const {
        return basic_function(mvars, std::make_shared<Fcosh<T>>(mp));
    }

    basic_function tanh() const {
        return basic_function(mvars, std::make_shared<Ftanh<T>>(mp));
    }

    basic_function sinint() const {
        return basic_function(mvars, std::make_shared<Fsinint<T>>(mp));
    }

    basic_function cosint() const {
        return basic_function(mvars, std::make_shared<Fcosint<T>>(mp));
    }

    basic_function sign() const {
        return basic_function(mvars, std::make_shared<Fsign<T>>(mp));
    }

    basic_function abs() const {
        return basic_function(mvars, std::make_shared<Fabs<T>>(mp));
    }

    basic_function delta(const basic_function& rf) const throw(cvmexception) {
        return basic_function(this->_check_vars(rf.mvars), std::make_shared<Fdelta<T>>(mp, rf.mp));
    }

    basic_function iif (const basic_function& rfNeg,
                        const basic_function& rfNotNeg) const throw(cvmexception)
    {
        const string_array& saVarsT = this->_check_vars(rfNeg.mvars);
        const string_array& saVarsF = this->_check_vars(rfNotNeg.mvars);
        const string_array& saVars = saVarsT.size() > saVarsF.size() ? saVarsT : saVarsF;
        return basic_function(saVars, std::make_shared<Fiif<T>>(mp, rfNeg.mp, rfNotNeg.mp));
    }

    friend basic_function<T> operator + (const T& d, const basic_function<T>& rf) {
        return rf + d;
    }

    friend basic_function<T> operator - (const T& d, const basic_function<T>& rf) {
        return basic_function(rf.mvars, std::make_shared<Fminus<T>>(std::make_shared<Fconst<T>>(d), rf.mp));
    }

    friend basic_function<T> operator *(const T& d, const basic_function<T>& rf) {
        return rf * d;
    }

    friend basic_function<T> operator / (const T& d, const basic_function<T>& rf) {
        return basic_function(rf.mvars, std::make_shared<Fdiv<T>>(std::make_shared<Fconst<T>>(d), rf.mp));
    }

    friend basic_function<T> operator ^(const T& d, const basic_function<T>& rf) {
        return basic_function(rf.mvars, std::make_shared<Fpower<T>>(std::make_shared<Fconst<T>>(d), rf.mp));
    }

    friend std::ostream& operator <<(std::ostream& os, const basic_function& rf)
    {
        std::string vars = __format_vars(rf.mvars, false);
        if (vars.size() > 0) {
            os << vars << CFUN_SSPACE;
        }
        os << rf.format(0);
        return os;
    }

    // sets caller to be equal to scalar product of fa1 and fa2
    basic_function& _scalar_product(size_t nSize,
                                    const basic_function* fa1, size_t incr1,
                                    const basic_function* fa2, size_t incr2) throw(cvmexception)
    {
        const string_array& saVars1 = this->_check_vars(fa1->vars());
        const string_array& saVars2 = this->_check_vars(fa2->vars());
        const string_array& saVars = saVars1.size() > saVars2.size() ? saVars1 : saVars2;

        mvars = saVars;
        mp = std::make_shared<Fconst<T>>();

        for (size_t i = 0; i < nSize; ++i, fa1 += incr1, fa2 += incr2)
        {
            *this += (*fa1) * (*fa2);
        }
        return this->simp();
    }
};


template<typename BaseFunction>
inline void _copy(size_t nSize,
                  const BaseFunction* faFrom, size_t incrFrom,
                  BaseFunction* faTo, size_t incrTo)
{
    for (size_t i = 0; i < nSize; ++i, faFrom += incrFrom, faTo += incrTo) {
        *faTo = *faFrom;
    }
}


template<typename T>
class FArray
{
protected:
    typedef basic_function<T> BaseFunction; 

    std::vector<BaseFunction> mv; 

    void _assign(const FArray& a)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i);
        }
    }

public:
    FArray()
      : mv()
    {
    }

    explicit FArray(size_t nSize)
      : mv(nSize)
    {
    }

    FArray(const FArray& a)
      : mv(a.size())
    {
        this->_assign(a);
    }

    FArray(FArray&& a)
      : mv(std::move(a.mv))
    {
    }

    explicit FArray(const string_array& saInput) 
      : mv(saInput.size())
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = BaseFunction(saInput[i]);
        }
    }

    FArray(const string_array& saVars, const string_array& saBodies,
           const string_array& saParameters, const string_array& saMeanings)
      : mv(saBodies.size())
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = BaseFunction(saVars, saBodies[i], saParameters, saMeanings);
        }
    }

    virtual ~FArray()
    {
    }

    std::vector<BaseFunction>& impl() {
        return mv;
    }

    const std::vector<BaseFunction>& impl() const {
        return mv;
    }

    size_t size() const {
        return mv.size();
    }

    const BaseFunction& operator [] (size_t n) const {
        return mv.at(n);
    }

    BaseFunction& operator [] (size_t n) {
        return mv.at(n);
    }

    void value(T* pv) const {
        value(nullptr, pv);
    }

    void value(T d, T* pv) const {
        value(&d, pv);
    }

    void value(T d1, T d2, T* pv) const {
        T d[2];
        d[0] = d1;
        d[1] = d2;
        value(d, pv);
    }

    void value(T d1, T d2, T d3, T* pv) const {
        T d[3];
        d[0] = d1;
        d[1] = d2;
        d[2] = d3;
        value(d, pv);
    }

    void value(const T* pd, T* pv) const
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            pv[i] = mv.at(i).value(pd);
        }
    }

protected:
    void _drv(const FArray& a, size_t nVarNum)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i).drv(nVarNum);
        }
    }

    void _simp()
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i).simp();
        }
    }

    bool _equals(const FArray& a) const
    {
        bool bResult = true;
        if (mv.size() != a.size()) {
            bResult = false;
        }
        else
            for (size_t i = 0; i < mv.size(); ++i) {
                if (mv.at(i) != a.mv.at(i)) {
                    bResult = false;
                    break;
                }
            }
        return bResult;
    }

    void _add(const FArray& a)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) += a.mv.at(i);
        }
    }

    void _subtract(const FArray& a)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) -= a.mv.at(i);
        }
    }

    void _add(const FArray& a1, const FArray& a2)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a1.mv.at(i) + a2.mv.at(i);
        }
    }

    void _subtract(const FArray& a1, const FArray& a2)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a1.mv.at(i) - a2.mv.at(i);
        }
    }

    void _mult(const BaseFunction& f)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) *= f;
        }
    }

    void _div(const BaseFunction& f)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) /= f;
        }
    }

    void _mult(const FArray& a, const BaseFunction& f)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i) * f;
        }
    }

    void _div(const FArray& a, const BaseFunction& f)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i) / f;
        }
    }

    void _mult(const T& t)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) *= t;
        }
    }

    void _div(const T& t)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) /= t;
        }
    }

    void _mult(const FArray& a, const T& t)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i) * t;
        }
    }

    void _div(const FArray& a, const T& t)
    {
        for (size_t i = 0; i < mv.size(); ++i) {
            mv.at(i) = a.mv.at(i) / t;
        }
    }
};


template<typename T> class basic_fvector;
template<typename T> class basic_fmatrix;


template<typename T>
class basic_fvector : public FArray<T>
{
    friend class basic_fmatrix<T>;  // for _mult

protected:
    typedef FArray<T> BaseFArray; 
    typedef typename BaseFArray::BaseFunction BaseFunction; 

    void _check_size(size_t size) const throw(cvmexception)
    {
        _check_ne(CVM_SIZESMISMATCH, size, this->size());
    }

public:
    basic_fvector()
      : BaseFArray()
    {
    }

    explicit basic_fvector(size_t nSize)
      : BaseFArray(nSize)
    {
    }

    explicit basic_fvector(const string_array& saInput)
      : BaseFArray(saInput)
    {
    }

    basic_fvector(const string_array& saVars, const string_array& saBodies,
                  const string_array& saParameters, const string_array& saMeanings)
      : BaseFArray(saVars, saBodies, saParameters, saMeanings)
    {
    }

    basic_fvector(const basic_fvector& fv)
      : BaseFArray(fv)
    {
    }

    basic_fvector(basic_fvector&& fv)
      : BaseFArray(std::move(fv))
    {
    }

    basic_fvector& operator = (const basic_fvector& fv) throw(cvmexception)
    {
        this->_check_size(fv.size());
        this->_assign(fv);
        return *this;
    }

    basic_fvector& operator = (basic_fvector&& fv) throw(cvmexception)
    {
        this->_check_size(fv.size());
        this->mv = std::move(fv.mv);
        return *this;
    }

    bool operator == (const basic_fvector& fv) const {
        return this->_equals(fv);
    }

    bool operator != (const basic_fvector& fv) const {
        return !this->_equals(fv);
    }

    basic_fvector& operator <<(const basic_fvector& fv)
    {
        this->mv.resize(fv.size());
        this->_assign(fv);
        return *this;
    }

    basic_fvector drv(size_t nVarNum) const
    {
        basic_fvector vRet(this->size());
        vRet._drv(*this, nVarNum);
        return vRet;
    }

    basic_fvector& simp()
    {
        this->_simp();
        return *this;
    }

    basic_fvector operator + (const basic_fvector& fv) const throw(cvmexception)
    {
        this->_check_size(fv.size());
        basic_fvector vRet(this->size());
        vRet._add(*this, fv);
        return vRet;
    }

    basic_fvector operator - (const basic_fvector& fv) const throw(cvmexception)
    {
        this->_check_size(fv.size());
        basic_fvector vRet(this->size());
        vRet._subtract(*this, fv);
        return vRet;
    }

    basic_fvector& operator += (const basic_fvector& fv) throw(cvmexception)
    {
        this->_check_size(fv.size());
        this->_add(fv);
        return *this;
    }

    basic_fvector& operator -= (const basic_fvector& fv) throw(cvmexception)
    {
        this->_check_size(fv.size());
        this->_subtract(fv);
        return *this;
    }

    basic_fvector operator *(const BaseFunction& f) const
    {
        basic_fvector vRet(this->size());
        vRet.BaseFArray::_mult(*this, f);
        return vRet;
    }

    basic_fvector operator / (const BaseFunction& f) const
    {
        basic_fvector vRet(this->size());
        vRet._div(*this, f);
        return vRet;
    }

    basic_fvector operator *(const T& d) const
    {
        basic_fvector vRet(this->size());
        vRet.BaseFArray::_mult(*this, d);
        return vRet;
    }

    basic_fvector operator / (const T& d) const
    {
        basic_fvector vRet(this->size());
        vRet._div(*this, d);
        return vRet;
    }

    basic_fvector& operator *= (const BaseFunction& f)
    {
        this->BaseFArray::_mult(f);
        return *this;
    }

    basic_fvector& operator /= (const BaseFunction& f)
    {
        this->BaseFArray::_div(f);
        return *this;
    }

    basic_fvector& operator *= (const T& d)
    {
        this->BaseFArray::_mult(d);
        return *this;
    }

    basic_fvector& operator /= (const T& d)
    {
        this->BaseFArray::_div(d);
        return *this;
    }

    BaseFunction operator *(const basic_fvector& fv) const throw(cvmexception)
    {
        this->_check_size(fv.size());
        BaseFunction fRet;
        fRet._scalar_product(this->mv.size(), &this->mv.at(0), 1, &fv.mv.at(0), 1);
        return fRet;
    }

    basic_fvector operator *(const basic_fmatrix<T>& fm) const throw(cvmexception)
    {
        this->_check_size(fm.msize());
        basic_fvector vRet(fm.nsize());
        vRet._mult(*this, fm);
        return vRet;
    }

    basic_fvector& mult(const basic_fvector& fv, const basic_fmatrix<T>& fm) throw(cvmexception)
    {
        fv._check_size(fm.msize());
        this->_check_size(fm.nsize());
        return this->_mult(fv, fm);
    }

    basic_fvector& mult(const basic_fmatrix<T>& fm, const basic_fvector& fv) throw(cvmexception)
    {
        fv._check_size(fm.nsize());
        this->_check_size(fm.msize());
        return this->_mult(fm, fv);
    }

    basic_fmatrix<T> jacobian(size_t nfrom = 0, size_t vars = 0) const
    {
        const size_t m = this->size();
        if (m <= 0) {
            return basic_fmatrix<T>();
        }
        const size_t sz = this->mv.at(0).vars().size();
        _check_ge(CVM_INDEX_GE, nfrom, sz);

        const size_t n = vars > 0 ? vars : sz - nfrom;
        _check_gt(CVM_INDEX_GT, nfrom + n, sz);

        basic_fmatrix<T> fmj(m, n);
        this->_jacobi(fmj, nfrom);
        return fmj;
    }

    void jacobian(basic_fmatrix<T>& fmj, size_t nfrom = 0, size_t vars = 0) const
    {
        const size_t m = this->size();
        _check_ne(CVM_SIZESMISMATCH, m, fmj.msize());
        if (m > 0) {
            const size_t sz = this->mv.at(0).vars().size();
            _check_ge(CVM_INDEX_GE, nfrom, sz);
            const size_t n = vars > 0 ? vars : sz - nfrom;
            _check_gt(CVM_INDEX_GT, nfrom + n, sz);
            _check_ne(CVM_SIZESMISMATCH, n, fmj.nsize());
            this->_jacobi(fmj, nfrom);
        }
    }

    friend std::ostream& operator <<(std::ostream& os, const basic_fvector<T>& fv)
    {
        for (size_t i = 0; i < fv.size(); ++i) {
            os << fv.mv.at(i) << " ";
        }
        os << std::endl;
        return os;
    }

protected:
    // internal version, it doesn't verify anything
    basic_fvector& _mult(const basic_fvector& fv, const basic_fmatrix<T>& fm) throw(cvmexception)
    {
        for (size_t i = 0; i < fm.nsize(); ++i) {
            this->mv.at(i)._scalar_product(fv.size(), &fv.mv.at(0), 1, &fm.mv.at(fm.msize() * i), 1);
        }
        return *this;
    }
    // internal version, it doesn't verify anything
    basic_fvector& _mult(const basic_fmatrix<T>& fm, const basic_fvector& fv) throw(cvmexception)
    {
        for (size_t i = 0; i < fm.msize(); ++i) {
            this->mv.at(i)._scalar_product(fm.nsize(), &fm.mv.at(i), fm.msize(), &fv.mv.at(0), 1);
        }
        return *this;
    }
    // internal version, it doesn't verify anything
    // starts from nfrom-th variable (0-based)
    void _jacobi(basic_fmatrix<T>& fmj, size_t nfrom) const
    {
        for (size_t i = 0; i < fmj.msize(); ++i) {
            for (size_t j = 0; j < fmj.nsize(); ++j) {
                fmj.mv.at(fmj.msize() * j + i) = this->mv.at(i).drv(j + nfrom);
            }
        }
    }
};


class rfvector : public basic_fvector<treal>
{
protected:
    typedef basic_fvector<treal> BaseFVector; 

public:
    rfvector()
      : BaseFVector()
    {
    }

    explicit rfvector(size_t nSize)
      : BaseFVector(nSize)
    {
    }

    explicit rfvector(const string_array& saInput)
      : BaseFVector(saInput)
    {
    }

    rfvector(const string_array& saVars, const string_array& saBodies,
                  const string_array& saParameters, const string_array& saMeanings)
      : BaseFVector(saVars, saBodies, saParameters, saMeanings)
    {
    }

    rfvector(const BaseFVector& fv)
      : BaseFVector(fv)
    {
    }

    rfvector(BaseFVector&& fv)
      : BaseFVector(std::move(fv))
    {
    }

    rvector operator () () const {
        rvector ret((tint)this->size());
        this->value(ret);
        return ret;
    }

    rvector operator () (treal d) const {
        rvector ret((tint)this->size());
        this->value(d, ret);
        return ret;
    }

    rvector operator () (treal d1, treal d2) const {
        rvector ret((tint)this->size());
        this->value(d1, d2, ret);
        return ret;
    }

    rvector operator () (treal d1, treal d2, treal d3) const {
        rvector ret((tint)this->size());
        this->value(d1, d2, d3, ret);
        return ret;
    }

    rvector operator () (const treal* pd) const {
        rvector ret((tint)this->size());
        this->value(pd, ret);
        return ret;
    }
};


class cfvector : public basic_fvector<tcomplex>
{
protected:
    typedef basic_fvector<tcomplex> BaseFVector; 

public:
    cfvector()
      : BaseFVector()
    {
    }

    explicit cfvector(size_t nSize)
      : BaseFVector(nSize)
    {
    }

    explicit cfvector(const string_array& saInput)
      : BaseFVector(saInput)
    {
    }

    cfvector(const string_array& saVars, const string_array& saBodies,
                  const string_array& saParameters, const string_array& saMeanings)
      : BaseFVector(saVars, saBodies, saParameters, saMeanings)
    {
    }

    cfvector(const BaseFVector& fv)
      : BaseFVector(fv)
    {
    }

    cfvector(BaseFVector&& fv)
      : BaseFVector(std::move(fv))
    {
    }

    cvector operator () () const {
        cvector ret((tint)this->size());
        this->value(ret);
        return ret;
    }

    cvector operator () (tcomplex d) const {
        cvector ret((tint)this->size());
        this->value(d, ret);
        return ret;
    }

    cvector operator () (tcomplex d1, tcomplex d2) const {
        cvector ret((tint)this->size());
        this->value(d1, d2, ret);
        return ret;
    }

    cvector operator () (tcomplex d1, tcomplex d2, tcomplex d3) const {
        cvector ret((tint)this->size());
        this->value(d1, d2, d3, ret);
        return ret;
    }

    cvector operator () (const tcomplex* pd) const {
        cvector ret((tint)this->size());
        this->value(pd, ret);
        return ret;
    }

};


template<typename T>
class basic_fmatrix : public FArray<T>
{
    friend class basic_fvector<T>;  // _mult

protected:
    typedef FArray<T> BaseFArray; 
    typedef typename BaseFArray::BaseFunction BaseFunction; 

    size_t mM; 
    size_t mN; 

    void _check_sizes(const basic_fmatrix& fm) const throw(cvmexception)
    {
        _check_ne(CVM_SIZESMISMATCH, fm.msize(), this->msize());
        _check_ne(CVM_SIZESMISMATCH, fm.nsize(), this->nsize());
    }
    void _check_sizes() const throw(cvmexception)
    {
        static const size_t ZERO(0);
        _check_le(CVM_WRONGSIZE_LE, this->msize(), ZERO);
        _check_le(CVM_WRONGSIZE_LE, this->nsize(), ZERO);
        _check_ne(CVM_SIZESMISMATCH, mM * mN, this->size());
    }

public:
    basic_fmatrix()
      : BaseFArray(),
        mM(0),
        mN(0)
    {
    }

    basic_fmatrix(size_t m, size_t n)
      : BaseFArray(m * n),
        mM(m),
        mN(n)
    {
        this->_check_sizes();
    }

    basic_fmatrix(size_t m, size_t n, const string_array& saInput)
      : BaseFArray(saInput),
        mM(m),
        mN(n)
    {
        this->_check_sizes();
    }

    basic_fmatrix(size_t m, size_t n, const string_array& saVars, const string_array& saBodies,
                  const string_array& saParameters, const string_array& saMeanings)
      : BaseFArray(saVars, saBodies, saParameters, saMeanings), mM(m), mN(n)
    {
        this->_check_sizes();
    }

    basic_fmatrix(const basic_fmatrix& fm)
      : BaseFArray(fm), mM(fm.msize()), mN(fm.nsize())
    {
    }

    basic_fmatrix(basic_fmatrix&& fm)
      : BaseFArray(std::move(fm)), mM(fm.msize()), mN(fm.nsize())
    {
    }

    size_t msize() const {
        return mM;
    }

    size_t nsize() const {
        return mN;
    }

    // no destructor here because non-virtual std::vector::~vector does the job

    basic_fmatrix& operator = (const basic_fmatrix& fm) throw(cvmexception)
    {
        this->_check_sizes(fm);
        this->_assign(fm);
        return *this;
    }

    basic_fmatrix& operator = (basic_fmatrix&& fm) throw(cvmexception)
    {
        this->_check_sizes(fm);
        this->mp = std::move(fm.mv);
        return *this;
    }

    BaseFunction& at(size_t nRow, size_t nCol) {
        return this->mv.at(mM * nCol + nRow);
    }

    const BaseFunction& at(size_t nRow, size_t nCol) const {
        return this->mv.at(mM * nCol + nRow);
    }

    bool operator == (const basic_fmatrix& fm) const {
        return this->msize() == fm.msize() && this->nsize() == fm.nsize() && this->_equals(fm);
    }

    bool operator != (const basic_fmatrix& fm) const {
        return !this->operator == (fm);
    }

    basic_fmatrix& operator<<(const basic_fmatrix& fm) {
        this->mv.resize(fm.size());
        this->_assign(fm);
        this->mM = fm.msize();
        this->mN = fm.nsize();
        return *this;
    }

    basic_fmatrix drv(size_t nVarNum) const
    {
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet._drv(*this, nVarNum);
        return fmRet;
    }

    basic_fmatrix& simp()
    {
        this->_simp();
        return *this;
    }

    basic_fmatrix operator + (const basic_fmatrix& fm) const throw(cvmexception)
    {
        this->_check_sizes(fm);
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet._add(*this, fm);
        return fmRet;
    }

    basic_fmatrix operator - (const basic_fmatrix& fm) const throw(cvmexception)
    {
        this->_check_sizes(fm);
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet._subtract(*this, fm);
        return fmRet;
    }

    basic_fmatrix& operator += (const basic_fmatrix& fm) throw(cvmexception)
    {
        this->_check_sizes(fm);
        this->_add(fm);
        return *this;
    }

    basic_fmatrix& operator -= (const basic_fmatrix& fm) throw(cvmexception)
    {
        this->_check_sizes(fm);
        this->_subtract(fm);
        return *this;
    }

    basic_fmatrix operator *(const BaseFunction& f) const
    {
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet.BaseFArray::_mult(*this, f);
        return fmRet;
    }

    basic_fmatrix operator / (const BaseFunction& f) const
    {
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet._div(*this, f);
        return fmRet;
    }

    basic_fmatrix operator *(const T& d) const
    {
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet.BaseFArray::_mult(*this, d);
        return fmRet;
    }

    basic_fmatrix operator / (const T& d) const
    {
        basic_fmatrix fmRet(this->mM, this->mN);
        fmRet._div(*this, d);
        return fmRet;
    }

    basic_fmatrix& operator *= (const BaseFunction& f) {
        this->BaseFArray::_mult(f);
        return *this;
    }

    basic_fmatrix& operator /= (const BaseFunction& f)
    {
        this->BaseFArray::_div(f);
        return *this;
    }

    basic_fmatrix& operator *= (const T& d)
    {
        this->BaseFArray::_mult(d);
        return *this;
    }

    basic_fmatrix& operator /= (const T& d)
    {
        this->BaseFArray::_div(d);
        return *this;
    }

    basic_fvector<T> get_row(size_t nRow) const throw(cvmexception)
    {
        _check_ge(CVM_INDEX_GE, nRow, this->mM);
        basic_fvector<T> fvRet(this->mN);
        _copy(this->mN, &this->mv.at(nRow), this->mM, &fvRet.mv.at(0), 1);
        return fvRet;
    }

    basic_fvector<T> get_col(size_t nCol) const throw(cvmexception)
    {
        _check_ge(CVM_INDEX_GE, nCol, this->mN);
        basic_fvector<T> fvRet(this->mM);
        _copy(this->mM, &this->mv.at(this->mM * nCol), 1, &fvRet.mv.at(0), 1);
        return fvRet;
    }

    basic_fmatrix& set_row(size_t nRow, const basic_fvector<T>& fv) throw(cvmexception)
    {
        _check_ge(CVM_INDEX_GE, nRow, this->mM);
        _check_ne(CVM_SIZESMISMATCH, fv.size(), this->nsize());
        _copy(this->mN, &fv.mv.at(0), 1, &this->mv.at(nRow), this->mM);
        return *this;
    }

    basic_fmatrix& set_col(size_t nCol, const basic_fvector<T>& fv) throw(cvmexception)
    {
        _check_ge(CVM_INDEX_GE, nCol, this->nsize());
        _check_ne(CVM_SIZESMISMATCH, fv.size(), this->msize());
        _copy(this->mM, &fv.mv.at(0), 1, &this->mv.at(this->mM * nCol), 1);
        return *this;
    }

    basic_fvector<T> operator *(const basic_fvector<T>& fv) const throw(cvmexception)
    {
        _check_ne(CVM_SIZESMISMATCH, fv.size(), this->nsize());
        basic_fvector<T> vRet(this->msize());
        vRet._mult(*this, fv);
        return vRet;
    }

    basic_fmatrix operator *(const basic_fmatrix& fm) const throw(cvmexception)
    {
        _check_ne(CVM_SIZESMISMATCH, fm.msize(), this->nsize());
        basic_fmatrix mRet(this->msize(), fm.nsize());
        mRet._mult(*this, fm);
        return mRet;
    }

    basic_fmatrix& mult(const basic_fmatrix& fmA, const basic_fmatrix& fmB) throw(cvmexception)
    {
        _check_ne(CVM_SIZESMISMATCH, fmA.msize(), this->msize());
        _check_ne(CVM_SIZESMISMATCH, fmB.nsize(), this->nsize());
        _check_ne(CVM_SIZESMISMATCH, fmB.msize(), fmA.nsize());
        return this->_mult(fmA, fmB);
    }

    friend std::ostream& operator <<(std::ostream& os, const basic_fmatrix<T>& fm)
    {
        for (size_t i = 0; i < fm.msize(); ++i) {
            for (size_t j = 0; j < fm.nsize(); ++j) {
                os << fm.mv.at(fm.msize() * j + i) << CFUN_SSPACE;
            }
            os << std::endl;
        }
        return os;
    }

protected:
    // internal version, it doesn't verify anything
    basic_fmatrix& _mult(const basic_fmatrix& fmA, const basic_fmatrix& fmB) throw(cvmexception)
    {
        for (size_t i = 0; i < this->mM; ++i) {
            for (size_t j = 0; j < this->mN; ++j) {
                this->mv.at(this->mM * j + i)._scalar_product(fmA.nsize(), &fmA.mv.at(i), this->mM, &fmB.mv.at(fmB.msize() * j), 1);
            }
        }
        return *this;
    }
};


class rfmatrix : public basic_fmatrix<treal>
{
protected:
    typedef basic_fmatrix<treal> BaseFMatrix; 

public:
    rfmatrix()
      : BaseFMatrix()
    {
    }

    rfmatrix(size_t m, size_t n)
      : BaseFMatrix(m, n)
    {
    }

    rfmatrix(size_t m, size_t n, const string_array& saInput)
      : BaseFMatrix(m, n, saInput)
    {
    }

    rfmatrix(size_t m, size_t n, const string_array& saVars, const string_array& saBodies,
             const string_array& saParameters, const string_array& saMeanings)
      : BaseFMatrix(m, n, saVars, saBodies, saParameters, saMeanings)
    {
    }

    rfmatrix(const BaseFMatrix& fm)
      : BaseFMatrix(fm)
    {
    }

    rfmatrix(BaseFMatrix&& fm)
      : BaseFMatrix(std::move(fm))
    {
    }

    rmatrix operator () () const {
        rmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(ret);
        return ret;
    }

    rmatrix operator () (treal d) const {
        rmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d, ret);
        return ret;
    }

    rmatrix operator () (treal d1, treal d2) const {
        rmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d1, d2, ret);
        return ret;
    }

    rmatrix operator () (treal d1, treal d2, treal d3) const {
        rmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d1, d2, d3, ret);
        return ret;
    }

    rmatrix operator () (const treal* pd) const {
        rmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(pd, ret);
        return ret;
    }

};



class cfmatrix : public basic_fmatrix<tcomplex>
{
protected:
    typedef basic_fmatrix<tcomplex> BaseFMatrix; 

public:
    cfmatrix()
      : BaseFMatrix()
    {
    }

    cfmatrix(size_t m, size_t n)
      : BaseFMatrix(m, n)
    {
    }

    cfmatrix(size_t m, size_t n, const string_array& saInput)
      : BaseFMatrix(m, n, saInput)
    {
    }

    cfmatrix(size_t m, size_t n, const string_array& saVars, const string_array& saBodies,
             const string_array& saParameters, const string_array& saMeanings)
      : BaseFMatrix(m, n, saVars, saBodies, saParameters, saMeanings)
    {
    }

    cfmatrix(const BaseFMatrix& fm)
      : BaseFMatrix(fm)
    {
    }

    cfmatrix(BaseFMatrix&& fm)
      : BaseFMatrix(std::move(fm))
    {
    }

    cmatrix operator () () const {
        cmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(ret);
        return ret;
    }

    cmatrix operator () (tcomplex d) const {
        cmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d, ret);
        return ret;
    }

    cmatrix operator () (tcomplex d1, tcomplex d2) const {
        cmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d1, d2, ret);
        return ret;
    }

    cmatrix operator () (tcomplex d1, tcomplex d2, tcomplex d3) const {
        cmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(d1, d2, d3, ret);
        return ret;
    }

    cmatrix operator () (const tcomplex* pd) const {
        cmatrix ret((tint)this->msize(), (tint)this->nsize());
        this->value(pd, ret);
        return ret;
    }

};





// specializations:
template<typename T>
class Comparator<T,T>
{
public:
    static bool lt(const T& l, const T& r)
    {
        return l < r;
    }
    static bool le(const T& l, const T& r)
    {
        return l <= r;
    }
    static bool gt(const T& l, const T& r)
    {
        return l > r;
    }
    static bool ge(const T& l, const T& r)
    {
        return l >= r;
    }
    static bool eq(const T& l, const T& r)
    {
        return cvm::_abs(l - r) <= cvm::basic_cvmMachMin<T>();
    }
};

template<typename T>
class Comparator<std::complex<T>, std::complex<T> >
{
public:
    static bool lt(const std::complex<T>& l, const std::complex<T>& r)
    {
        return l.real() < r.real();
    }
    static bool le(const std::complex<T>& l, const std::complex<T>& r)
    {
        return l.real() <= r.real();
    }
    static bool gt(const std::complex<T>& l, const std::complex<T>& r)
    {
        return l.real() > r.real();
    }
    static bool ge(const std::complex<T>& l, const std::complex<T>& r)
    {
        return l.real() >= r.real();
    }
    static bool eq(const std::complex<T>& l, const std::complex<T>& r)
    {
        return cvm::_abs(l.real() - r.real()) <= cvm::basic_cvmMachMin<T>() &&
               cvm::_abs(l.imag() - r.imag()) <= cvm::basic_cvmMachMin<T>();
    }
};

template<typename T>
class Comparator<T,std::complex<T> >
{
public:
    static bool lt(const T& l, const std::complex<T>& r)
    {
        return l < r.real();
    }
    static bool le(const T& l, const std::complex<T>& r)
    {
        return l <= r.real();
    }
    static bool gt(const T& l, const std::complex<T>& r)
    {
        return l > r.real();
    }
    static bool ge(const T& l, const std::complex<T>& r)
    {
        return l >= r.real();
    }
    static bool eq(const T& l, const std::complex<T>& r)
    {
        return cvm::_abs(l - r.real()) <= cvm::basic_cvmMachMin<T>() &&
               cvm::_abs(r.imag()) <= cvm::basic_cvmMachMin<T>();
    }
};

template<typename T>
class Comparator<std::complex<T>, T>
{
public:
    static bool lt(const std::complex<T>& l, const T& r)
    {
        return l.real() < r;
    }
    static bool le(const std::complex<T>& l, const T& r)
    {
        return l.real() <= r;
    }
    static bool gt(const std::complex<T>& l, const T& r)
    {
        return l.real() > r;
    }
    static bool ge(const std::complex<T>& l, const T& r)
    {
        return l.real() >= r;
    }
    static bool eq(const std::complex<T>& l, const T& r)
    {
        return cvm::_abs(l.real() - r) <= cvm::basic_cvmMachMin<T>() &&
               cvm::_abs(l.imag()) <= cvm::basic_cvmMachMin<T>();
    }
};

// end-user classes
typedef basic_function<treal>    rfunction; 
typedef basic_function<tcomplex> cfunction; 

CVM_NAMESPACE_END

#endif  // _CFUN_H