1 /* 2 __________ 3 _____ __ __\______ \_____ _______ ______ ____ _______ 4 / \ | | \| ___/\__ \ \_ __ \/ ___/_/ __ \\_ __ \ 5 | Y Y \| | /| | / __ \_| | \/\___ \ \ ___/ | | \/ 6 |__|_| /|____/ |____| (____ /|__| /____ > \___ >|__| 7 \/ \/ \/ \/ 8 9 Copyright (C) 2011 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 28 //--- Standard includes ------------------------------------------------------------------------ 29 #include <cmath> 30 #include <algorithm> 31 #include <numeric> 32 33 /** \brief Pi (what else?). */ 34 #define PARSER_CONST_PI 3.141592653589793238462643 35 36 /** \brief The eulerian number. */ 37 #define PARSER_CONST_E 2.718281828459045235360287 38 39 using namespace std; 40 41 /** \file 42 \brief Implementation of the standard floating point parser. 43 */ 44 45 46 /** \brief Namespace for mathematical applications. */ 47 namespace mu 48 { 49 //--------------------------------------------------------------------------- 50 // Trigonometric function Sin(value_type v)51 value_type Parser::Sin(value_type v) { return sin(v); } Cos(value_type v)52 value_type Parser::Cos(value_type v) { return cos(v); } Tan(value_type v)53 value_type Parser::Tan(value_type v) { return tan(v); } ASin(value_type v)54 value_type Parser::ASin(value_type v) { return asin(v); } ACos(value_type v)55 value_type Parser::ACos(value_type v) { return acos(v); } ATan(value_type v)56 value_type Parser::ATan(value_type v) { return atan(v); } Sinh(value_type v)57 value_type Parser::Sinh(value_type v) { return sinh(v); } Cosh(value_type v)58 value_type Parser::Cosh(value_type v) { return cosh(v); } Tanh(value_type v)59 value_type Parser::Tanh(value_type v) { return tanh(v); } ASinh(value_type v)60 value_type Parser::ASinh(value_type v) { return log(v + sqrt(v * v + 1)); } ACosh(value_type v)61 value_type Parser::ACosh(value_type v) { return log(v + sqrt(v * v - 1)); } ATanh(value_type v)62 value_type Parser::ATanh(value_type v) { return ((value_type)0.5 * log((1 + v) / (1 - v))); } 63 64 //--------------------------------------------------------------------------- 65 // Logarithm functions Log2(value_type v)66 value_type Parser::Log2(value_type v) { return log(v)/log((value_type)2); } // Logarithm base 2 Log10(value_type v)67 value_type Parser::Log10(value_type v) { return log10(v); } // Logarithm base 10 Ln(value_type v)68 value_type Parser::Ln(value_type v) { return log(v); } // Logarithm base e (natural logarithm) 69 70 //--------------------------------------------------------------------------- 71 // misc Exp(value_type v)72 value_type Parser::Exp(value_type v) { return exp(v); } Abs(value_type v)73 value_type Parser::Abs(value_type v) { return fabs(v); } Sqrt(value_type v)74 value_type Parser::Sqrt(value_type v) { return sqrt(v); } Rint(value_type v)75 value_type Parser::Rint(value_type v) { return floor(v + (value_type)0.5); } Sign(value_type v)76 value_type Parser::Sign(value_type v) { return (value_type)((v<0) ? -1 : (v>0) ? 1 : 0); } 77 78 //--------------------------------------------------------------------------- 79 /** \brief Callback for the unary minus operator. 80 \param v The value to negate 81 \return -v 82 */ UnaryMinus(value_type v)83 value_type Parser::UnaryMinus(value_type v) 84 { 85 return -v; 86 } 87 88 //--------------------------------------------------------------------------- 89 /** \brief Callback for adding multiple values. 90 \param [in] a_afArg Vector with the function arguments 91 \param [in] a_iArgc The size of a_afArg 92 */ Sum(const value_type * a_afArg,int a_iArgc)93 value_type Parser::Sum(const value_type *a_afArg, int a_iArgc) 94 { 95 if (!a_iArgc) 96 throw exception_type(_T("too few arguments for function sum.")); 97 98 value_type fRes=0; 99 for (int i=0; i<a_iArgc; ++i) fRes += a_afArg[i]; 100 return fRes; 101 } 102 103 //--------------------------------------------------------------------------- 104 /** \brief Callback for averaging multiple values. 105 \param [in] a_afArg Vector with the function arguments 106 \param [in] a_iArgc The size of a_afArg 107 */ Avg(const value_type * a_afArg,int a_iArgc)108 value_type Parser::Avg(const value_type *a_afArg, int a_iArgc) 109 { 110 if (!a_iArgc) 111 throw exception_type(_T("too few arguments for function sum.")); 112 113 value_type fRes=0; 114 for (int i=0; i<a_iArgc; ++i) fRes += a_afArg[i]; 115 return fRes/(value_type)a_iArgc; 116 } 117 118 119 //--------------------------------------------------------------------------- 120 /** \brief Callback for determining the minimum value out of a vector. 121 \param [in] a_afArg Vector with the function arguments 122 \param [in] a_iArgc The size of a_afArg 123 */ Min(const value_type * a_afArg,int a_iArgc)124 value_type Parser::Min(const value_type *a_afArg, int a_iArgc) 125 { 126 if (!a_iArgc) 127 throw exception_type(_T("too few arguments for function min.")); 128 129 value_type fRes=a_afArg[0]; 130 for (int i=0; i<a_iArgc; ++i) fRes = std::min(fRes, a_afArg[i]); 131 132 return fRes; 133 } 134 135 136 //--------------------------------------------------------------------------- 137 /** \brief Callback for determining the maximum value out of a vector. 138 \param [in] a_afArg Vector with the function arguments 139 \param [in] a_iArgc The size of a_afArg 140 */ Max(const value_type * a_afArg,int a_iArgc)141 value_type Parser::Max(const value_type *a_afArg, int a_iArgc) 142 { 143 if (!a_iArgc) 144 throw exception_type(_T("too few arguments for function min.")); 145 146 value_type fRes=a_afArg[0]; 147 for (int i=0; i<a_iArgc; ++i) fRes = std::max(fRes, a_afArg[i]); 148 149 return fRes; 150 } 151 152 153 //--------------------------------------------------------------------------- 154 /** \brief Default value recognition callback. 155 \param [in] a_szExpr Pointer to the expression 156 \param [in, out] a_iPos Pointer to an index storing the current position within the expression 157 \param [out] a_fVal Pointer where the value should be stored in case one is found. 158 \return 1 if a value was found 0 otherwise. 159 */ IsVal(const char_type * a_szExpr,int * a_iPos,value_type * a_fVal)160 int Parser::IsVal(const char_type* a_szExpr, int *a_iPos, value_type *a_fVal) 161 { 162 value_type fVal(0); 163 164 stringstream_type stream(a_szExpr); 165 stream.seekg(0); // todo: check if this really is necessary 166 stream.imbue(Parser::s_locale); 167 stream >> fVal; 168 stringstream_type::pos_type iEnd = stream.tellg(); // Position after reading 169 170 if (iEnd==(stringstream_type::pos_type)-1) 171 return 0; 172 173 *a_iPos += (int)iEnd; 174 *a_fVal = fVal; 175 return 1; 176 } 177 178 179 //--------------------------------------------------------------------------- 180 /** \brief Constructor. 181 182 Call ParserBase class constructor and trigger Function, Operator and Constant initialization. 183 */ Parser()184 Parser::Parser() 185 :ParserBase() 186 { 187 AddValIdent(IsVal); 188 189 InitCharSets(); 190 InitFun(); 191 InitConst(); 192 InitOprt(); 193 } 194 195 //--------------------------------------------------------------------------- 196 /** \brief Define the character sets. 197 \sa DefineNameChars, DefineOprtChars, DefineInfixOprtChars 198 199 This function is used for initializing the default character sets that define 200 the characters to be useable in function and variable names and operators. 201 */ InitCharSets()202 void Parser::InitCharSets() 203 { 204 DefineNameChars( _T("0123456789_abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ") ); 205 DefineOprtChars( _T("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ+-*^/?<>=#!$�%&|~'_{}") ); 206 DefineInfixOprtChars( _T("/+-*^?<>=#!$%&|~'_") ); 207 } 208 209 //--------------------------------------------------------------------------- 210 /** \brief Initialize the default functions. */ InitFun()211 void Parser::InitFun() 212 { 213 // trigonometric functions 214 DefineFun(_T("sin"), Sin); 215 DefineFun(_T("cos"), Cos); 216 DefineFun(_T("tan"), Tan); 217 // arcus functions 218 DefineFun(_T("asin"), ASin); 219 DefineFun(_T("acos"), ACos); 220 DefineFun(_T("atan"), ATan); 221 // hyperbolic functions 222 DefineFun(_T("sinh"), Sinh); 223 DefineFun(_T("cosh"), Cosh); 224 DefineFun(_T("tanh"), Tanh); 225 // arcus hyperbolic functions 226 DefineFun(_T("asinh"), ASinh); 227 DefineFun(_T("acosh"), ACosh); 228 DefineFun(_T("atanh"), ATanh); 229 // Logarithm functions 230 DefineFun(_T("log2"), Log2); 231 DefineFun(_T("log10"), Log10); 232 DefineFun(_T("log"), Log10); 233 DefineFun(_T("ln"), Ln); 234 // misc 235 DefineFun(_T("exp"), Exp); 236 DefineFun(_T("sqrt"), Sqrt); 237 DefineFun(_T("sign"), Sign); 238 DefineFun(_T("rint"), Rint); 239 DefineFun(_T("abs"), Abs); 240 // Functions with variable number of arguments 241 DefineFun(_T("sum"), Sum); 242 DefineFun(_T("avg"), Avg); 243 DefineFun(_T("min"), Min); 244 DefineFun(_T("max"), Max); 245 } 246 247 //--------------------------------------------------------------------------- 248 /** \brief Initialize constants. 249 250 By default the parser recognizes two constants. Pi ("pi") and the eulerian 251 number ("_e"). 252 */ InitConst()253 void Parser::InitConst() 254 { 255 DefineConst(_T("_pi"), (value_type)PARSER_CONST_PI); 256 DefineConst(_T("_e"), (value_type)PARSER_CONST_E); 257 } 258 259 //--------------------------------------------------------------------------- 260 /** \brief Initialize operators. 261 262 By default only the unary minus operator is added. 263 */ InitOprt()264 void Parser::InitOprt() 265 { 266 DefineInfixOprt(_T("-"), UnaryMinus); 267 } 268 269 //--------------------------------------------------------------------------- OnDetectVar(string_type *,int &,int &)270 void Parser::OnDetectVar(string_type * /*pExpr*/, int & /*nStart*/, int & /*nEnd*/) 271 { 272 // this is just sample code to illustrate modifying variable names on the fly. 273 // I'm not sure anyone really needs such a feature... 274 /* 275 276 277 string sVar(pExpr->begin()+nStart, pExpr->begin()+nEnd); 278 string sRepl = std::string("_") + sVar + "_"; 279 280 int nOrigVarEnd = nEnd; 281 cout << "variable detected!\n"; 282 cout << " Expr: " << *pExpr << "\n"; 283 cout << " Start: " << nStart << "\n"; 284 cout << " End: " << nEnd << "\n"; 285 cout << " Var: \"" << sVar << "\"\n"; 286 cout << " Repl: \"" << sRepl << "\"\n"; 287 nEnd = nStart + sRepl.length(); 288 cout << " End: " << nEnd << "\n"; 289 pExpr->replace(pExpr->begin()+nStart, pExpr->begin()+nOrigVarEnd, sRepl); 290 cout << " New expr: " << *pExpr << "\n"; 291 */ 292 } 293 294 //--------------------------------------------------------------------------- 295 /** \brief Numerically differentiate with regard to a variable. 296 \param [in] a_Var Pointer to the differentiation variable. 297 \param [in] a_fPos Position at which the differentiation should take place. 298 \param [in] a_fEpsilon Epsilon used for the numerical differentiation. 299 300 Numerical differentiation uses a 5 point operator yielding a 4th order 301 formula. The default value for epsilon is 0.00074 which is 302 numeric_limits<double>::epsilon() ^ (1/5) as suggested in the muparser 303 forum: 304 305 http://sourceforge.net/forum/forum.php?thread_id=1994611&forum_id=462843 306 */ Diff(value_type * a_Var,value_type a_fPos,value_type a_fEpsilon) const307 value_type Parser::Diff(value_type *a_Var, 308 value_type a_fPos, 309 value_type a_fEpsilon) const 310 { 311 value_type fRes(0), 312 fBuf(*a_Var), 313 f[4] = {0,0,0,0}, 314 fEpsilon(a_fEpsilon); 315 316 // Backwards compatible calculation of epsilon inc case the user doesnt provide 317 // his own epsilon 318 if (fEpsilon==0) 319 fEpsilon = (a_fPos==0) ? (value_type)1e-10 : (value_type)1e-7 * a_fPos; 320 321 *a_Var = a_fPos+2 * fEpsilon; f[0] = Eval(); 322 *a_Var = a_fPos+1 * fEpsilon; f[1] = Eval(); 323 *a_Var = a_fPos-1 * fEpsilon; f[2] = Eval(); 324 *a_Var = a_fPos-2 * fEpsilon; f[3] = Eval(); 325 *a_Var = fBuf; // restore variable 326 327 fRes = (-f[0] + 8*f[1] - 8*f[2] + f[3]) / (12*fEpsilon); 328 return fRes; 329 } 330 } // namespace mu 331