muParser.cpp

Go to the documentation of this file.

/*
                 __________
    _____   __ __\______   \_____  _______  ______  ____ _______
   /     \ |  |  \|     ___/\__  \ \_  __ \/  ___/_/ __ \\_  __ \
  |  Y Y  \|  |  /|    |     / __ \_|  | \/\___ \ \  ___/ |  | \/
  |__|_|  /|____/ |____|    (____  /|__|  /____  > \___  >|__|
        \/                       \/            \/      \/
 
  Copyright (C) 2012 Ingo Berg
 
  Permission is hereby granted, free of charge, to any person obtaining a copy of this
  software and associated documentation files (the "Software"), to deal in the Software
  without restriction, including without limitation the rights to use, copy, modify,
  merge, publish, distribute, sublicense, and/or sell copies of the Software, and to
  permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all copies or
  substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
  NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
  DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "muParser.h"
#include "muParserTemplateMagic.h"
 
//--- Standard includes ------------------------------------------------------------------------
#include <cmath>
#include <algorithm>
#include <numeric>
 
#define PARSER_CONST_PI  3.141592653589793238462643
 
#define PARSER_CONST_E   2.718281828459045235360287
 
using namespace std;
 
mu::value_type parser_Sum(const mu::value_type*, int);
mu::value_type parser_Avg(const mu::value_type*, int);
mu::value_type parser_Min(const mu::value_type*, int);
mu::value_type parser_Max(const mu::value_type*, int);
 
namespace mu
{
 
 
    //---------------------------------------------------------------------------
    // Trigonometric function
    value_type Parser::Sin(const value_type& v)
    {
        return MathImpl<value_type>::Sin(v);
    }
    value_type Parser::Cos(const value_type& v)
    {
        return MathImpl<value_type>::Cos(v);
    }
    value_type Parser::Tan(const value_type& v)
    {
        return MathImpl<value_type>::Tan(v);
    }
    value_type Parser::ASin(const value_type& v)
    {
        return MathImpl<value_type>::ASin(v);
    }
    value_type Parser::ACos(const value_type& v)
    {
        return MathImpl<value_type>::ACos(v);
    }
    value_type Parser::ATan(const value_type& v)
    {
        return MathImpl<value_type>::ATan(v);
    }
    value_type Parser::ATan2(const value_type& v1, const value_type& v2)
    {
        return MathImpl<double>::ATan2(v1.real(), v2.real());
    }
    value_type Parser::Sinh(const value_type& v)
    {
        return MathImpl<value_type>::Sinh(v);
    }
    value_type Parser::Cosh(const value_type& v)
    {
        return MathImpl<value_type>::Cosh(v);
    }
    value_type Parser::Tanh(const value_type& v)
    {
        return MathImpl<value_type>::Tanh(v);
    }
    value_type Parser::ASinh(const value_type& v)
    {
        return MathImpl<value_type>::ASinh(v);
    }
    value_type Parser::ACosh(const value_type& v)
    {
        return MathImpl<value_type>::ACosh(v);
    }
    value_type Parser::ATanh(const value_type& v)
    {
        return MathImpl<value_type>::ATanh(v);
    }
 
    //---------------------------------------------------------------------------
    // Logarithm functions
    value_type Parser::Log2(const value_type& v)
    {
        return MathImpl<value_type>::Log2(v);     // Logarithm base 2
    }
    value_type Parser::Log10(const value_type& v)
    {
        return MathImpl<value_type>::Log10(v);    // Logarithm base 10
    }
    value_type Parser::Ln(const value_type& v)
    {
        return MathImpl<value_type>::Log(v);      // Logarithm base e (natural logarithm)
    }
 
    //---------------------------------------------------------------------------
    //  misc
    value_type Parser::Exp(const value_type& v)
    {
        return MathImpl<value_type>::Exp(v);
    }
    value_type Parser::Abs(const value_type& v)
    {
        return (v.real() == 0.0 || v.imag() == 0.0) ? (std::abs(v.real()) + std::abs(v.imag())) : std::abs(v); //MathImpl<value_type>::Abs(v);
    }
    value_type Parser::Sqrt(const value_type& v)
    {
        return MathImpl<value_type>::Sqrt(v);
    }
    value_type Parser::Rint(const value_type& v)
    {
        return value_type(MathImpl<double>::Rint(v.real()), MathImpl<double>::Rint(v.imag()));
    }
    value_type Parser::Sign(const value_type& v)
    {
        return value_type(MathImpl<double>::Sign(v.real()), MathImpl<double>::Sign(v.imag()));
    }
 
    //---------------------------------------------------------------------------
    value_type Parser::UnaryMinus(const value_type& v)
    {
        return -v;
    }
 
    //---------------------------------------------------------------------------
    value_type Parser::Sum(const value_type* a_afArg, int a_iArgc)
    {
        if (!a_iArgc)
            throw exception_type(_nrT("too few arguments for function sum."));
        return parser_Sum(a_afArg, a_iArgc);
    }
 
    //---------------------------------------------------------------------------
    value_type Parser::Avg(const value_type* a_afArg, int a_iArgc)
    {
        if (!a_iArgc)
            throw exception_type(_nrT("too few arguments for function avg."));
        return parser_Avg(a_afArg, a_iArgc);
    }
 
 
    //---------------------------------------------------------------------------
    value_type Parser::Min(const value_type* a_afArg, int a_iArgc)
    {
        if (!a_iArgc)
            throw exception_type(_nrT("too few arguments for function min."));
        return parser_Min(a_afArg, a_iArgc);
    }
 
 
    //---------------------------------------------------------------------------
    value_type Parser::Max(const value_type* a_afArg, int a_iArgc)
    {
        if (!a_iArgc)
            throw exception_type(_nrT("too few arguments for function max."));
        return parser_Max(a_afArg, a_iArgc);
    }
 
 
    //---------------------------------------------------------------------------
    int Parser::IsVal(const char_type* a_szExpr, int* a_iPos, value_type* a_fVal)
    {
        value_type fVal(0);
 
        stringstream_type stream(a_szExpr);
        stream.seekg(0);        // todo:  check if this really is necessary
        stream.imbue(Parser::s_locale);
        stream >> fVal;
        stringstream_type::pos_type iEnd = stream.tellg(); // Position after reading
 
        if (iEnd == (stringstream_type::pos_type) - 1)
            return 0;
 
        *a_iPos += (int)iEnd;
        *a_fVal = fVal;
        return 1;
    }
 
 
    //---------------------------------------------------------------------------
    Parser::Parser()
        : ParserBase()
    {
        AddValIdent(IsVal);
 
        InitCharSets();
        InitFun();
        InitConst();
        InitOprt();
    }
 
    //---------------------------------------------------------------------------
    void Parser::InitCharSets()
    {
        DefineNameChars( _nrT("0123456789_abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ[]~\\") );
        DefineOprtChars( _nrT("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ+-*^/?<>=#!$%&|~'_{}") );
        DefineInfixOprtChars( _nrT("/+-*^?<>=#!$%&|~'_") );
    }
 
    //---------------------------------------------------------------------------
    void Parser::InitFun()
    {
        if (mu::TypeInfo<mu::value_type>::IsInteger())
        {
            // When setting MUP_BASETYPE to an integer type
            // Place functions for dealing with integer values here
            // ...
            // ...
            // ...
        }
        else
        {
            // trigonometric functions
            DefineFun("sin", Sin);
            DefineFun("cos", Cos);
            DefineFun("tan", Tan);
            // arcus functions
            DefineFun("asin", ASin);
            DefineFun("arcsin", ASin);
            DefineFun("acos", ACos);
            DefineFun("arccos", ACos);
            DefineFun("atan", ATan);
            DefineFun("arctan", ATan);
            DefineFun("atan2", ATan2);
            // hyperbolic functions
            DefineFun("sinh", Sinh);
            DefineFun("cosh", Cosh);
            DefineFun("tanh", Tanh);
            // arcus hyperbolic functions
            DefineFun("asinh", ASinh);
            DefineFun("arsinh", ASinh);
            DefineFun("acosh", ACosh);
            DefineFun("arcosh", ACosh);
            DefineFun("atanh", ATanh);
            DefineFun("artanh", ATanh);
            // Logarithm functions
            DefineFun("log2", Log2);
            DefineFun("log10", Log10);
            DefineFun("log", Log10);
            DefineFun("ln", Ln);
            // misc
            DefineFun("exp", Exp);
            DefineFun("sqrt", Sqrt);
            DefineFun("sign", Sign);
            DefineFun("rint", Rint);
            DefineFun("abs", Abs);
            // Functions with variable number of arguments
            DefineFun("sum", Sum);
            DefineFun("avg", Avg);
            DefineFun("min", Min);
            DefineFun("max", Max);
        }
    }
 
    //---------------------------------------------------------------------------
    void Parser::InitConst()
    {
        DefineConst("_pi", (value_type)PARSER_CONST_PI);
    }
 
    //---------------------------------------------------------------------------
    void Parser::InitOprt()
    {
        DefineInfixOprt("-", UnaryMinus);
    }
 
    //---------------------------------------------------------------------------
    void Parser::OnDetectVar(string_type* /*pExpr*/, int& /*nStart*/, int& /*nEnd*/)
    {
    }
    // this is just sample code to illustrate modifying variable names on the fly.
    // I'm not sure anyone really needs such a feature...
    /*
 
 
    string sVar(pExpr->begin()+nStart, pExpr->begin()+nEnd);
    string sRepl = std::string("_") + sVar + "_";
 
    int nOrigVarEnd = nEnd;
    cout << "variable detected!\n";
    cout << "  Expr: " << *pExpr << "\n";
    cout << "  Start: " << nStart << "\n";
    cout << "  End: " << nEnd << "\n";
    cout << "  Var: \"" << sVar << "\"\n";
    cout << "  Repl: \"" << sRepl << "\"\n";
    nEnd = nStart + sRepl.length();
    cout << "  End: " << nEnd << "\n";
    pExpr->replace(pExpr->begin()+nStart, pExpr->begin()+nOrigVarEnd, sRepl);
    cout << "  New expr: " << *pExpr << "\n";
    */
//}
 
    //---------------------------------------------------------------------------
    value_type Parser::Diff(value_type* a_Var,
                            value_type  a_fPos,
                            value_type  a_fEpsilon,
                            size_t order)
    {
        value_type fBuf(*a_Var),
                   fEpsilon(a_fEpsilon),
                   fRes;
        std::array<value_type, 5> f;
        std::array<double, 5> factors = {-2, -1, 0, 1, 2};
 
        // Backwards compatible calculation of epsilon inc case the user doesnt provide
        // his own epsilon
        if (fEpsilon == 0.0)
            fEpsilon = (a_fPos == 0.0) ? (value_type)1e-10 : (value_type)(1e-7*std::abs(a_fPos)*intPower(10, 2*(order-1)));
 
 
        for (size_t i = 0; i < f.size(); i++)
        {
            *a_Var = a_fPos + factors[i] * fEpsilon;
            f[i] = Eval();
        }
 
        // Reference: https://web.media.mit.edu/~crtaylor/calculator.html
        if (order == 1)
            fRes = ( f[0]  - 8.0 * f[1]              + 8.0 * f[3] - f[4]) / (12.0 * fEpsilon);
        else if (order == 2)
            fRes = (-f[0] + 16.0 * f[1] - 30.0*f[2] + 16.0 * f[3] - f[4]) / (12.0 * fEpsilon * fEpsilon);
        else if (order == 3)
            fRes = (-f[0]  + 2.0 * f[1]              - 2.0 * f[3] + f[4]) / (2.0 * fEpsilon * fEpsilon * fEpsilon);
        else
            return NAN;
 
        *a_Var = fBuf; // restore variable
        return fRes;
    }
} // namespace mu