1 /* 2 __________ 3 _____ __ __\______ \_____ _______ ______ ____ _______ 4 / \ | | \| ___/\__ \ \_ __ \/ ___/_/ __ \\_ __ \ 5 | Y Y \| | /| | / __ \_| | \/\___ \ \ ___/ | | \/ 6 |__|_| /|____/ |____| (____ /|__| /____ > \___ >|__| 7 \/ \/ \/ \/ 8 9 Copyright (C) 2013 Ingo Berg 10 11 Permission is hereby granted, free of charge, to any person obtaining a copy of this 12 software and associated documentation files (the "Software"), to deal in the Software 13 without restriction, including without limitation the rights to use, copy, modify, 14 merge, publish, distribute, sublicense, and/or sell copies of the Software, and to 15 permit persons to whom the Software is furnished to do so, subject to the following conditions: 16 17 The above copyright notice and this permission notice shall be included in all copies or 18 substantial portions of the Software. 19 20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT 21 NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 22 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 23 DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 24 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 25 */ 26 #include "muParser.h" 27 #include "muParserTemplateMagic.h" 28 29 //--- Standard includes ------------------------------------------------------------------------ 30 #include <cmath> 31 #include <algorithm> 32 #include <numeric> 33 34 /** \brief Pi (what else?). */ 35 #define PARSER_CONST_PI 3.141592653589793238462643 36 37 /** \brief The eulerian number. */ 38 #define PARSER_CONST_E 2.718281828459045235360287 39 40 using namespace std; 41 42 /** \file 43 \brief Implementation of the standard floating point parser. 44 */ 45 46 47 /** \brief Namespace for mathematical applications. */ 48 namespace mu 49 { 50 51 52 //--------------------------------------------------------------------------- 53 // Trigonometric function Sin(value_type v)54 value_type Parser::Sin(value_type v) { return MathImpl<value_type>::Sin(v); } Cos(value_type v)55 value_type Parser::Cos(value_type v) { return MathImpl<value_type>::Cos(v); } Tan(value_type v)56 value_type Parser::Tan(value_type v) { return MathImpl<value_type>::Tan(v); } ASin(value_type v)57 value_type Parser::ASin(value_type v) { return MathImpl<value_type>::ASin(v); } ACos(value_type v)58 value_type Parser::ACos(value_type v) { return MathImpl<value_type>::ACos(v); } ATan(value_type v)59 value_type Parser::ATan(value_type v) { return MathImpl<value_type>::ATan(v); } ATan2(value_type v1,value_type v2)60 value_type Parser::ATan2(value_type v1, value_type v2) { return MathImpl<value_type>::ATan2(v1, v2); } Sinh(value_type v)61 value_type Parser::Sinh(value_type v) { return MathImpl<value_type>::Sinh(v); } Cosh(value_type v)62 value_type Parser::Cosh(value_type v) { return MathImpl<value_type>::Cosh(v); } Tanh(value_type v)63 value_type Parser::Tanh(value_type v) { return MathImpl<value_type>::Tanh(v); } ASinh(value_type v)64 value_type Parser::ASinh(value_type v) { return MathImpl<value_type>::ASinh(v); } ACosh(value_type v)65 value_type Parser::ACosh(value_type v) { return MathImpl<value_type>::ACosh(v); } ATanh(value_type v)66 value_type Parser::ATanh(value_type v) { return MathImpl<value_type>::ATanh(v); } 67 68 //--------------------------------------------------------------------------- 69 // Logarithm functions 70 71 // Logarithm base 2 Log2(value_type v)72 value_type Parser::Log2(value_type v) 73 { 74 #ifdef MUP_MATH_EXCEPTIONS 75 if (v<=0) 76 throw ParserError(ecDOMAIN_ERROR, _T("Log2")); 77 #endif 78 79 return MathImpl<value_type>::Log2(v); 80 } 81 82 // Logarithm base 10 Log10(value_type v)83 value_type Parser::Log10(value_type v) 84 { 85 #ifdef MUP_MATH_EXCEPTIONS 86 if (v<=0) 87 throw ParserError(ecDOMAIN_ERROR, _T("Log10")); 88 #endif 89 90 return MathImpl<value_type>::Log10(v); 91 } 92 93 // Logarithm base e (natural logarithm) Ln(value_type v)94 value_type Parser::Ln(value_type v) 95 { 96 #ifdef MUP_MATH_EXCEPTIONS 97 if (v<=0) 98 throw ParserError(ecDOMAIN_ERROR, _T("Ln")); 99 #endif 100 101 return MathImpl<value_type>::Log(v); 102 } 103 104 //--------------------------------------------------------------------------- 105 // misc Exp(value_type v)106 value_type Parser::Exp(value_type v) { return MathImpl<value_type>::Exp(v); } Abs(value_type v)107 value_type Parser::Abs(value_type v) { return MathImpl<value_type>::Abs(v); } Sqrt(value_type v)108 value_type Parser::Sqrt(value_type v) 109 { 110 #ifdef MUP_MATH_EXCEPTIONS 111 if (v<0) 112 throw ParserError(ecDOMAIN_ERROR, _T("sqrt")); 113 #endif 114 115 return MathImpl<value_type>::Sqrt(v); 116 } Rint(value_type v)117 value_type Parser::Rint(value_type v) { return MathImpl<value_type>::Rint(v); } Sign(value_type v)118 value_type Parser::Sign(value_type v) { return MathImpl<value_type>::Sign(v); } 119 120 //--------------------------------------------------------------------------- 121 /** \brief Callback for the unary minus operator. 122 \param v The value to negate 123 \return -v 124 */ UnaryMinus(value_type v)125 value_type Parser::UnaryMinus(value_type v) 126 { 127 return -v; 128 } 129 130 //--------------------------------------------------------------------------- 131 /** \brief Callback for the unary minus operator. 132 \param v The value to negate 133 \return -v 134 */ UnaryPlus(value_type v)135 value_type Parser::UnaryPlus(value_type v) 136 { 137 return v; 138 } 139 140 //--------------------------------------------------------------------------- 141 /** \brief Callback for adding multiple values. 142 \param [in] a_afArg Vector with the function arguments 143 \param [in] a_iArgc The size of a_afArg 144 */ Sum(const value_type * a_afArg,int a_iArgc)145 value_type Parser::Sum(const value_type *a_afArg, int a_iArgc) 146 { 147 if (!a_iArgc) 148 throw exception_type(_T("too few arguments for function sum.")); 149 150 value_type fRes=0; 151 for (int i=0; i<a_iArgc; ++i) fRes += a_afArg[i]; 152 return fRes; 153 } 154 155 //--------------------------------------------------------------------------- 156 /** \brief Callback for averaging multiple values. 157 \param [in] a_afArg Vector with the function arguments 158 \param [in] a_iArgc The size of a_afArg 159 */ Avg(const value_type * a_afArg,int a_iArgc)160 value_type Parser::Avg(const value_type *a_afArg, int a_iArgc) 161 { 162 if (!a_iArgc) 163 throw exception_type(_T("too few arguments for function sum.")); 164 165 value_type fRes=0; 166 for (int i=0; i<a_iArgc; ++i) fRes += a_afArg[i]; 167 return fRes/(value_type)a_iArgc; 168 } 169 170 171 //--------------------------------------------------------------------------- 172 /** \brief Callback for determining the minimum value out of a vector. 173 \param [in] a_afArg Vector with the function arguments 174 \param [in] a_iArgc The size of a_afArg 175 */ Min(const value_type * a_afArg,int a_iArgc)176 value_type Parser::Min(const value_type *a_afArg, int a_iArgc) 177 { 178 if (!a_iArgc) 179 throw exception_type(_T("too few arguments for function min.")); 180 181 value_type fRes=a_afArg[0]; 182 for (int i=0; i<a_iArgc; ++i) 183 fRes = std::min(fRes, a_afArg[i]); 184 185 return fRes; 186 } 187 188 189 //--------------------------------------------------------------------------- 190 /** \brief Callback for determining the maximum value out of a vector. 191 \param [in] a_afArg Vector with the function arguments 192 \param [in] a_iArgc The size of a_afArg 193 */ Max(const value_type * a_afArg,int a_iArgc)194 value_type Parser::Max(const value_type *a_afArg, int a_iArgc) 195 { 196 if (!a_iArgc) 197 throw exception_type(_T("too few arguments for function min.")); 198 199 value_type fRes=a_afArg[0]; 200 for (int i=0; i<a_iArgc; ++i) fRes = std::max(fRes, a_afArg[i]); 201 202 return fRes; 203 } 204 205 206 207 //--------------------------------------------------------------------------- 208 /** \brief Default value recognition callback. 209 \param [in] a_szExpr Pointer to the expression 210 \param [in, out] a_iPos Pointer to an index storing the current position within the expression 211 \param [out] a_fVal Pointer where the value should be stored in case one is found. 212 \return 1 if a value was found 0 otherwise. 213 */ IsVal(const char_type * a_szExpr,int * a_iPos,value_type * a_fVal)214 int Parser::IsVal(const char_type* a_szExpr, int *a_iPos, value_type *a_fVal) 215 { 216 value_type fVal(0); 217 218 stringstream_type stream(a_szExpr); 219 std::locale loc (Parser::s_locale,new space_out); 220 stream.seekg(0); // todo: check if this really is necessary 221 stream.imbue(loc); 222 stream >> fVal; 223 stringstream_type::pos_type iEnd = stream.tellg(); // Position after reading 224 225 if (iEnd==(stringstream_type::pos_type)-1) 226 return 0; 227 228 *a_iPos += (int)iEnd; 229 *a_fVal = fVal; 230 return 1; 231 } 232 233 234 //--------------------------------------------------------------------------- 235 /** \brief Constructor. 236 237 Call ParserBase class constructor and trigger Function, Operator and Constant initialization. 238 */ Parser()239 Parser::Parser() 240 :ParserBase() 241 { 242 AddValIdent(IsVal); 243 244 InitCharSets(); 245 InitFun(); 246 InitConst(); 247 InitOprt(); 248 } 249 250 //--------------------------------------------------------------------------- 251 /** \brief Define the character sets. 252 \sa DefineNameChars, DefineOprtChars, DefineInfixOprtChars 253 254 This function is used for initializing the default character sets that define 255 the characters to be useable in function and variable names and operators. 256 */ InitCharSets()257 void Parser::InitCharSets() 258 { 259 DefineNameChars( _T("0123456789_abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ") ); 260 DefineOprtChars( _T("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ+-*^/?<>=#!$%&|~'_{}") ); 261 DefineInfixOprtChars( _T("/+-*^?<>=#!$%&|~'_") ); 262 } 263 264 //--------------------------------------------------------------------------- 265 /** \brief Initialize the default functions. */ InitFun()266 void Parser::InitFun() 267 { 268 if (mu::TypeInfo<mu::value_type>::IsInteger()) 269 { 270 // When setting MUP_BASETYPE to an integer type 271 // Place functions for dealing with integer values here 272 // ... 273 // ... 274 // ... 275 } 276 else 277 { 278 // trigonometric functions 279 DefineFun(_T("sin"), Sin); 280 DefineFun(_T("cos"), Cos); 281 DefineFun(_T("tan"), Tan); 282 // arcus functions 283 DefineFun(_T("asin"), ASin); 284 DefineFun(_T("acos"), ACos); 285 DefineFun(_T("atan"), ATan); 286 DefineFun(_T("atan2"), ATan2); 287 // hyperbolic functions 288 DefineFun(_T("sinh"), Sinh); 289 DefineFun(_T("cosh"), Cosh); 290 DefineFun(_T("tanh"), Tanh); 291 // arcus hyperbolic functions 292 DefineFun(_T("asinh"), ASinh); 293 DefineFun(_T("acosh"), ACosh); 294 DefineFun(_T("atanh"), ATanh); 295 // Logarithm functions 296 DefineFun(_T("log2"), Log2); 297 DefineFun(_T("log10"), Log10); 298 DefineFun(_T("log"), Ln); 299 DefineFun(_T("ln"), Ln); 300 // misc 301 DefineFun(_T("exp"), Exp); 302 DefineFun(_T("sqrt"), Sqrt); 303 DefineFun(_T("sign"), Sign); 304 DefineFun(_T("rint"), Rint); 305 DefineFun(_T("abs"), Abs); 306 // Functions with variable number of arguments 307 DefineFun(_T("sum"), Sum); 308 DefineFun(_T("avg"), Avg); 309 DefineFun(_T("min"), Min); 310 DefineFun(_T("max"), Max); 311 } 312 } 313 314 //--------------------------------------------------------------------------- 315 /** \brief Initialize constants. 316 317 By default the parser recognizes two constants. Pi ("pi") and the eulerian 318 number ("_e"). 319 */ InitConst()320 void Parser::InitConst() 321 { 322 DefineConst(_T("_pi"), (value_type)PARSER_CONST_PI); 323 DefineConst(_T("_e"), (value_type)PARSER_CONST_E); 324 } 325 326 //--------------------------------------------------------------------------- 327 /** \brief Initialize operators. 328 329 By default only the unary minus operator is added. 330 */ InitOprt()331 void Parser::InitOprt() 332 { 333 DefineInfixOprt(_T("-"), UnaryMinus); 334 DefineInfixOprt(_T("+"), UnaryPlus); 335 } 336 337 //--------------------------------------------------------------------------- OnDetectVar(string_type *,int &,int &)338 void Parser::OnDetectVar(string_type * /*pExpr*/, int & /*nStart*/, int & /*nEnd*/) 339 { 340 // this is just sample code to illustrate modifying variable names on the fly. 341 // I'm not sure anyone really needs such a feature... 342 /* 343 344 345 string sVar(pExpr->begin()+nStart, pExpr->begin()+nEnd); 346 string sRepl = std::string("_") + sVar + "_"; 347 348 int nOrigVarEnd = nEnd; 349 cout << "variable detected!\n"; 350 cout << " Expr: " << *pExpr << "\n"; 351 cout << " Start: " << nStart << "\n"; 352 cout << " End: " << nEnd << "\n"; 353 cout << " Var: \"" << sVar << "\"\n"; 354 cout << " Repl: \"" << sRepl << "\"\n"; 355 nEnd = nStart + sRepl.length(); 356 cout << " End: " << nEnd << "\n"; 357 pExpr->replace(pExpr->begin()+nStart, pExpr->begin()+nOrigVarEnd, sRepl); 358 cout << " New expr: " << *pExpr << "\n"; 359 */ 360 } 361 362 //--------------------------------------------------------------------------- 363 /** \brief Numerically differentiate with regard to a variable. 364 \param [in] a_Var Pointer to the differentiation variable. 365 \param [in] a_fPos Position at which the differentiation should take place. 366 \param [in] a_fEpsilon Epsilon used for the numerical differentiation. 367 368 Numerical differentiation uses a 5 point operator yielding a 4th order 369 formula. The default value for epsilon is 0.00074 which is 370 numeric_limits<double>::epsilon() ^ (1/5) as suggested in the muparser 371 forum: 372 373 http://sourceforge.net/forum/forum.php?thread_id=1994611&forum_id=462843 374 */ Diff(value_type * a_Var,value_type a_fPos,value_type a_fEpsilon) const375 value_type Parser::Diff(value_type *a_Var, 376 value_type a_fPos, 377 value_type a_fEpsilon) const 378 { 379 value_type fRes(0), 380 fBuf(*a_Var), 381 f[4] = {0,0,0,0}, 382 fEpsilon(a_fEpsilon); 383 384 // Backwards compatible calculation of epsilon in case the user doesn't provide 385 // his own epsilon 386 if (fEpsilon==0) 387 fEpsilon = (a_fPos==0) ? (value_type)1e-10 : (value_type)1e-7 * a_fPos; 388 389 *a_Var = a_fPos+2 * fEpsilon; f[0] = Eval(); 390 *a_Var = a_fPos+1 * fEpsilon; f[1] = Eval(); 391 *a_Var = a_fPos-1 * fEpsilon; f[2] = Eval(); 392 *a_Var = a_fPos-2 * fEpsilon; f[3] = Eval(); 393 *a_Var = fBuf; // restore variable 394 395 fRes = (-f[0] + 8*f[1] - 8*f[2] + f[3]) / (12*fEpsilon); 396 return fRes; 397 } 398 } // namespace mu 399