From c26501501ea71d6166909914a3c5455eedf11e87 Mon Sep 17 00:00:00 2001 From: heltai Date: Sun, 7 Dec 2008 21:47:06 +0000 Subject: [PATCH] Updated function parser version git-svn-id: https://svn.dealii.org/trunk@17874 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/contrib/functionparser/example.cc | 50 + deal.II/contrib/functionparser/fparser.cc | 2351 +++-------------- deal.II/contrib/functionparser/fparser.h | 48 +- deal.II/contrib/functionparser/fparser.txt | 292 +- deal.II/contrib/functionparser/fpconfig.h | 67 + deal.II/contrib/functionparser/fpoptimizer.cc | 2094 +++++++++++++++ deal.II/contrib/functionparser/fptypes.h | 117 + deal.II/doc/news/changes.h | 7 + 8 files changed, 2892 insertions(+), 2134 deletions(-) create mode 100644 deal.II/contrib/functionparser/example.cc create mode 100644 deal.II/contrib/functionparser/fpconfig.h create mode 100644 deal.II/contrib/functionparser/fpoptimizer.cc create mode 100644 deal.II/contrib/functionparser/fptypes.h diff --git a/deal.II/contrib/functionparser/example.cc b/deal.II/contrib/functionparser/example.cc new file mode 100644 index 0000000000..608492e475 --- /dev/null +++ b/deal.II/contrib/functionparser/example.cc @@ -0,0 +1,50 @@ +// Simple example file for the function parser +// =========================================== + +/* When running the program, try for example with these values: + +f(x) = x^2 +min x: -5 +max x: 5 +step: 1 + +*/ + +#include "fparser.hh" + +#include +#include + +int main() +{ + std::string function; + double minx, maxx, step; + FunctionParser fparser; + + fparser.AddConstant("pi", 3.1415926535897932); + + while(true) + { + std::cout << "f(x) = "; + std::getline(std::cin, function); + int res = fparser.Parse(function, "x"); + if(res < 0) break; + + std::cout << std::string(res+7, ' ') << "^\n" + << fparser.ErrorMsg() << "\n\n"; + } + + std::cout << "min x: "; + std::cin >> minx; + std::cout << "max x: "; + std::cin >> maxx; + std::cout << "step: "; + std::cin >> step; + + double vals[1]; + for(vals[0] = minx; vals[0] <= maxx; vals[0] += step) + std::cout << "f(" << vals[0] << ") = " << fparser.Eval(vals) + << std::endl; + + return 0; +} diff --git a/deal.II/contrib/functionparser/fparser.cc b/deal.II/contrib/functionparser/fparser.cc index 456970953a..d71f7e5f3a 100644 --- a/deal.II/contrib/functionparser/fparser.cc +++ b/deal.II/contrib/functionparser/fparser.cc @@ -1,31 +1,12 @@ //=============================== -// Function parser v2.71 by Warp +// Function parser v2.83 by Warp //=============================== -// Comment out the following line if your compiler supports the (non-standard) -// asinh, acosh and atanh functions and you want them to be supported. If -// you are not sure, just leave it (those function will then not be supported). -#define NO_ASINH - - -// Uncomment the following line to disable the eval() function if it could -// be too dangerous in the target application: -//#define DISABLE_EVAL - - -// Comment this line out if you are not going to use the optimizer and want -// a slightly smaller library. The Optimize() method can still be called, -// but it will not do anything. -// If you are unsure, just leave it. It won't slow down the other parts of -// the library. -#ifndef NO_SUPPORT_OPTIMIZER -#define SUPPORT_OPTIMIZER -#endif - - -//============================================================================ - +#include "fpconfig.h" #include "fparser.h" +#include "fptypes.h" +using namespace FUNCTIONPARSERTYPES; +using namespace fparser; #include #include @@ -33,7 +14,12 @@ #include using namespace std; -using namespace fparser; + +#ifdef FP_USE_THREAD_SAFE_EVAL_WITH_ALLOCA +#ifndef FP_USE_THREAD_SAFE_EVAL +#define FP_USE_THREAD_SAFE_EVAL +#endif +#endif #ifndef M_PI #define M_PI 3.1415926535897932384626433832795 @@ -41,109 +27,6 @@ using namespace fparser; namespace { -// The functions must be in alphabetical order: - enum OPCODE - { - cAbs, cAcos, -#ifndef NO_ASINH - cAcosh, -#endif - cAsin, -#ifndef NO_ASINH - cAsinh, -#endif - cAtan, - cAtan2, -#ifndef NO_ASINH - cAtanh, -#endif - cCeil, cCos, cCosh, cCot, cCsc, -#ifndef DISABLE_EVAL - cEval, -#endif - cExp, cFloor, cIf, cInt, cLog, cLog10, cMax, cMin, - cSec, cSin, cSinh, cSqrt, cTan, cTanh, - -// These do not need any ordering: - cImmed, cJump, - cNeg, cAdd, cSub, cMul, cDiv, cMod, cPow, - cEqual, cLess, cGreater, cAnd, cOr, - - cDeg, cRad, - - cFCall, cPCall, - -#ifdef SUPPORT_OPTIMIZER - cVar, cDup, cInv, -#endif - - VarBegin - }; - - struct FuncDefinition - { - const char* name; - unsigned nameLength; - unsigned opcode; - unsigned params; - - // This is basically strcmp(), but taking 'nameLength' as string - // length (not ending '\0'): - bool operator<(const FuncDefinition& rhs) const - { - for(unsigned i = 0; i < nameLength; ++i) - { - if(i == rhs.nameLength) return false; - const char c1 = name[i], c2 = rhs.name[i]; - if(c1 < c2) return true; - if(c2 < c1) return false; - } - return nameLength < rhs.nameLength; - } - }; - - -// This list must be in alphabetical order: - const FuncDefinition Functions[]= - { - { "abs", 3, cAbs, 1 }, - { "acos", 4, cAcos, 1 }, -#ifndef NO_ASINH - { "acosh", 5, cAcosh, 1 }, -#endif - { "asin", 4, cAsin, 1 }, -#ifndef NO_ASINH - { "asinh", 5, cAsinh, 1 }, -#endif - { "atan", 4, cAtan, 1 }, - { "atan2", 5, cAtan2, 2 }, -#ifndef NO_ASINH - { "atanh", 5, cAtanh, 1 }, -#endif - { "ceil", 4, cCeil, 1 }, - { "cos", 3, cCos, 1 }, - { "cosh", 4, cCosh, 1 }, - { "cot", 3, cCot, 1 }, - { "csc", 3, cCsc, 1 }, -#ifndef DISABLE_EVAL - { "eval", 4, cEval, 0 }, -#endif - { "exp", 3, cExp, 1 }, - { "floor", 5, cFloor, 1 }, - { "if", 2, cIf, 0 }, - { "int", 3, cInt, 1 }, - { "log", 3, cLog, 1 }, - { "log10", 5, cLog10, 1 }, - { "max", 3, cMax, 2 }, - { "min", 3, cMin, 2 }, - { "sec", 3, cSec, 1 }, - { "sin", 3, cSin, 1 }, - { "sinh", 4, cSinh, 1 }, - { "sqrt", 4, cSqrt, 1 }, - { "tan", 3, cTan, 1 }, - { "tanh", 4, cTanh, 1 } - }; - const unsigned FUNC_AMOUNT = sizeof(Functions)/sizeof(Functions[0]); @@ -186,7 +69,7 @@ namespace //--------------------------------------------------------------------------- // Copy-on-write method //--------------------------------------------------------------------------- -inline void FunctionParser::copyOnWrite() +void FunctionParser::CopyOnWrite() { if(data->referenceCounter > 1) { @@ -204,7 +87,8 @@ inline void FunctionParser::copyOnWrite() //=========================================================================== FunctionParser::FunctionParser(): parseErrorType(FP_NO_ERROR), evalErrorType(0), - data(new Data) + data(new Data), + evalRecursionLevel(0) { data->referenceCounter = 1; } @@ -220,7 +104,8 @@ FunctionParser::~FunctionParser() FunctionParser::FunctionParser(const FunctionParser& cpy): parseErrorType(cpy.parseErrorType), evalErrorType(cpy.evalErrorType), - data(cpy.data) + data(cpy.data), + evalRecursionLevel(0) { ++(data->referenceCounter); } @@ -234,6 +119,7 @@ FunctionParser& FunctionParser::operator=(const FunctionParser& cpy) parseErrorType = cpy.parseErrorType; evalErrorType = cpy.evalErrorType; data = cpy.data; + evalRecursionLevel = cpy.evalRecursionLevel; ++(data->referenceCounter); } @@ -241,9 +127,14 @@ FunctionParser& FunctionParser::operator=(const FunctionParser& cpy) return *this; } +void FunctionParser::ForceDeepCopy() +{ + CopyOnWrite(); +} FunctionParser::Data::Data(): useDegreeConversion(false), + isOptimized(false), ByteCode(0), ByteCodeSize(0), Immed(0), ImmedSize(0), Stack(0), StackSize(0) @@ -259,7 +150,7 @@ FunctionParser::Data::~Data() // Makes a deep-copy of Data: FunctionParser::Data::Data(const Data& cpy): varAmount(cpy.varAmount), useDegreeConversion(cpy.useDegreeConversion), - Variables(cpy.Variables), Constants(cpy.Constants), + Variables(cpy.Variables), Constants(cpy.Constants), Units(cpy.Units), FuncPtrNames(cpy.FuncPtrNames), FuncPtrs(cpy.FuncPtrs), FuncParserNames(cpy.FuncParserNames), FuncParsers(cpy.FuncParsers), ByteCode(0), ByteCodeSize(cpy.ByteCodeSize), @@ -292,8 +183,8 @@ namespace "Empty parentheses", // 3 "Syntax error: Operator expected", // 4 "Not enough memory", // 5 - "An unexpected error ocurred. Please make a full bug report " - "to warp@iki.fi", // 6 + "An unexpected error occurred. Please make a full bug report " + "to the author", // 6 "Syntax error in parameter 'Vars' given to " "FunctionParser::Parse()", // 7 "Illegal number of parameters to function", // 8 @@ -303,76 +194,102 @@ namespace }; - // Parse variables - bool ParseVars(const string& Vars, map& dest) + enum VarNameType + { NOT_VALID_NAME, VALID_NEW_NAME, RESERVED_FUNCTION, + USER_DEF_CONST, USER_DEF_UNIT, + USER_DEF_FUNC_PTR, USER_DEF_FUNC_PARSER + }; +} + + +// Parse variables +bool FunctionParser::ParseVars(const string& Vars, map& dest) +{ + unsigned varNumber = VarBegin; + unsigned ind1 = 0, ind2; + + while(ind1 < Vars.size()) { - unsigned varNumber = VarBegin; - unsigned ind1 = 0, ind2; + if(!isalpha(Vars[ind1]) && Vars[ind1]!='_') return false; + for(ind2=ind1+1; ind2Constants) != data->Constants.end()) + return USER_DEF_CONST; - return true; + // Units are independent from the rest and thus the following check + // is actually not needed: + //if(FindConstant(n, data->Units) != data->Units.end()) + // return USER_DEF_UNIT; + + if(FindVariable(n, data->FuncParserNames) != data->FuncParserNames.end()) + return USER_DEF_FUNC_PARSER; + if(FindVariable(n, data->FuncPtrNames) != data->FuncPtrNames.end()) + return USER_DEF_FUNC_PTR; + + return VALID_NEW_NAME; } // Constants: bool FunctionParser::AddConstant(const string& name, double value) { - if(isValidName(name)) + int nameType = VarNameType(name); + if(nameType == VALID_NEW_NAME || nameType == USER_DEF_CONST) { - const char* n = name.c_str(); - if(FindVariable(n, data->FuncParserNames) != - data->FuncParserNames.end() || - FindVariable(n, data->FuncPtrNames) != - data->FuncPtrNames.end()) - return false; - - copyOnWrite(); - + CopyOnWrite(); data->Constants[name] = value; return true; } return false; } +// Units: +bool FunctionParser::AddUnit(const std::string& name, double value) +{ + if(!varNameHasOnlyValidChars(name)) return false; + CopyOnWrite(); + data->Units[name] = value; + return true; +} + // Function pointers bool FunctionParser::AddFunction(const std::string& name, FunctionPtr func, unsigned paramsAmount) { - if(paramsAmount == 0) return false; // Currently must be at least one - - if(isValidName(name)) + int nameType = VarNameType(name); + if(nameType == VALID_NEW_NAME) { - const char* n = name.c_str(); - if(FindVariable(n, data->FuncParserNames) != - data->FuncParserNames.end() || - FindConstant(n) != data->Constants.end()) - return false; - - copyOnWrite(); - + CopyOnWrite(); data->FuncPtrNames[name] = data->FuncPtrs.size(); data->FuncPtrs.push_back(Data::FuncPtrData(func, paramsAmount)); return true; @@ -380,30 +297,23 @@ bool FunctionParser::AddFunction(const std::string& name, return false; } -bool FunctionParser::checkRecursiveLinking(const FunctionParser* fp) const +bool FunctionParser::CheckRecursiveLinking(const FunctionParser* fp) const { if(fp == this) return true; for(unsigned i=0; idata->FuncParsers.size(); ++i) - if(checkRecursiveLinking(fp->data->FuncParsers[i])) return true; + if(CheckRecursiveLinking(fp->data->FuncParsers[i])) return true; return false; } bool FunctionParser::AddFunction(const std::string& name, FunctionParser& parser) { - if(parser.data->varAmount == 0) // Currently must be at least one - return false; - - if(isValidName(name)) + int nameType = VarNameType(name); + if(nameType == VALID_NEW_NAME) { - const char* n = name.c_str(); - if(FindVariable(n, data->FuncPtrNames) != data->FuncPtrNames.end() || - FindConstant(n) != data->Constants.end()) - return false; + if(CheckRecursiveLinking(&parser)) return false; - if(checkRecursiveLinking(&parser)) return false; - - copyOnWrite(); + CopyOnWrite(); data->FuncParserNames[name] = data->FuncParsers.size(); data->FuncParsers.push_back(&parser); @@ -420,7 +330,7 @@ int FunctionParser::Parse(const std::string& Function, const std::string& Vars, bool useDegrees) { - copyOnWrite(); + CopyOnWrite(); data->Variables.clear(); @@ -436,10 +346,11 @@ int FunctionParser::Parse(const std::string& Function, parseErrorType = FP_NO_ERROR; int Result = CheckSyntax(Func); - if(Result>=0) return Result; + if(Result >= 0) return Result; data->useDegreeConversion = useDegrees; - if(!Compile(Func)) return Function.size(); + Result = Compile(Func); + if(Result >= 0) return Result; data->Variables.clear(); @@ -449,10 +360,27 @@ int FunctionParser::Parse(const std::string& Function, namespace { + const char* const fpOperators[] = + { + "+", "-", "*", "/", "%", "^", + "=", "!=", "<=", "<", ">=", ">", "&", "|", + 0 + }; + // Is given char an operator? - inline bool IsOperator(int c) + // (Returns 0 if not, else the size of the operator) + inline int IsOperator(const char* F) { - return strchr("+-*/%^=<>&|",c)!=NULL; + for(unsigned opInd = 0; fpOperators[opInd]; ++opInd) + { + const char* op = fpOperators[opInd]; + for(unsigned n = 0; F[n] == *op; ++n) + { + ++op; + if(*op == 0) return op-fpOperators[opInd]; + } + } + return 0; } // skip whitespace @@ -481,19 +409,36 @@ FunctionParser::FindVariable(const char* F, const Data::VarMap_t& vars) const } inline FunctionParser::Data::ConstMap_t::const_iterator -FunctionParser::FindConstant(const char* F) const +FunctionParser::FindConstant(const char* F, + const Data::ConstMap_t& consts) const { - if(data->Constants.size()) + if(consts.size()) { unsigned ind = 0; while(isalnum(F[ind]) || F[ind] == '_') ++ind; if(ind) { string name(F, ind); - return data->Constants.find(name); + return consts.find(name); + } + } + return consts.end(); +} + +int FunctionParser::CheckForUnit(const char* Function, int Ind) const +{ + int c = Function[Ind]; + if(isalpha(c) || c == '_') + { + Data::ConstMap_t::const_iterator uIter = + FindConstant(&Function[Ind], data->Units); + if(uIter != data->Units.end()) + { + Ind += uIter->first.size(); + sws(Function, Ind); } } - return data->Constants.end(); + return Ind; } //--------------------------------------------------------------------------- @@ -518,8 +463,8 @@ int FunctionParser::CheckSyntax(const char* Function) // Check for valid operand (must appear) - // Check for leading - - if(c=='-') { sws(Function, ++Ind); c=Function[Ind]; } + // Check for leading - or ! + if(c=='-' || c=='!') { sws(Function, ++Ind); c=Function[Ind]; } if(c==0) { parseErrorType=PREMATURE_EOS; return Ind; } // Check for math function @@ -557,6 +502,18 @@ int FunctionParser::CheckSyntax(const char* Function) sws(Function, Ind); c = Function[Ind]; if(c!='(') { parseErrorType=EXPECT_PARENTH_FUNC; return Ind; } + + int Ind2 = Ind+1; + sws(Function, Ind2); + if(Function[Ind2] == ')') + { + Ind = Ind2+1; + sws(Function, Ind); + c = Function[Ind]; + // Ugly, but other methods would just be uglier... + goto CheckOperator; + } + functionParenthDepth.push_back(ParenthCnt+1); } @@ -587,7 +544,7 @@ int FunctionParser::CheckSyntax(const char* Function) { // Check for constant Data::ConstMap_t::const_iterator cIter = - FindConstant(&Function[Ind]); + FindConstant(&Function[Ind], Constants); if(cIter != Constants.end()) Ind += cIter->first.size(); else @@ -597,6 +554,10 @@ int FunctionParser::CheckSyntax(const char* Function) c = Function[Ind]; } + // Check for unit + Ind = CheckForUnit(Function, Ind); + c = Function[Ind]; + // Check for closing parenthesis while(c==')') { @@ -608,20 +569,30 @@ int FunctionParser::CheckSyntax(const char* Function) c=Function[Ind]; } -// If we get here, we have a legal operand and now a legal operator or + CheckOperator: +// If we get here, we have a legal operand and now a legal unit, operator or // end of string must follow + // Check for unit + Ind = CheckForUnit(Function, Ind); + c = Function[Ind]; + // Check for EOS if(c==0) break; // The only way to end the checking loop without error + // Check for operator - if(!IsOperator(c) && - (c != ',' || functionParenthDepth.empty() || - functionParenthDepth.back() != ParenthCnt)) + int opSize = 0; + if(c == ',' && !functionParenthDepth.empty() && + functionParenthDepth.back() == ParenthCnt) + opSize = 1; + else + opSize = IsOperator(Function+Ind); + if(opSize == 0) { parseErrorType=EXPECT_OPERATOR; return Ind; } // If we get here, we have an operand and an operator; the next loop will // check for another operand (must appear) - ++Ind; + Ind += opSize; } // while // Check that all opened parentheses are also closed @@ -635,11 +606,12 @@ int FunctionParser::CheckSyntax(const char* Function) // Compile function string to bytecode // ----------------------------------- -bool FunctionParser::Compile(const char* Function) +int FunctionParser::Compile(const char* Function) { if(data->ByteCode) { delete[] data->ByteCode; data->ByteCode=0; } if(data->Immed) { delete[] data->Immed; data->Immed=0; } if(data->Stack) { delete[] data->Stack; data->Stack=0; } + data->isOptimized = false; vector byteCode; byteCode.reserve(1024); tempByteCode = &byteCode; @@ -649,8 +621,8 @@ bool FunctionParser::Compile(const char* Function) data->StackSize = StackPtr = 0; - CompileExpression(Function, 0); - if(parseErrorType != FP_NO_ERROR) return false; + int ind = CompileExpression(Function, 0); + if(parseErrorType != FP_NO_ERROR) return ind; data->ByteCodeSize = byteCode.size(); data->ImmedSize = immed.size(); @@ -667,10 +639,12 @@ bool FunctionParser::Compile(const char* Function) memcpy(data->Immed, &immed[0], sizeof(double)*data->ImmedSize); } +#ifndef FP_USE_THREAD_SAFE_EVAL if(data->StackSize) data->Stack = new double[data->StackSize]; +#endif - return true; + return -1; } @@ -765,13 +739,22 @@ int FunctionParser::CompileIf(const char* F, int ind) int FunctionParser::CompileFunctionParams(const char* F, int ind, unsigned requiredParams) { - unsigned curStackPtr = StackPtr; - int ind2 = CompileExpression(F, ind); + int ind2 = ind; + if(requiredParams > 0) + { + unsigned curStackPtr = StackPtr; + ind2 = CompileExpression(F, ind); + + if(StackPtr != curStackPtr+requiredParams) + { parseErrorType=ILL_PARAMS_AMOUNT; return ind; } - if(StackPtr != curStackPtr+requiredParams) - { parseErrorType=ILL_PARAMS_AMOUNT; return ind; } + StackPtr -= requiredParams - 1; + } + else + { + incStackPtr(); + } - StackPtr -= requiredParams - 1; sws(F, ind2); return ind2+1; // F[ind2] is ')' } @@ -833,7 +816,8 @@ int FunctionParser::CompileElement(const char* F, int ind) return ind + vIter->first.size(); } - Data::ConstMap_t::const_iterator cIter = FindConstant(F+ind); + Data::ConstMap_t::const_iterator cIter = + FindConstant(F+ind, data->Constants); if(cIter != data->Constants.end()) // is constant { AddImmediate(cIter->second); @@ -881,15 +865,38 @@ int FunctionParser::CompileElement(const char* F, int ind) return ind; } +// Compiles a unit factor possibly appearing after an element +int FunctionParser::CompilePossibleUnit(const char* F, int ind) +{ + if(isalpha(F[ind]) || F[ind] == '_') + { + // If the syntax checker is bug-free, this should always work: + Data::ConstMap_t::const_iterator uIter = + FindConstant(F+ind, data->Units); + + AddImmediate(uIter->second); + AddCompiledByte(cImmed); + incStackPtr(); + AddCompiledByte(cMul); + --StackPtr; + ind += uIter->first.size(); + sws(F, ind); + } + return ind; +} + // Compiles '^' int FunctionParser::CompilePow(const char* F, int ind) { int ind2 = CompileElement(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); + ind2 = CompilePossibleUnit(F, ind2); while(F[ind2] == '^') { ind2 = CompileUnaryMinus(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); AddCompiledByte(cPow); --StackPtr; @@ -902,28 +909,30 @@ int FunctionParser::CompilePow(const char* F, int ind) int FunctionParser::CompileUnaryMinus(const char* F, int ind) { sws(F, ind); - if(F[ind] == '-') + if(F[ind] == '-' || F[ind] == '!') { int ind2 = ind+1; sws(F, ind2); ind2 = CompilePow(F, ind2); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); // if we are negating a constant, negate the constant itself: - if(tempByteCode->back() == cImmed) + if(F[ind] == '-' && tempByteCode->back() == cImmed) tempImmed->back() = -tempImmed->back(); // if we are negating a negation, we can remove both: - else if(tempByteCode->back() == cNeg) + else if((F[ind] == '-' && tempByteCode->back() == cNeg)) tempByteCode->pop_back(); else - AddCompiledByte(cNeg); + AddCompiledByte(F[ind] == '-' ? cNeg : cNot); return ind2; } int ind2 = CompilePow(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); return ind2; } @@ -932,12 +941,14 @@ int FunctionParser::CompileUnaryMinus(const char* F, int ind) int FunctionParser::CompileMult(const char* F, int ind) { int ind2 = CompileUnaryMinus(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); char op; while((op = F[ind2]) == '*' || op == '/' || op == '%') { ind2 = CompileUnaryMinus(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); switch(op) { @@ -955,12 +966,14 @@ int FunctionParser::CompileMult(const char* F, int ind) int FunctionParser::CompileAddition(const char* F, int ind) { int ind2 = CompileMult(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); char op; while((op = F[ind2]) == '+' || op == '-') { ind2 = CompileMult(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); AddCompiledByte(op=='+' ? cAdd : cSub); --StackPtr; @@ -973,18 +986,26 @@ int FunctionParser::CompileAddition(const char* F, int ind) int FunctionParser::CompileComparison(const char* F, int ind) { int ind2 = CompileAddition(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); char op; - while((op = F[ind2]) == '=' || op == '<' || op == '>') + while((op = F[ind2]) == '=' || op == '<' || op == '>' || op == '!') { - ind2 = CompileAddition(F, ind2+1); + int opSize = (F[ind2+1] == '=' ? 2 : 1); + ind2 = CompileAddition(F, ind2+opSize); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); switch(op) { - case '=': AddCompiledByte(cEqual); break; - case '<': AddCompiledByte(cLess); break; - case '>': AddCompiledByte(cGreater); break; + case '=': + AddCompiledByte(cEqual); break; + case '<': + AddCompiledByte(opSize == 1 ? cLess : cLessOrEq); break; + case '>': + AddCompiledByte(opSize == 1 ? cGreater : cGreaterOrEq); break; + case '!': + AddCompiledByte(cNEqual); break; } --StackPtr; } @@ -996,11 +1017,13 @@ int FunctionParser::CompileComparison(const char* F, int ind) int FunctionParser::CompileAnd(const char* F, int ind) { int ind2 = CompileComparison(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); while(F[ind2] == '&') { ind2 = CompileComparison(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); AddCompiledByte(cAnd); --StackPtr; @@ -1013,11 +1036,13 @@ int FunctionParser::CompileAnd(const char* F, int ind) int FunctionParser::CompileOr(const char* F, int ind) { int ind2 = CompileAnd(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); while(F[ind2] == '|') { ind2 = CompileAnd(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); AddCompiledByte(cOr); --StackPtr; @@ -1030,6 +1055,7 @@ int FunctionParser::CompileOr(const char* F, int ind) int FunctionParser::CompileExpression(const char* F, int ind, bool stopAtComma) { int ind2 = CompileOr(F, ind); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); if(stopAtComma) return ind2; @@ -1037,6 +1063,7 @@ int FunctionParser::CompileExpression(const char* F, int ind, bool stopAtComma) while(F[ind2] == ',') { ind2 = CompileOr(F, ind2+1); + if(parseErrorType != FP_NO_ERROR) return ind2; sws(F, ind2); } @@ -1087,11 +1114,20 @@ double FunctionParser::Eval(const double* Vars) { const unsigned* const ByteCode = data->ByteCode; const double* const Immed = data->Immed; - double* const Stack = data->Stack; const unsigned ByteCodeSize = data->ByteCodeSize; unsigned IP, DP=0; int SP=-1; +#ifdef FP_USE_THREAD_SAFE_EVAL +#ifdef FP_USE_THREAD_SAFE_EVAL_WITH_ALLOCA + double* const Stack = (double*)alloca(data->StackSize*sizeof(double)); +#else + std::vector Stack(data->StackSize); +#endif +#else + double* const Stack = data->Stack; +#endif + for(IP=0; IPStack = new double[data->StackSize]; - double retVal = Eval(&Stack[SP-data->varAmount+1]); - delete[] data->Stack; - data->Stack = Stack; + double retVal = 0; + if(evalRecursionLevel == EVAL_MAX_REC_LEVEL) + { + evalErrorType = 5; + } + else + { +#ifndef FP_USE_THREAD_SAFE_EVAL + data->Stack = new double[data->StackSize]; +#endif + ++evalRecursionLevel; + retVal = Eval(&Stack[SP-data->varAmount+1]); + --evalRecursionLevel; +#ifndef FP_USE_THREAD_SAFE_EVAL + delete[] data->Stack; + data->Stack = Stack; +#endif + } SP -= data->varAmount-1; Stack[SP] = retVal; break; @@ -1204,12 +1254,37 @@ double FunctionParser::Eval(const double* Vars) case cPow: Stack[SP-1] = pow(Stack[SP-1], Stack[SP]); --SP; break; +#ifdef FP_EPSILON + case cEqual: Stack[SP-1] = + (fabs(Stack[SP-1]-Stack[SP]) <= FP_EPSILON); + --SP; break; + case cNEqual: Stack[SP-1] = + (fabs(Stack[SP-1] - Stack[SP]) >= FP_EPSILON); + --SP; break; + case cLess: Stack[SP-1] = (Stack[SP-1] < Stack[SP]-FP_EPSILON); + --SP; break; + case cLessOrEq: Stack[SP-1] = (Stack[SP-1] <= Stack[SP]+FP_EPSILON); + --SP; break; + case cGreater: Stack[SP-1] = (Stack[SP-1]-FP_EPSILON > Stack[SP]); + --SP; break; + case cGreaterOrEq: Stack[SP-1] = + (Stack[SP-1]+FP_EPSILON >= Stack[SP]); + --SP; break; +#else case cEqual: Stack[SP-1] = (Stack[SP-1] == Stack[SP]); --SP; break; + case cNEqual: Stack[SP-1] = (Stack[SP-1] != Stack[SP]); + --SP; break; case cLess: Stack[SP-1] = (Stack[SP-1] < Stack[SP]); --SP; break; + case cLessOrEq: Stack[SP-1] = (Stack[SP-1] <= Stack[SP]); + --SP; break; case cGreater: Stack[SP-1] = (Stack[SP-1] > Stack[SP]); --SP; break; + case cGreaterOrEq: Stack[SP-1] = (Stack[SP-1] >= Stack[SP]); + --SP; break; +#endif + case cAnd: Stack[SP-1] = (doubleToInt(Stack[SP-1]) && doubleToInt(Stack[SP])); @@ -1218,6 +1293,7 @@ double FunctionParser::Eval(const double* Vars) (doubleToInt(Stack[SP-1]) || doubleToInt(Stack[SP])); --SP; break; + case cNot: Stack[SP] = !doubleToInt(Stack[SP]); break; // Degrees-radians conversion: case cDeg: Stack[SP] = RadiansToDegrees(Stack[SP]); break; @@ -1230,7 +1306,7 @@ double FunctionParser::Eval(const double* Vars) unsigned params = data->FuncPtrs[index].params; double retVal = data->FuncPtrs[index].ptr(&Stack[SP-params+1]); - SP -= params-1; + SP -= int(params)-1; Stack[SP] = retVal; break; } @@ -1241,7 +1317,7 @@ double FunctionParser::Eval(const double* Vars) unsigned params = data->FuncParsers[index]->data->varAmount; double retVal = data->FuncParsers[index]->Eval(&Stack[SP-params+1]); - SP -= params-1; + SP -= int(params)-1; Stack[SP] = retVal; break; } @@ -1268,17 +1344,20 @@ double FunctionParser::Eval(const double* Vars) #ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT +#include namespace { inline void printHex(std::ostream& dest, unsigned n) { - dest.width(8); dest.fill('0'); hex(dest); //uppercase(dest); + dest.width(8); dest.fill('0'); std::hex(dest); //uppercase(dest); dest << n; } } void FunctionParser::PrintByteCode(std::ostream& dest) const { + dest << "Size of stack: " << data->StackSize << "\n"; + const unsigned* const ByteCode = data->ByteCode; const double* const Immed = data->Immed; @@ -1315,7 +1394,8 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const Data::VarMap_t::const_iterator iter = data->FuncPtrNames.begin(); while(iter->second != index) ++iter; - dest << "call\t" << iter->first << endl; + dest << "fcall\t" << iter->first + << " (" << data->FuncPtrs[index].params << ")" << endl; break; } @@ -1325,7 +1405,9 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const Data::VarMap_t::const_iterator iter = data->FuncParserNames.begin(); while(iter->second != index) ++iter; - dest << "call\t" << iter->first << endl; + dest << "pcall\t" << iter->first + << " (" << data->FuncParsers[index]->data->varAmount + << ")" << endl; break; } @@ -1333,6 +1415,7 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const if(opcode < VarBegin) { string n; + unsigned params = 1; switch(opcode) { case cNeg: n = "neg"; break; @@ -1343,10 +1426,14 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const case cMod: n = "mod"; break; case cPow: n = "pow"; break; case cEqual: n = "eq"; break; + case cNEqual: n = "neq"; break; case cLess: n = "lt"; break; + case cLessOrEq: n = "le"; break; case cGreater: n = "gt"; break; + case cGreaterOrEq: n = "ge"; break; case cAnd: n = "and"; break; case cOr: n = "or"; break; + case cNot: n = "not"; break; case cDeg: n = "deg"; break; case cRad: n = "rad"; break; @@ -1360,9 +1447,13 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const case cInv: n = "inv"; break; #endif - default: n = Functions[opcode-cAbs].name; + default: + n = Functions[opcode-cAbs].name; + params = Functions[opcode-cAbs].params; } - dest << n << endl; + dest << n; + if(params != 1) dest << " (" << params << ")"; + dest << endl; } else { @@ -1374,1795 +1465,9 @@ void FunctionParser::PrintByteCode(std::ostream& dest) const #endif -//======================================================================== -// Optimization code was contributed by Bisqwit (http://iki.fi/bisqwit/) -//======================================================================== -#ifdef SUPPORT_OPTIMIZER - -#include -#include - -#define CONSTANT_E 2.71828182845904509080 // exp(1) -#define CONSTANT_PI M_PI // atan2(0,-1) -#define CONSTANT_L10 2.30258509299404590109 // log(10) -#define CONSTANT_L10I 0.43429448190325176116 // 1/log(10) -#define CONSTANT_L10E CONSTANT_L10I // log10(e) -#define CONSTANT_L10EI CONSTANT_L10 // 1/log10(e) -#define CONSTANT_DR (180.0 / M_PI) // 180/pi -#define CONSTANT_RD (M_PI / 180.0) // pi/180 - -namespace { -class compres -{ - // states: 0=false, 1=true, 2=unknown -public: - compres(bool b) : state(b) {} - compres(char v) : state(v) {} - // is it? - operator bool() const { return state != 0; } -private: - char state; -}; - -const compres maybe = (char)2; - -struct CodeTree; - -class SubTree -{ - CodeTree *tree; - bool sign; // Only possible when parent is cAdd or cMul - - inline void flipsign() { sign = !sign; } -public: - SubTree(); - SubTree(double value); - SubTree(const SubTree &b); - SubTree(const CodeTree &b); - - ~SubTree(); - const SubTree &operator= (const SubTree &b); - const SubTree &operator= (const CodeTree &b); - - bool getsign() const { return sign; } - - const CodeTree* operator-> () const { return tree; } - const CodeTree& operator* () const { return *tree; } - struct CodeTree* operator-> () { return tree; } - struct CodeTree& operator* () { return *tree; } - - bool operator< (const SubTree& b) const; - bool operator== (const SubTree& b) const; - void Negate(); // Note: Parent must be cAdd - void Invert(); // Note: Parent must be cMul - - void CheckConstNeg(); - void CheckConstInv(); -}; - -bool IsNegate(const SubTree &p1, const SubTree &p2); -bool IsInverse(const SubTree &p1, const SubTree &p2); - -typedef list paramlist; - -struct CodeTreeData -{ - paramlist args; - -private: - unsigned op; // Operation - double value; // In case of cImmed - unsigned var; // In case of cVar - unsigned funcno; // In case of cFCall, cPCall - -public: - CodeTreeData() : op(cAdd) {} - ~CodeTreeData() {} - - void SetOp(unsigned newop) { op=newop; } - void SetFuncNo(unsigned newno) { funcno=newno; } - unsigned GetFuncNo() const { return funcno; } - - bool IsFunc() const { return op == cFCall || op == cPCall; } - bool IsImmed() const { return op == cImmed; } - bool IsVar() const { return op == cVar; } - inline unsigned GetOp() const { return op; } - inline double GetImmed() const - { - return value; - } - inline unsigned GetVar() const - { - return var; - } - - void AddParam(const SubTree &p) - { - args.push_back(p); - } - void SetVar(unsigned v) - { - args.clear(); - op = cVar; - var = v; - } - void SetImmed(double v) - { - args.clear(); - op = cImmed; - value = orig = v; - inverted = negated = false; - } - void NegateImmed() - { - negated = !negated; - UpdateValue(); - } - void InvertImmed() - { - inverted = !inverted; - UpdateValue(); - } - - bool IsInverted() const { return inverted; } - bool IsNegated() const { return negated; } - bool IsNegatedOriginal() const { return !IsInverted() && IsNegated(); } - -private: - void UpdateValue() - { - value = orig; - if(IsInverted()) { value = 1.0 / value; - // FIXME: potential divide by zero. - } - if(IsNegated()) value = -value; - } - - double orig; - bool inverted; - bool negated; -protected: - // Ensure we don't accidentally copy this - void operator=(const CodeTreeData &b); -}; - - -class CodeTreeDataPtr -{ - typedef pair p_t; - typedef p_t* pp; - mutable pp p; - - void Alloc() const { ++p->second; } - void Dealloc() const { if(!--p->second) delete p; p = 0; } - - void PrepareForWrite() - { - // We're ready if we're the only owner. - if(p->second == 1) return; - - // Then make a clone. - p_t *newtree = new p_t(p->first, 1); - // Forget the old - Dealloc(); - // Keep the new - p = newtree; - } - -public: - CodeTreeDataPtr() : p(new p_t) { p->second = 1; } - CodeTreeDataPtr(const CodeTreeDataPtr &b): p(b.p) { Alloc(); } - ~CodeTreeDataPtr() { Dealloc(); } - const CodeTreeDataPtr &operator= (const CodeTreeDataPtr &b) - { - b.Alloc(); - Dealloc(); - p = b.p; - return *this; - } - const CodeTreeData *operator-> () const { return &p->first; } - CodeTreeData *operator-> () { PrepareForWrite(); return &p->first; } - - void Shock(); -}; - - -#define CHECKCONSTNEG(item, op) \ - ((op)==cMul) \ - ? (item).CheckConstInv() \ - : (item).CheckConstNeg() - -struct CodeTree -{ - CodeTreeDataPtr data; - -private: - typedef paramlist::iterator pit; - typedef paramlist::const_iterator pcit; - - /* - template inline void chk() const - { - } - */ - -public: - const pcit GetBegin() const { return data->args.begin(); } - const pcit GetEnd() const { return data->args.end(); } - const pit GetBegin() { return data->args.begin(); } - const pit GetEnd() { return data->args.end(); } - const SubTree& getp0() const { /*chk<1>();*/pcit tmp=GetBegin(); return *tmp; } - const SubTree& getp1() const { /*chk<2>();*/pcit tmp=GetBegin(); ++tmp; return *tmp; } - const SubTree& getp2() const { /*chk<3>();*/pcit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; } - unsigned GetArgCount() const { return data->args.size(); } - void Erase(const pit p) { data->args.erase(p); } - - SubTree& getp0() { /*chk<1>();*/pit tmp=GetBegin(); return *tmp; } - SubTree& getp1() { /*chk<2>();*/pit tmp=GetBegin(); ++tmp; return *tmp; } - - // set - void SetImmed(double v) { data->SetImmed(v); } - void SetOp(unsigned op) { data->SetOp(op); } - void SetVar(unsigned v) { data->SetVar(v); } - // get - double GetImmed() const { return data->GetImmed(); } - unsigned GetVar() const { return data->GetVar(); } - unsigned GetOp() const { return data->GetOp(); } - // test - bool IsImmed() const { return data->IsImmed(); } - bool IsVar() const { return data->IsVar(); } - // act - void AddParam(const SubTree &p) { data->AddParam(p); } - void NegateImmed() { data->NegateImmed(); } // don't use when op!=cImmed - void InvertImmed() { data->InvertImmed(); } // don't use when op!=cImmed - - compres NonZero() const { if(!IsImmed()) return maybe; - return GetImmed() != 0.0; } - -private: - struct ConstList - { - double voidvalue; - list cp; - double value; - unsigned size() const { return cp.size(); } - }; - struct ConstList BuildConstList(); - void KillConst(const ConstList &cl) - { - for(list::const_iterator i=cl.cp.begin(); i!=cl.cp.end(); ++i) - Erase(*i); - } - void FinishConst(const ConstList &cl) - { - if(cl.value != cl.voidvalue && cl.size() > 1) AddParam(cl.value); - if(cl.value == cl.voidvalue || cl.size() > 1) KillConst(cl); - } - -public: - CodeTree() {} - CodeTree(double v) { SetImmed(v); } - - CodeTree(unsigned op, const SubTree &p) - { - SetOp(op); - AddParam(p); - } - CodeTree(unsigned op, const SubTree &p1, const SubTree &p2) - { - SetOp(op); - AddParam(p1); - AddParam(p2); - } - - bool operator== (const CodeTree& b) const; - bool operator< (const CodeTree& b) const; - -private: - bool IsSortable() const - { - switch(GetOp()) - { - case cAdd: case cMul: - case cEqual: - case cAnd: case cOr: - case cMax: case cMin: - return true; - default: - return false; - } - } - void SortIfPossible() - { - if(IsSortable()) - { - data->args.sort(); - } - } - - void ReplaceWithConst(double value) - { - SetImmed(value); - - /* REMEMBER TO CALL CheckConstInv / CheckConstNeg - * FOR PARENT SubTree, OR MAYHEM HAPPENS - */ - } - - void ReplaceWith(const CodeTree &b) - { - // If b is child of *this, mayhem - // happens. So we first make a clone - // and then proceed with copy. - CodeTreeDataPtr tmp = b.data; - tmp.Shock(); - data = tmp; - } - - void ReplaceWith(unsigned op, const SubTree &p) - { - ReplaceWith(CodeTree(op, p)); - } - - void ReplaceWith(unsigned op, const SubTree &p1, const SubTree &p2) - { - ReplaceWith(CodeTree(op, p1, p2)); - } - - void OptimizeConflict() - { - // This optimization does this: x-x = 0, x/x = 1, a+b-a = b. - - if(GetOp() == cAdd || GetOp() == cMul) - { - Redo: - pit a, b; - for(a=GetBegin(); a!=GetEnd(); ++a) - { - for(b=GetBegin(); ++b != GetEnd(); ) - { - const SubTree &p1 = *a; - const SubTree &p2 = *b; - - if(GetOp() == cMul ? IsInverse(p1,p2) - : IsNegate(p1,p2)) - { - // These parameters complement each others out - Erase(b); - Erase(a); - goto Redo; - } - } - } - } - OptimizeRedundant(); - } - - void OptimizeRedundant() - { - // This optimization does this: min()=0, max()=0, add()=0, mul()=1 - - if(!GetArgCount()) - { - if(GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax) - ReplaceWithConst(0); - else if(GetOp() == cMul) - ReplaceWithConst(1); - return; - } - - // And this: mul(x) = x, min(x) = x, max(x) = x, add(x) = x - - if(GetArgCount() == 1) - { - if(GetOp() == cMul || GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax) - if(!getp0().getsign()) - ReplaceWith(*getp0()); - } - - OptimizeDoubleNegations(); - } - - void OptimizeDoubleNegations() - { - if(GetOp() == cAdd) - { - // Eschew double negations - - // If any of the elements is cMul - // and has a numeric constant, negate - // the constant and negate sign. - - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - SubTree &pa = *a; - if(pa.getsign() - && pa->GetOp() == cMul) - { - CodeTree &p = *pa; - for(pit b=p.GetBegin(); - b!=p.GetEnd(); ++b) - { - SubTree &pb = *b; - if(pb->IsImmed()) - { - pb.Negate(); - pa.Negate(); - break; - } - } - } - } - } - - if(GetOp() == cMul) - { - // If any of the elements is cPow - // and has a numeric exponent, negate - // the exponent and negate sign. - - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - SubTree &pa = *a; - if(pa.getsign() && pa->GetOp() == cPow) - { - CodeTree &p = *pa; - if(p.getp1()->IsImmed()) - { - // negate ok for pow when op=cImmed - p.getp1().Negate(); - pa.Negate(); - } - } - } - } - } - - void OptimizeConstantMath1() - { - // This optimization does three things: - // - For adding groups: - // Constants are added together. - // - For multiplying groups: - // Constants are multiplied together. - // - For function calls: - // If all parameters are constants, - // the call is replaced with constant value. - - // First, do this: - OptimizeAddMulFlat(); - - switch(GetOp()) - { - case cAdd: - { - ConstList cl = BuildConstList(); - FinishConst(cl); - break; - } - case cMul: - { - ConstList cl = BuildConstList(); - - if(cl.value == 0.0) ReplaceWithConst(0.0); - else FinishConst(cl); - - break; - } - #define ConstantUnaryFun(token, fun) \ - case token: { const SubTree &p0 = getp0(); \ - if(p0->IsImmed()) ReplaceWithConst(fun(p0->GetImmed())); \ - break; } - #define ConstantBinaryFun(token, fun) \ - case token: { const SubTree &p0 = getp0(); \ - const SubTree &p1 = getp1(); \ - if(p0->IsImmed() && \ - p1->IsImmed()) ReplaceWithConst(fun(p0->GetImmed(), p1->GetImmed())); \ - break; } - - // FIXME: potential invalid parameters for functions - // can cause exceptions here - - ConstantUnaryFun(cAbs, fabs); - ConstantUnaryFun(cAcos, acos); - ConstantUnaryFun(cAsin, asin); - ConstantUnaryFun(cAtan, atan); - ConstantUnaryFun(cCeil, ceil); - ConstantUnaryFun(cCos, cos); - ConstantUnaryFun(cCosh, cosh); - ConstantUnaryFun(cFloor, floor); - ConstantUnaryFun(cLog, log); - ConstantUnaryFun(cSin, sin); - ConstantUnaryFun(cSinh, sinh); - ConstantUnaryFun(cTan, tan); - ConstantUnaryFun(cTanh, tanh); - ConstantBinaryFun(cAtan2, atan2); - ConstantBinaryFun(cMax, Max); - ConstantBinaryFun(cMin, Min); - ConstantBinaryFun(cMod, fmod); // not a func, but belongs here too - ConstantBinaryFun(cPow, pow); - - case cNeg: - case cSub: - case cDiv: - /* Unreached (nonexistent operator) - * TODO: internal error here? - */ - break; - - case cCot: - case cCsc: - case cSec: - case cDeg: - case cRad: - case cLog10: - case cSqrt: - case cExp: - /* Unreached (nonexistent function) - * TODO: internal error here? - */ - break; - } - - OptimizeConflict(); - } - - void OptimizeAddMulFlat() - { - // This optimization flattens the topography of the tree. - // Examples: - // x + (y+z) = x+y+z - // x * (y/z) = x*y/z - // x / (y/z) = x/y*z - - if(GetOp() == cAdd || GetOp() == cMul) - { - // If children are same type as parent add them here - for(pit b, a=GetBegin(); a!=GetEnd(); a=b) - { - const SubTree &pa = *a; b=a; ++b; - if(pa->GetOp() != GetOp()) continue; - - // Child is same type - for(pcit c=pa->GetBegin(); - c!=pa->GetEnd(); - ++c) - { - const SubTree &pb = *c; - if(pa.getsign()) - { - // +a -(+b +c) - // means b and c will be negated - - SubTree tmp = pb; - if(GetOp() == cMul) - tmp.Invert(); - else - tmp.Negate(); - AddParam(tmp); - } - else - AddParam(pb); - } - Erase(a); - - // Note: OptimizeConstantMath1() would be a good thing to call next. - } - } - } - - void OptimizeLinearCombine() - { - // This optimization does the following: - // - // x*x*x*x -> x^4 - // x+x+x+x -> x*4 - // x*x -> x^2 - // x/z/z -> - // - - // Remove conflicts first, so we don't have to worry about signs. - OptimizeConflict(); - - bool didchanges = false; - if(GetOp() == cAdd || GetOp() == cMul) - { - Redo: - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - - list poslist; - - for(pit b=a; ++b!=GetEnd(); ) - { - const SubTree &pb = *b; - if(*pa == *pb) - poslist.push_back(b); - } - - unsigned min = 2; - if(poslist.size() >= min) - { - SubTree arvo = pa; - bool negate = arvo.getsign(); - - double factor = poslist.size() + 1; - - if(negate) - { - arvo.Negate(); - factor = -factor; - } - - CodeTree tmp(GetOp()==cAdd ? cMul : cPow, - arvo, - factor); - - list::const_iterator j; - for(j=poslist.begin(); j!=poslist.end(); ++j) - Erase(*j); - poslist.clear(); - - *a = tmp; - didchanges = true; - goto Redo; - } - } - } - if(didchanges) - { - // As a result, there might be need for this: - OptimizeAddMulFlat(); - // And this: - OptimizeRedundant(); - } - } - - void OptimizeLogarithm() - { - /* - This is basic logarithm math: - pow(X,Y)/log(Y) = X - log(X)/log(Y) = logY(X) - log(X^Y) = log(X)*Y - log(X*Y) = log(X)+log(Y) - exp(log(X)*Y) = X^Y - - This function does these optimizations: - pow(const_E, log(x)) = x - pow(const_E, log(x)*y) = x^y - pow(10, log(x)*const_L10I*y) = x^y - pow(z, log(x)/log(z)*y) = x^y - - And this: - log(x^z) = z * log(x) - Which automatically causes these too: - log(pow(const_E, x)) = x - log(pow(y, x)) = x * log(y) - log(pow(pow(const_E, y), x)) = x*y - - And it does this too: - log(x) + log(y) + log(z) = log(x * y * z) - log(x * exp(y)) = log(x) + y - - */ - - // Must be already in exponential form. - - // Optimize exponents before doing something. - OptimizeExponents(); - - if(GetOp() == cLog) - { - // We should have one parameter for log() function. - // If we don't, we're screwed. - - const SubTree &p = getp0(); - - if(p->GetOp() == cPow) - { - // Found log(x^y) - SubTree p0 = p->getp0(); // x - SubTree p1 = p->getp1(); // y - - // Build the new logarithm. - CodeTree tmp(GetOp(), p0); // log(x) - - // Become log(x) * y - ReplaceWith(cMul, tmp, p1); - } - else if(p->GetOp() == cMul) - { - // Redefine &p nonconst - SubTree &p = getp0(); - - p->OptimizeAddMulFlat(); - p->OptimizeExponents(); - CHECKCONSTNEG(p, p->GetOp()); - - list adds; - - for(pit b, a = p->GetBegin(); - a != p->GetEnd(); a=b) - { - SubTree &pa = *a; b=a; ++b; - if(pa->GetOp() == cPow - && pa->getp0()->IsImmed() - && pa->getp0()->GetImmed() == CONSTANT_E) - { - adds.push_back(pa->getp1()); - p->Erase(a); - continue; - } - } - if(adds.size()) - { - CodeTree tmp(cAdd, *this); - - list::const_iterator i; - for(i=adds.begin(); i!=adds.end(); ++i) - tmp.AddParam(*i); - - ReplaceWith(tmp); - } - } - } - if(GetOp() == cAdd) - { - // Check which ones are logs. - list poslist; - - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - if(pa->GetOp() == cLog) - poslist.push_back(a); - } - - if(poslist.size() >= 2) - { - CodeTree tmp(cMul, 1.0); // eek - - list::const_iterator j; - for(j=poslist.begin(); j!=poslist.end(); ++j) - { - const SubTree &pb = **j; - // Take all of its children - for(pcit b=pb->GetBegin(); - b!=pb->GetEnd(); - ++b) - { - SubTree tmp2 = *b; - if(pb.getsign()) tmp2.Negate(); - tmp.AddParam(tmp2); - } - Erase(*j); - } - poslist.clear(); - - AddParam(CodeTree(cLog, tmp)); - } - // Done, hopefully - } - if(GetOp() == cPow) - { - const SubTree &p0 = getp0(); - SubTree &p1 = getp1(); - - if(p0->IsImmed() && p0->GetImmed() == CONSTANT_E - && p1->GetOp() == cLog) - { - // pow(const_E, log(x)) = x - ReplaceWith(*(p1->getp0())); - } - else if(p1->GetOp() == cMul) - { - //bool didsomething = true; - - pit poslogpos; bool foundposlog = false; - pit neglogpos; bool foundneglog = false; - - ConstList cl = p1->BuildConstList(); - - for(pit a=p1->GetBegin(); a!=p1->GetEnd(); ++a) - { - const SubTree &pa = *a; - if(pa->GetOp() == cLog) - { - if(!pa.getsign()) - { - foundposlog = true; - poslogpos = a; - } - else if(*p0 == *(pa->getp0())) - { - foundneglog = true; - neglogpos = a; - } - } - } - - if(p0->IsImmed() - && p0->GetImmed() == 10.0 - && cl.value == CONSTANT_L10I - && foundposlog) - { - SubTree base = (*poslogpos)->getp0(); - p1->KillConst(cl); - p1->Erase(poslogpos); - p1->OptimizeRedundant(); - SubTree mul = p1; - - ReplaceWith(cPow, base, mul); - - // FIXME: what optimizations should be done now? - return; - } - - // Put back the constant - FinishConst(cl); - - if(p0->IsImmed() - && p0->GetImmed() == CONSTANT_E - && foundposlog) - { - SubTree base = (*poslogpos)->getp0(); - p1->Erase(poslogpos); - - p1->OptimizeRedundant(); - SubTree mul = p1; - - ReplaceWith(cPow, base, mul); - - // FIXME: what optimizations should be done now? - return; - } - - if(foundposlog - && foundneglog - && *((*neglogpos)->getp0()) == *p0) - { - SubTree base = (*poslogpos)->getp0(); - p1->Erase(poslogpos); - p1->Erase(neglogpos); - - p1->OptimizeRedundant(); - SubTree mul = p1; - - ReplaceWith(cPow, base, mul); - - // FIXME: what optimizations should be done now? - return; - } - } - } - } - - void OptimizeFunctionCalls() - { - /* Goals: sin(asin(x)) = x - * cos(acos(x)) = x - * tan(atan(x)) = x - * NOTE: - * Do NOT do these: - * asin(sin(x)) - * acos(cos(x)) - * atan(tan(x)) - * Because someone might want to wrap the angle. - */ - // FIXME: TODO - } - - void OptimizePowMulAdd() - { - // x^3 * x -> x^4 - // x*3 + x -> x*4 - // FIXME: Do those - - // x^1 -> x - if(GetOp() == cPow) - { - const SubTree &base = getp0(); - const SubTree &exponent = getp1(); - - if(exponent->IsImmed()) - { - if(exponent->GetImmed() == 1.0) - ReplaceWith(*base); - else if(exponent->GetImmed() == 0.0 - && base->NonZero()) - ReplaceWithConst(1.0); - } - } - } - - void OptimizeExponents() - { - /* Goals: - * (x^y)^z -> x^(y*z) - * x^y * x^z -> x^(y+z) - */ - // First move to exponential form. - OptimizeLinearCombine(); - - bool didchanges = false; - - Redo: - if(GetOp() == cPow) - { - // (x^y)^z -> x^(y*z) - - const SubTree &p0 = getp0(); - const SubTree &p1 = getp1(); - if(p0->GetOp() == cPow) - { - CodeTree tmp(cMul, p0->getp1(), p1); - tmp.Optimize(); - - ReplaceWith(cPow, p0->getp0(), tmp); - - didchanges = true; - goto Redo; - } - } - if(GetOp() == cMul) - { - // x^y * x^z -> x^(y+z) - - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - - if(pa->GetOp() != cPow) continue; - - list poslist; - - for(pit b=a; ++b != GetEnd(); ) - { - const SubTree &pb = *b; - if(pb->GetOp() == cPow - && *(pa->getp0()) - == *(pb->getp0())) - { - poslist.push_back(b); - } - } - - if(poslist.size() >= 1) - { - poslist.push_back(a); - - CodeTree base = *(pa->getp0()); - - CodeTree exponent(cAdd, 0.0); //eek - - // Collect all exponents to cAdd - list::const_iterator i; - for(i=poslist.begin(); i!=poslist.end(); ++i) - { - const SubTree &pb = **i; - - SubTree tmp2 = pb->getp1(); - if(pb.getsign()) tmp2.Invert(); - - exponent.AddParam(tmp2); - } - - exponent.Optimize(); - - CodeTree result(cPow, base, exponent); - - for(i=poslist.begin(); i!=poslist.end(); ++i) - Erase(*i); - poslist.clear(); - - AddParam(result); // We're cMul, remember - - didchanges = true; - goto Redo; - } - } - } - - OptimizePowMulAdd(); - - if(didchanges) - { - // As a result, there might be need for this: - OptimizeConflict(); - } - } - - void OptimizeLinearExplode() - { - // x^2 -> x*x - // But only if x is just a simple thing - - // Won't work on anything else. - if(GetOp() != cPow) return; - - // TODO TODO TODO - } - - void OptimizePascal() - { -#if 0 // Too big, too specific, etc - - // Won't work on anything else. - if(GetOp() != cAdd) return; - - // Must be done after OptimizeLinearCombine(); - - // Don't need pascal triangle - // Coefficient for x^a * y^b * z^c = 3! / (a! * b! * c!) - - // We are greedy and want other than just binomials - // FIXME - - // note: partial ones are also nice - // x*x + x*y + y*y - // = (x+y)^2 - x*y - // - // x x * x y * + y y * + - // -> x y + dup * x y * - -#endif - } - -public: - - void Optimize(); - - void Assemble(vector &byteCode, - vector &immed) const; - - void FinalOptimize() - { - // First optimize each parameter. - for(pit a=GetBegin(); a!=GetEnd(); ++a) - (*a)->FinalOptimize(); - - /* These things are to be done: - * - * x * CONSTANT_DR -> cDeg(x) - * x * CONSTANT_RD -> cRad(x) - * pow(x, 0.5) -> sqrt(x) - * log(x) * CONSTANT_L10I -> log10(x) - * pow(CONSTANT_E, x) -> exp(x) - * inv(sin(x)) -> csc(x) - * inv(cos(x)) -> sec(x) - * inv(tan(x)) -> cot(x) - */ - - - if(GetOp() == cPow) - { - const SubTree &p0 = getp0(); - const SubTree &p1 = getp1(); - if(p0->GetOp() == cImmed - && p0->GetImmed() == CONSTANT_E) - { - ReplaceWith(cExp, p1); - } - else if(p1->GetOp() == cImmed - && p1->GetImmed() == 0.5) - { - ReplaceWith(cSqrt, p0); - } - } - if(GetOp() == cMul) - { - if(GetArgCount() == 1 && getp0().getsign()) - { - /***/if(getp0()->GetOp() == cSin)ReplaceWith(cCsc, getp0()->getp0()); - else if(getp0()->GetOp() == cCos)ReplaceWith(cSec, getp0()->getp0()); - else if(getp0()->GetOp() == cTan)ReplaceWith(cCot, getp0()->getp0()); - } - } - // Separate "if", because op may have just changed - if(GetOp() == cMul) - { - CodeTree *found_log = 0; - - ConstList cl = BuildConstList(); - - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - SubTree &pa = *a; - if(pa->GetOp() == cLog && !pa.getsign()) - found_log = &*pa; - } - if(cl.value == CONSTANT_L10I && found_log) - { - // Change the log() to log10() - found_log->SetOp(cLog10); - // And forget the constant - KillConst(cl); - } - else if(cl.value == CONSTANT_DR) - { - OptimizeRedundant(); - ReplaceWith(cDeg, *this); - } - else if(cl.value == CONSTANT_RD) - { - OptimizeRedundant(); - ReplaceWith(cRad, *this); - } - else FinishConst(cl); - } - - SortIfPossible(); - } -}; - -void CodeTreeDataPtr::Shock() -{ - /* - PrepareForWrite(); - paramlist &p2 = (*this)->args; - for(paramlist::iterator i=p2.begin(); i!=p2.end(); ++i) - { - (*i)->data.Shock(); - } - */ -} - -CodeTree::ConstList CodeTree::BuildConstList() -{ - ConstList result; - result.value = - result.voidvalue = GetOp()==cMul ? 1.0 : 0.0; - - list &cp = result.cp; - for(pit b, a=GetBegin(); a!=GetEnd(); a=b) - { - SubTree &pa = *a; b=a; ++b; - if(!pa->IsImmed()) continue; - - double thisvalue = pa->GetImmed(); - if(thisvalue == result.voidvalue) - { - // This value is no good, forget it - Erase(a); - continue; - } - if(GetOp() == cMul) - result.value *= thisvalue; - else - result.value += thisvalue; - cp.push_back(a); - } - if(GetOp() == cMul) - { - /* - Jos joku niistä arvoista on -1 eikä se ole ainoa arvo, - niin joku muu niistä arvoista negatoidaan. - */ - for(bool done=false; cp.size() > 1 && !done; ) - { - done = true; - for(list::iterator b,a=cp.begin(); a!=cp.end(); a=b) - { - b=a; ++b; - if((**a)->GetImmed() == -1.0) - { - Erase(*a); - cp.erase(a); - - // take randomly something - (**cp.begin())->data->NegateImmed(); - if(cp.size() < 2)break; - done = false; - } - } - } - } - return result; -} - -void CodeTree::Assemble - (vector &byteCode, - vector &immed) const -{ - #define AddCmd(op) byteCode.push_back((op)) - #define AddConst(v) do { \ - byteCode.push_back(cImmed); \ - immed.push_back((v)); \ - } while(0) - - if(IsVar()) - { - AddCmd(GetVar()); - return; - } - if(IsImmed()) - { - AddConst(GetImmed()); - return; - } - - switch(GetOp()) - { - case cAdd: - case cMul: - { - unsigned opcount = 0; - for(pcit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - - if(opcount < 2) ++opcount; - - bool pnega = pa.getsign(); - - bool done = false; - if(pa->IsImmed()) - { - if(GetOp() == cMul - && pa->data->IsInverted() - && (pnega || opcount==2) - ) - { - CodeTree tmp = *pa; - tmp.data->InvertImmed(); - tmp.Assemble(byteCode, immed); - pnega = !pnega; - done = true; - } - else if(GetOp() == cAdd - && (pa->data->IsNegatedOriginal() - // || pa->GetImmed() < 0 - ) - && (pnega || opcount==2) - ) - { - CodeTree tmp = *pa; - tmp.data->NegateImmed(); - tmp.Assemble(byteCode, immed); - pnega = !pnega; - done = true; - } - } - if(!done) - pa->Assemble(byteCode, immed); - - if(opcount == 2) - { - unsigned tmpop = GetOp(); - if(pnega) // negate - { - tmpop = (tmpop == cMul) ? cDiv : cSub; - } - AddCmd(tmpop); - } - else if(pnega) - { - if(GetOp() == cMul) AddCmd(cInv); - else AddCmd(cNeg); - } - } - break; - } - case cIf: - { - // If the parameter amount is != 3, we're screwed. - getp0()->Assemble(byteCode, immed); - - unsigned ofs = byteCode.size(); - AddCmd(cIf); - AddCmd(0); // code index - AddCmd(0); // immed index - - getp1()->Assemble(byteCode, immed); - - byteCode[ofs+1] = byteCode.size()+2; - byteCode[ofs+2] = immed.size(); - - ofs = byteCode.size(); - AddCmd(cJump); - AddCmd(0); // code index - AddCmd(0); // immed index - - getp2()->Assemble(byteCode, immed); - - byteCode[ofs+1] = byteCode.size()-1; - byteCode[ofs+2] = immed.size(); - - break; - } - case cFCall: - { - // If the parameter count is invalid, we're screwed. - for(pcit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - pa->Assemble(byteCode, immed); - } - AddCmd(GetOp()); - AddCmd(data->GetFuncNo()); - break; - } - case cPCall: - { - // If the parameter count is invalid, we're screwed. - for(pcit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - pa->Assemble(byteCode, immed); - } - AddCmd(GetOp()); - AddCmd(data->GetFuncNo()); - break; - } - default: - { - // If the parameter count is invalid, we're screwed. - for(pcit a=GetBegin(); a!=GetEnd(); ++a) - { - const SubTree &pa = *a; - pa->Assemble(byteCode, immed); - } - AddCmd(GetOp()); - break; - } - } -} - -void CodeTree::Optimize() -{ - // Phase: - // Phase 0: Do local optimizations. - // Phase 1: Optimize each. - // Phase 2: Do local optimizations again. - - for(unsigned phase=0; phase<=2; ++phase) - { - if(phase == 1) - { - // Optimize each parameter. - for(pit a=GetBegin(); a!=GetEnd(); ++a) - { - (*a)->Optimize(); - CHECKCONSTNEG(*a, GetOp()); - } - continue; - } - if(phase == 0 || phase == 2) - { - // Do local optimizations. - - OptimizeConstantMath1(); - OptimizeLogarithm(); - OptimizeFunctionCalls(); - OptimizeExponents(); - OptimizeLinearExplode(); - OptimizePascal(); - - /* Optimization paths: - - doublenegations= - redundant= * doublenegations - conflict= * redundant - addmulflat= - constantmath1= addmulflat * conflict - linearcombine= conflict * addmulflat¹ redundant¹ - powmuladd= - exponents= linearcombine * powmuladd conflict¹ - logarithm= exponents * - functioncalls= IDLE - linearexplode= IDLE - pascal= IDLE - - * = actions here - ¹ = only if made changes - */ - } - } -} - - -bool CodeTree::operator== (const CodeTree& b) const -{ - if(GetOp() != b.GetOp()) return false; - if(IsImmed()) if(GetImmed() != b.GetImmed()) return false; - if(IsVar()) if(GetVar() != b.GetVar()) return false; - if(data->IsFunc()) - if(data->GetFuncNo() != b.data->GetFuncNo()) return false; - return data->args == b.data->args; -} - -bool CodeTree::operator< (const CodeTree& b) const -{ - if(GetArgCount() != b.GetArgCount()) - return GetArgCount() > b.GetArgCount(); - - if(GetOp() != b.GetOp()) - { - // sort immeds last - if(IsImmed() != b.IsImmed()) return IsImmed() < b.IsImmed(); - - return GetOp() < b.GetOp(); - } - - if(IsImmed()) - { - if(GetImmed() != b.GetImmed()) return GetImmed() < b.GetImmed(); - } - if(IsVar() && GetVar() != b.GetVar()) - { - return GetVar() < b.GetVar(); - } - if(data->IsFunc() && data->GetFuncNo() != b.data->GetFuncNo()) - { - return data->GetFuncNo() < b.data->GetFuncNo(); - } - - pcit i = GetBegin(), j = b.GetBegin(); - for(; i != GetEnd(); ++i, ++j) - { - const SubTree &pa = *i, &pb = *j; - - if(!(pa == pb)) - return pa < pb; - } - return false; -} - - -bool IsNegate(const SubTree &p1, const SubTree &p2) /*const */ -{ - if(p1->IsImmed() && p2->IsImmed()) - { - return p1->GetImmed() == -p2->GetImmed(); - } - if(p1.getsign() == p2.getsign()) return false; - return *p1 == *p2; -} -bool IsInverse(const SubTree &p1, const SubTree &p2) /*const*/ -{ - if(p1->IsImmed() && p2->IsImmed()) - { - // FIXME: potential divide by zero. - return p1->GetImmed() == 1.0 / p2->GetImmed(); - } - if(p1.getsign() == p2.getsign()) return false; - return *p1 == *p2; -} - -SubTree::SubTree() : tree(new CodeTree), sign(false) -{ -} - -SubTree::SubTree(const SubTree &b) : tree(new CodeTree(*b.tree)), sign(b.sign) -{ -} - -#define SubTreeDecl(p1, p2) \ - SubTree::SubTree p1 : tree(new CodeTree p2), sign(false) { } - -SubTreeDecl( (const CodeTree &b), (b) ) -SubTreeDecl( (double value), (value) ) - -#undef SubTreeDecl - -SubTree::~SubTree() -{ - delete tree; tree=0; -} - -const SubTree &SubTree::operator= (const SubTree &b) -{ - sign = b.sign; - CodeTree *oldtree = tree; - tree = new CodeTree(*b.tree); - delete oldtree; - return *this; -} -const SubTree &SubTree::operator= (const CodeTree &b) -{ - sign = false; - CodeTree *oldtree = tree; - tree = new CodeTree(b); - delete oldtree; - return *this; -} - -bool SubTree::operator< (const SubTree& b) const -{ - if(getsign() != b.getsign()) return getsign() < b.getsign(); - return *tree < *b.tree; -} -bool SubTree::operator== (const SubTree& b) const -{ - return sign == b.sign && *tree == *b.tree; -} -void SubTree::Negate() // Note: Parent must be cAdd -{ - flipsign(); - CheckConstNeg(); -} -void SubTree::CheckConstNeg() -{ - if(tree->IsImmed() && getsign()) - { - tree->NegateImmed(); - sign = false; - } -} -void SubTree::Invert() // Note: Parent must be cMul -{ - flipsign(); - CheckConstInv(); -} -void SubTree::CheckConstInv() -{ - if(tree->IsImmed() && getsign()) - { - tree->InvertImmed(); - sign = false; - } -} - -}//namespace - -void FunctionParser::MakeTree(void *r) const -{ - // Dirty hack. Should be fixed. - CodeTree* result = static_cast(r); - - vector stack(1); - - #define GROW(n) do { \ - stacktop += n; \ - if(stack.size() <= stacktop) stack.resize(stacktop+1); \ - } while(0) - - #define EAT(n, opcode) do { \ - unsigned newstacktop = stacktop-n; \ - stack[stacktop].SetOp((opcode)); \ - for(unsigned a=0, b=(n); a labels; - - const unsigned* const ByteCode = data->ByteCode; - const unsigned ByteCodeSize = data->ByteCodeSize; - const double* const Immed = data->Immed; - - for(unsigned IP=0, DP=0; ; ++IP) - { - while(labels.size() > 0 - && *labels.begin() == IP) - { - // The "else" of an "if" ends here - EAT(3, cIf); - labels.erase(labels.begin()); - } - - if(IP >= ByteCodeSize) - { - break; - } - - unsigned opcode = ByteCode[IP]; - - if(opcode == cIf) - { - IP += 2; - } - else if(opcode == cJump) - { - labels.push_front(ByteCode[IP+1]+1); - IP += 2; - } - else if(opcode == cImmed) - { - ADDCONST(Immed[DP++]); - } - else if(opcode < VarBegin) - { - switch(opcode) - { - // Unary operators - case cNeg: - { - EAT(1, cAdd); // Unary minus is negative adding. - stack[stacktop-1].getp0().Negate(); - break; - } - // Binary operators - case cSub: - { - EAT(2, cAdd); // Minus is negative adding - stack[stacktop-1].getp1().Negate(); - break; - } - case cDiv: - { - EAT(2, cMul); // Divide is inverse multiply - stack[stacktop-1].getp1().Invert(); - break; - } - - // ADD ALL TWO PARAMETER NON-FUNCTIONS HERE - case cAdd: case cMul: - case cMod: case cPow: - case cEqual: case cLess: case cGreater: - case cAnd: case cOr: - EAT(2, opcode); - break; - - case cFCall: - { - unsigned index = ByteCode[++IP]; - unsigned params = data->FuncPtrs[index].params; - EAT(params, opcode); - stack[stacktop-1].data->SetFuncNo(index); - break; - } - case cPCall: - { - unsigned index = ByteCode[++IP]; - unsigned params = - data->FuncParsers[index]->data->varAmount; - EAT(params, opcode); - stack[stacktop-1].data->SetFuncNo(index); - break; - } - - // Converted to cMul on fly - case cDeg: - ADDCONST(CONSTANT_DR); - EAT(2, cMul); - break; - - // Converted to cMul on fly - case cRad: - ADDCONST(CONSTANT_RD); - EAT(2, cMul); - break; - - // Functions - default: - { - const FuncDefinition& func = Functions[opcode-cAbs]; - - unsigned paramcount = func.params; -#ifndef DISABLE_EVAL - if(opcode == cEval) paramcount = data->varAmount; -#endif - if(opcode == cSqrt) - { - // Converted on fly: sqrt(x) = x^0.5 - opcode = cPow; - paramcount = 2; - ADDCONST(0.5); - } - if(opcode == cExp) - { - // Converted on fly: exp(x) = CONSTANT_E^x - - opcode = cPow; - paramcount = 2; - // reverse the parameters... kludgey - stack[stacktop] = stack[stacktop-1]; - stack[stacktop-1].SetImmed(CONSTANT_E); - GROW(1); - } - bool do_inv = false; - if(opcode == cCot) { do_inv = true; opcode = cTan; } - if(opcode == cCsc) { do_inv = true; opcode = cSin; } - if(opcode == cSec) { do_inv = true; opcode = cCos; } - - bool do_log10 = false; - if(opcode == cLog10) - { - // Converted on fly: log10(x) = log(x) * CONSTANT_L10I - opcode = cLog; - do_log10 = true; - } - EAT(paramcount, opcode); - if(do_log10) - { - ADDCONST(CONSTANT_L10I); - EAT(2, cMul); - } - if(do_inv) - { - // Unary cMul, inverted. No need for "1.0" - EAT(1, cMul); - stack[stacktop-1].getp0().Invert(); - } - break; - } - } - } - else - { - stack[stacktop].SetVar(opcode); - GROW(1); - } - } - - if(!stacktop) - { - // ERROR: Stack does not have any values! - return; - } - - --stacktop; // Ignore the last element, it is always nop (cAdd). - - if(stacktop > 0) - { - // ERROR: Stack has too many values! - return; - } - - // Okay, the tree is now stack[0] - *result = stack[0]; -} - -void FunctionParser::Optimize() -{ - copyOnWrite(); - - CodeTree tree; - MakeTree(&tree); - - // Do all sorts of optimizations - tree.Optimize(); - // Last changes before assembly - tree.FinalOptimize(); - - // Now rebuild from the tree. - - vector byteCode; - vector immed; - -#if 0 - byteCode.resize(Comp.ByteCodeSize); - for(unsigned a=0; aByteCode; data->ByteCode = 0; - if((data->ByteCodeSize = byteCode.size()) > 0) - { - data->ByteCode = new unsigned[data->ByteCodeSize]; - for(unsigned a=0; aByteCode[a] = byteCode[a]; - } - - delete[] data->Immed; data->Immed = 0; - if((data->ImmedSize = immed.size()) > 0) - { - data->Immed = new double[data->ImmedSize]; - for(unsigned a=0; aImmed[a] = immed[a]; - } -} - -#else /* !SUPPORT_OPTIMIZER */ -/* keep the linker happy */ +#ifndef SUPPORT_OPTIMIZER void FunctionParser::MakeTree(void *) const {} void FunctionParser::Optimize() { diff --git a/deal.II/contrib/functionparser/fparser.h b/deal.II/contrib/functionparser/fparser.h index c28cf5defa..9ef0d3dc85 100644 --- a/deal.II/contrib/functionparser/fparser.h +++ b/deal.II/contrib/functionparser/fparser.h @@ -1,7 +1,8 @@ /***************************************************************************\ -|* Function parser v2.71 by Warp *| +|* Function parser v2.83 by Warp *| |* ----------------------------- *| |* Parses and evaluates the given function with the given variable values. *| +|* See fparser.txt for details. *| |* *| \***************************************************************************/ @@ -15,10 +16,8 @@ #ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT #include #endif +namespace fparser { -namespace fparser -{ - class FunctionParser { public: @@ -40,6 +39,7 @@ public: inline int EvalError() const { return evalErrorType; } bool AddConstant(const std::string& name, double value); + bool AddUnit(const std::string& name, double value); typedef double (*FunctionPtr)(const double*); @@ -59,18 +59,23 @@ public: FunctionParser& operator=(const FunctionParser&); + void ForceDeepCopy(); + + #ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT // For debugging purposes only: void PrintByteCode(std::ostream& dest) const; #endif - - + // Added by Luca Heltai. For consistency checking. unsigned int NVars() { return data->varAmount; } + + + //======================================================================== private: //======================================================================== @@ -85,13 +90,14 @@ private: unsigned referenceCounter; int varAmount; - bool useDegreeConversion; + bool useDegreeConversion, isOptimized; typedef std::map VarMap_t; VarMap_t Variables; typedef std::map ConstMap_t; ConstMap_t Constants; + ConstMap_t Units; VarMap_t FuncPtrNames; struct FuncPtrData @@ -119,6 +125,7 @@ private: }; Data* data; + unsigned evalRecursionLevel; // Temp data needed in Compile(): unsigned StackPtr; @@ -128,21 +135,21 @@ private: // Private methods: // --------------- - inline void copyOnWrite(); - + void CopyOnWrite(); - bool checkRecursiveLinking(const FunctionParser*) const; - bool isValidName(const std::string&) const; + bool CheckRecursiveLinking(const FunctionParser*) const; - inline - Data::VarMap_t::const_iterator FindVariable(const char*, - const Data::VarMap_t&) const; - - inline - Data::ConstMap_t::const_iterator FindConstant(const char*) const; + bool ParseVars(const std::string& Vars, + std::map& dest); + int VarNameType(const std::string&) const; + Data::VarMap_t::const_iterator + FindVariable(const char*, const Data::VarMap_t&) const; + Data::ConstMap_t::const_iterator + FindConstant(const char*, const Data::ConstMap_t&) const; + int CheckForUnit(const char*, int) const; int CheckSyntax(const char*); - bool Compile(const char*); + int Compile(const char*); bool IsVariable(int); void AddCompiledByte(unsigned); void AddImmediate(double); @@ -151,6 +158,7 @@ private: int CompileIf(const char*, int); int CompileFunctionParams(const char*, int, unsigned); int CompileElement(const char*, int); + int CompilePossibleUnit(const char*, int); int CompilePow(const char*, int); int CompileUnaryMinus(const char*, int); int CompileMult(const char*, int); @@ -163,8 +171,6 @@ private: void MakeTree(void*) const; }; - -} - +} #endif diff --git a/deal.II/contrib/functionparser/fparser.txt b/deal.II/contrib/functionparser/fparser.txt index d54e5ec260..27803faa99 100644 --- a/deal.II/contrib/functionparser/fparser.txt +++ b/deal.II/contrib/functionparser/fparser.txt @@ -1,4 +1,4 @@ - Function parser for C++ v2.71 by Warp. + Function parser for C++ v2.83 by Warp. ====================================== Optimization code contributed by Bisqwit (http://iki.fi/bisqwit/) @@ -7,31 +7,34 @@ The usage license of this library is located at the end of this text file. + What's new in v2.83 + ------------------- + - Fixed bug in the optimizer which caused functions like "x-y*1" to crash. + + What's new in v2.82 + ------------------- + - Fixed bug in the optimizer which caused a crash with some functions. - What's new in v2.71 + What's new in v2.81 ------------------- - - Fixed a couple of syntax errors which slipped through in the v2.7 - distribution version: - * non-standard semicolons at the end of namespace blocks - * a compilation error which happens if SUPPORT_OPTIMIZER is disabled - - Included a simple example program. - - What's new in v2.7 - ------------------ - - Changed precedence rules for unary minus and the power operator (^) to - make it closer in functionality to the power "operator" in mathematics - (ie. superscript): - * Consecutive power operators at the same precedence level are - evaluated from right to left. That is, for example "2^3^4" is - now evaluated as if it had been written as "2^(3^4)" (and not - as "(2^3)^4" like in previous versions). - * The unary minus in the base of the power has now a lower precedence - than the power operator. That is, "-2^3" will be evaluated as if - written as "-(2^3)", ie. the result is "-8" (and not "8" like in - previous versions). The unary minus in the exponent is still - evaluated first because of the right-left precedence change above - (that is, "-2^-3" is evaluated as "-(2^(-3))"). - - Fixed a bug in the copy-on-write engine. + - Added support for units (see help on the AddUnit() method). + - Added support for optionally making the Eval() function thread-safe + with a precompiler constant. + + + Table of contents: + ----------------- + - Preface + - Usage + * Configuring the compilation + * Copying and assignment + * Short descriptions of FunctionParser methods + * Long descriptions of FunctionParser methods + - About thread safety + - The function string + - Contacting the author + - The algorithm used in the library + - Usage license @@ -63,31 +66,44 @@ compiler. - Usage ============================================================================= - To use the FunctionParser class, you have to include "fparser.hh". When -compiling, you have to compile fparser.cc and link it to the main program. -You can also make a library from the fparser.cc (see the help on your -compiler to see how this is done). + To use the FunctionParser class, you have to include "fparser.hh" in +your source code files which use the FunctionParser class. + + When compiling, you have to compile fparser.cc and fpoptimizer.cc and +link them to the main program. In some developement environments it's +enough to add those two files to your current project. + + If you are not going to use the optimizer (ie. you have commented out +SUPPORT_OPTIMIZER in fpconfig.hh), you can leave the latter file out. - * Conditional compiling: - --------------------- + * Configuring the compilation: + --------------------------- - There is a set of precompiler options at the beginning of fparser.cc - which can be used for setting certain features on or off. These lines - can be commented or uncommented depending on the desired behaviour: + There is a set of precompiler options in the fpconfig.hh file + which can be used for setting certain features on or off: NO_ASINH : (Default on) By default the library does not support the asinh(), acosh() and atanh() functions because they are not part of the ISO C++ standard. If your compiler supports them and you want the - parser to support them as well, comment this line. + parser to support them as well, comment out this line. DISABLE_EVAL : (Default off) - The eval() function can be dangerous because it can cause an - infinite recursion in the parser when not used properly (which - causes the function stack created by the compiler to overflow). - If this possibility should be prevented then the eval() function - can be disabled completely by uncommenting this line. + Even though the maximum recursion level of the eval() function + is limited, it is still possible to write functions which never + reach this maximum recursion level but take enormous amounts of + time to evaluate (this can be undesirable eg. in web server-side + applications). + Uncommenting this line will disable the eval() function completely, + thus removing the danger of exploitation. + + Note that you can also disable eval() by specifying the + DISABLE_EVAL precompiler constant in your compiler (eg. + with -DDISABLE_EVAL in gcc). + + EVAL_MAX_REC_LEVEL : (Default 1000) + Sets the maximum recursion level allowed for eval(). SUPPORT_OPTIMIZER : (Default on) If you are not going to use the Optimize() method, you can comment @@ -95,6 +111,25 @@ compiler to see how this is done). well as making the binary a bit smaller. (Optimize() can still be called, but it will not do anything.) + You can also disable the optimizer by specifying the + NO_SUPPORT_OPTIMIZER precompiler constant in your compiler + (eg. with -DNO_SUPPORT_OPTIMIZER in gcc). + + FP_EPSILON : (Default 1e-14) + Epsilon value used in comparison operators. + If this line is commented out, then no epsilon will be used. + + FP_USE_THREAD_SAFE_EVAL : (Default off) + Define this precompiler constant to make Eval() thread-safe. + Refer to the thread safety section later in this document for + more information. + Note that defining this may make Eval() slightly slower. + + FP_USE_THREAD_SAFE_EVAL_WITH_ALLOCA : (Default off) + This like the previous, but makes Eval() use the alloca() function + (instead of std::vector). This will make it faster, but the alloca() + function is not standard and thus not supported by all compilers. + * Copying and assignment: ---------------------- @@ -173,6 +208,12 @@ bool AddConstant(const std::string& name, double value); is invalid, else true. +bool AddUnit(const std::string& name, double value); + + Add a new unit to the parser. Returns false if the name of the unit + is invalid, else true. + + bool AddFunction(const std::string& name, double (*functionPtr)(const double*), unsigned paramsAmount); @@ -212,7 +253,7 @@ int Parse(const std::string& Function, const std::string& Vars, "VAR" and "Var" are *different* variable names). Letters, digits and underscores can be used in variable names, but the name of a variable can't begin with a digit. Each variable name can appear only once in - the string. Function names are not legal variable names. + the 'Vars' string. Function names are not legal variable names. Using longer variable names causes no overhead whatsoever to the Eval() method, so it's completely safe to use variable names of any size. @@ -265,8 +306,8 @@ MISSING_PARENTH : "Missing ')'" EMPTY_PARENTH : "Empty parentheses" EXPECT_OPERATOR : "Syntax error: Operator expected" OUT_OF_MEMORY : "Not enough memory" -UNEXPECTED_ERROR : "An unexpected error ocurred. Please make a full bug " - "report to warp@iki.fi" +UNEXPECTED_ERROR : "An unexpected error occurred. Please make a full bug " + "report to the author" INVALID_VARS : "Syntax error in parameter 'Vars' given to " "FunctionParser::Parse()" ILL_PARAMS_AMOUNT : "Illegal number of parameters to function" @@ -308,6 +349,7 @@ int EvalError(void) const; 2: sqrt error (sqrt of a negative value) 3: log error (logarithm of a negative value) 4: trigonometric error (asin or acos of illegal value) + 5: maximum recursion level in eval() reached --------------------------------------------------------------------------- @@ -341,9 +383,9 @@ void Optimize(); Optimize() and checking EvalError(). If the destination application is not going to use this method, - the compiler constant SUPPORT_OPTIMIZER can be undefined at the beginning - of fparser.cc to make the library smaller (Optimize() can still be called, - but it will not do anything). + the compiler constant SUPPORT_OPTIMIZER can be undefined in fpconfig.hh + to make the library smaller (Optimize() can still be called, but it will + not do anything). (If you are interested in seeing how this method optimizes the opcode, you can call the PrintByteCode() method before and after the call to @@ -374,10 +416,53 @@ bool AddConstant(const std::string& name, double value); The return value will be false if the 'name' of the constant was illegal, else true. If the name was illegal, the method does nothing. - Example: parser.AddConstant("pi", 3.14159265); + Example: parser.AddConstant("pi", 3.1415926535897932); Now for example parser.Parse("x*pi", "x"); will be identical to the - call parser.Parse("x*3.14159265", "x"); + call parser.Parse("x*3.1415926535897932", "x"); + + +--------------------------------------------------------------------------- +bool AddUnit(const std::string& name, double value); +--------------------------------------------------------------------------- + + In some applications it is desirable to have units of measurement. + A typical example is an application which creates a page layout to be + printed. When printing, distances are usually measured in points + (defined by the resolution of the printer). However, it is often more + useful for the user to be able to specify measurements in other units + such as centimeters or inches. + + A unit is simply a value by which the preceding element is multiplied. + For example, if the printing has been set up to 300 DPI, one inch is + then 300 points (dots). Thus saying eg. "5in" is the same as saying + "5*300" or "1500" (assuming "in" has been added as a unit with the + value 300). + + Note that units are slightly different from a multiplication in + that they have a higher precedence than any other operator (except + parentheses). Thus for example "5/2in" is parsed as "5/(2*300)". + (If 5/2 inches is what one wants, it has to be written "(5/2)in".) + + You can use the AddUnit() method to add a new unit. The unit can then + be used after any element in the function (and will work as a multiplier + for that element). An element is a float literal, a constant, a variable, + a function or any expression in parentheses. When the element is not a + float literal nor an expression in parentheses, there has to naturally + be at least one whitespace between the element and the unit (eg. "x in"). + To change the value of a unit, call AddUnit() again with the same unit + name and the new value. + + Note that units are independent from variables, constants and functions. + This means that a unit can have the same name as those (there isn't any + syntactical problem in that), thus for example "x x" is a valid function + if "x" is a variable or constant and "x" is a unit, although for the sake + of clarity that is usually not recommended. + + Example: parser.AddUnit("in", 300); + + Now for example the function "5in" will be identical to "(5*300)". + Other usage examples include "x in", "3in+2", "pow(x,2)in", "(x+2)in". --------------------------------------------------------------------------- @@ -401,9 +486,9 @@ bool AddFunction(const std::string& name, must have the form: double functionName(const double* params); - - The number of parameters the function takes. NOTE: Currently this - value must be at least 1; the parser does not support functions which - take no parameters (this problem may be fixed in the future). + - The number of parameters the function takes. 0 is a valid value + in which case the function takes no parameters (such function + should simply ignore the double* it gets as a parameter). The return value will be false if the given name was invalid (either it did not follow the variable naming conventions, or the name was already @@ -424,10 +509,22 @@ bool AddFunction(const std::string& name, parser.Parse("2*sqr(x)", "x"); + An example of a useful function taking no parameters is a function + returning a random value. For example: + + double Rand(const double*) + { + return drand48(); + } + + parser.AddFunction("rand", Rand, 0); + + IMPORTANT NOTE: If you use the Optimize() method, it will assume that the user-given function has no side-effects, that is, it always returns the same value for the same parameters. The optimizer will optimize the function call away in some cases, making this assumption. + (The Rand() function given as example above is one such problematic case.) --------------------------------------------------------------------------- @@ -446,48 +543,57 @@ bool AddFunction(const std::string& name, FunctionParser&); want to use the parser A in the parser B, you must call A.Parse() before you can call B.AddFunction("name", A). - - The amount of parameters in the FunctionParser instance given as + - The amount of variables in the FunctionParser instance given as parameter must not change after it has been given to the AddFunction() - of another instance. Changing the number of parameters will result in + of another instance. Changing the number of variables will result in malfunction. - AddFunction() will fail (ie. return false) if a recursive loop is formed. The method specifically checks that no such loop is built. - - As with the other AddFunction(), the number of parameters taken by - the user-defined function must be at least 1 (this may be fixed in - the future). - Example: FunctionParser f1, f2; f1.Parse("x*x", "x"); f2.AddFunction("sqr", f1); + This version of the AddFunction() method can be useful to eg. chain + user-given functions. For example, ask the user for a function F1, + and then ask the user another function F2, but now the user can + call F1 in this second function if he wants (and so on with a third + function F3, where he can call F1 and F2, etc). ---------------------------------------------------------------------------- - - Example program: - -#include "fparser.hh" -#include -int main() -{ - FunctionParser fp; - - int ret = fp.Parse("x+y-1", "x,y"); - if(ret >= 0) - { - std::cerr << "At col " << ret << ": " << fp.ErrorMsg() << std::endl; - return 1; - } - double vals[] = { 4, 8 }; - std::cout << fp.Eval(vals) << std::endl; -} +============================================================================= + - About thread safety +============================================================================= + None of the member functions of the FunctionParser class are thread-safe. +Most prominently, the Eval() function is not thread-safe. (In other words, +the Eval() function of a single FunctionParser instance cannot be safely +called simultaneously by two threads.) + + There are ways to use this library in a thread-safe way, though. If each +thread uses its own FunctionParser instance, no problems will obviously +happen. Note, however, that if these instances need to be a copy of a given +FunctionParser instance (eg. one where the user has entered a function), +a deep copy of this instance has to be performed for each thread. By +default FunctionParser uses shallow-copying (copy-on-write), which means +that a simple assignment of copy construction will not copy the data itself. +To force a deep copy you can all the ForceDeepCopy() function on each of +the instances of each thread after the assignment or copying has been done. + + Another possibility is to compile the FunctionParser library so that +its Eval() function will be thread-safe. (This can be done by defining +the FP_USE_THREAD_SAFE_EVAL or the FP_USE_THREAD_SAFE_EVAL_WITH_ALLOCA +precompiler constant.) As long as only one thread calls the other functions +of FunctionParser, the other threads can safely call the Eval() of this +one instance. + Note, however, that compiling the library like this can make Eval() +slightly slower. (The alloca version, if supported by the compiler, +will not be as slow.) ============================================================================= @@ -499,18 +605,26 @@ int main() or functions using the following operators in this order of precedence: () expressions in parentheses first - -A unary minus + A unit a unit multiplier (if one has been added) A^B exponentiation (A raised to the power B) + -A unary minus + !A unary logical not (result is 1 if int(A) is 0, else 0) A*B A/B A%B multiplication, division and modulo A+B A-B addition and subtraction - A=B AB comparison between A and B (result is either 0 or 1) - A&B result is 1 if int(A) and int(B) differ from 0, else 0. - A|B result is 1 if int(A) or int(B) differ from 0, else 0. + A=B A!=B AB A>=B + comparison between A and B (result is either 0 or 1) + A&B result is 1 if int(A) and int(B) differ from 0, else 0 + A|B result is 1 if int(A) or int(B) differ from 0, else 0 Since the unary minus has higher precedence than any other operator, for example the following expression is valid: x*-y - Note that the '=' comparison can be inaccurate due to floating point - precision problems (eg. "sqrt(100)=10" probably returns 0, not 1). + + The comparison operators use an epsilon value, so expressions which may + differ in very least-significant digits should work correctly. For example, + "0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1 = 1" should always return 1, and + the same comparison done with ">" or "<" should always return 0. + (The epsilon value can be configured in the fpconfig.hh file.) + Without epsilon this comparison probably returns the wrong value. The class supports these functions: @@ -538,7 +652,7 @@ or functions using the following operators in this order of precedence: csc(A) : Cosecant of A (equivalent to 1/sin(A)). eval(...) : This a recursive call to the function to be evaluated. The number of parameters must be the same as the number of parameters - taken by the function. Usually called inside if() to avoid + taken by the function. Must be called inside if() to avoid infinite recursion. exp(A) : Exponential of A. Returns the value of e raised to the power A where e is the base of the natural logarithm, i.e. the @@ -581,13 +695,11 @@ or functions using the following operators in this order of precedence: extremely deep recursions should be avoided (eg. millions of nested recursive calls). - Also note that the if() function is the only place where making a recursive - call is safe. In any other place it will cause an infinite recursion (which - will make the program eventually run out of memory). If this is something - which should be avoided, it may be a good idea to disable the eval() - function completely. - The eval() function can be disabled with the DISABLE_EVAL precompiler - constant (see the beginning of fparser.cc). + Also note that even though the maximum recursion level of eval() is + limited, it is possible to write functions which never reach that level + but still take enormous amounts of time to evaluate. + This can sometimes be undesirable because it is prone to exploitation, + but you can disable the eval() function completely in the fpconfig.hh file. ============================================================================= @@ -676,7 +788,7 @@ make loops or anything. Usage license: ============================================================================= -Copyright © 2003 Juha Nieminen, Joel Yliluoma +Copyright © 2003-2005 Juha Nieminen, Joel Yliluoma This library is distributed under two distinct usage licenses depending on the software ("Software" below) which uses the Function Parser library diff --git a/deal.II/contrib/functionparser/fpconfig.h b/deal.II/contrib/functionparser/fpconfig.h new file mode 100644 index 0000000000..1eda03e54c --- /dev/null +++ b/deal.II/contrib/functionparser/fpconfig.h @@ -0,0 +1,67 @@ +//=============================== +// Function parser v2.83 by Warp +//=============================== + +// Configuration file +// ------------------ + +// NOTE: +// This file is for the internal use of the function parser only. +// You don't need to include this file in your source files, just +// include "fparser.hh". + +/* + Comment out the following line if your compiler supports the (non-standard) + asinh, acosh and atanh functions and you want them to be supported. If + you are not sure, just leave it (those function will then not be supported). +*/ +#define NO_ASINH + + +/* + Uncomment the following line to disable the eval() function if it could + be too dangerous in the target application. + Note that even though the maximum recursion level of eval() is limited, + it is still possible to write functions using it which take enormous + amounts of time to evaluate even though the maximum recursion is never + reached. This may be undesirable in some applications. +*/ +//#define DISABLE_EVAL + + +/* + Maximum recursion level for eval() calls: +*/ +#define EVAL_MAX_REC_LEVEL 1000 + + +/* + Comment out the following lines out if you are not going to use the + optimizer and want a slightly smaller library. The Optimize() method + can still be called, but it will not do anything. + If you are unsure, just leave it. It won't slow down the other parts of + the library. +*/ +#ifndef NO_SUPPORT_OPTIMIZER +#define SUPPORT_OPTIMIZER +#endif + + +/* + Epsilon value used with the comparison operators (must be non-negative): + (Comment it out if you don't want to use epsilon in comparisons. Might + lead to marginally faster evaluation of the comparison operators, but + can introduce inaccuracies in comparisons.) +*/ +#define FP_EPSILON 1e-14 + + +/* + No member function of FunctionParser is thread-safe. Most prominently, + Eval() is not thread-safe. By uncommenting one of these lines the Eval() + function can be made thread-safe at the cost of a possible small overhead. + The second version requires that the compiler supports the alloca() function, + which is not standard, but is faster. + */ +//#define FP_USE_THREAD_SAFE_EVAL +//#define FP_USE_THREAD_SAFE_EVAL_WITH_ALLOCA diff --git a/deal.II/contrib/functionparser/fpoptimizer.cc b/deal.II/contrib/functionparser/fpoptimizer.cc new file mode 100644 index 0000000000..66207296c4 --- /dev/null +++ b/deal.II/contrib/functionparser/fpoptimizer.cc @@ -0,0 +1,2094 @@ +//============================================ +// Function parser v2.82 optimizer by Bisqwit +//============================================ + +/* + NOTE! + ---- + Everything that goes into the #ifndef SUPPORT_OPTIMIZER part + (ie. when SUPPORT_OPTIMIZER is not defined) should be put in + the end of fparser.cc file, not in this file. + + Everything in this file should be inside the #ifdef SUPPORT_OPTIMIZER + block (except the #include "fpconfig.h" line). +*/ + + +#include "fpconfig.h" + +#ifdef SUPPORT_OPTIMIZER + +#include "fparser.h" +#include "fptypes.h" +using namespace FUNCTIONPARSERTYPES; + +#include +#include +#include +#include +using namespace std; + +#ifndef M_PI +#define M_PI 3.1415926535897932384626433832795 +#endif + +#define CONSTANT_E 2.71828182845904509080 // exp(1) +#define CONSTANT_PI M_PI // atan2(0,-1) +#define CONSTANT_L10 2.30258509299404590109 // log(10) +#define CONSTANT_L10I 0.43429448190325176116 // 1/log(10) +#define CONSTANT_L10E CONSTANT_L10I // log10(e) +#define CONSTANT_L10EI CONSTANT_L10 // 1/log10(e) +#define CONSTANT_DR (180.0 / M_PI) // 180/pi +#define CONSTANT_RD (M_PI / 180.0) // pi/180 + +// Debugging in optimizer? Not recommended for production use. +//#define TREE_DEBUG +#ifdef TREE_DEBUG +// Debug verbosely? +#include +#endif + +namespace { +inline double Min(double d1, double d2) +{ + return d1d2 ? d1 : d2; +} + +class compres +{ + // states: 0=false, 1=true, 2=unknown +public: + compres(bool b) : state(b) {} + compres(char v) : state(v) {} + // is it? + operator bool() const { return state != 0; } + // is it not? + bool operator! () const { return state != 1; } + bool operator==(bool b) const { return state != !b; } + bool operator!=(bool b) const { return state != b; } +private: + char state; +}; + +const compres maybe = (char)2; + +struct CodeTree; + +class SubTree +{ + CodeTree *tree; + bool sign; // Only possible when parent is cAdd or cMul + + inline void flipsign() { sign = !sign; } +public: + SubTree(); + SubTree(double value); + SubTree(const SubTree &b); + SubTree(const CodeTree &b); + + ~SubTree(); + const SubTree &operator= (const SubTree &b); + const SubTree &operator= (const CodeTree &b); + + bool getsign() const { return sign; } + + const CodeTree* operator-> () const { return tree; } + const CodeTree& operator* () const { return *tree; } + struct CodeTree* operator-> () { return tree; } + struct CodeTree& operator* () { return *tree; } + + bool operator< (const SubTree& b) const; + bool operator== (const SubTree& b) const; + void Negate(); // Note: Parent must be cAdd + void Invert(); // Note: Parent must be cMul + + void CheckConstNeg(); + void CheckConstInv(); +}; + +bool IsNegate(const SubTree &p1, const SubTree &p2); +bool IsInverse(const SubTree &p1, const SubTree &p2); + +typedef list paramlist; + +#ifdef TREE_DEBUG +struct CodeTreeData; +std::ostream& operator << (std::ostream& str, const CodeTreeData& tree); +std::ostream& operator << (std::ostream& str, const CodeTree& tree); +#endif + +struct CodeTreeData +{ + paramlist args; + +private: + unsigned op; // Operation + double value; // In case of cImmed + unsigned var; // In case of cVar + unsigned funcno; // In case of cFCall, cPCall + +public: + CodeTreeData() : op(cAdd) {} + ~CodeTreeData() {} + + void SetOp(unsigned newop) { op=newop; } + void SetFuncNo(unsigned newno) { funcno=newno; } + + inline unsigned GetOp() const { return op; } + inline double GetImmed() const { return value; } + inline unsigned GetVar() const { return var; } + inline unsigned GetFuncNo() const { return funcno; } + + bool IsFunc() const { return op == cFCall || op == cPCall; } + bool IsImmed() const { return op == cImmed; } + bool IsVar() const { return op == cVar; } + + void AddParam(const SubTree &p) + { + args.push_back(p); + } + void SetVar(unsigned v) + { + args.clear(); + op = cVar; + var = v; + } + void SetImmed(double v) + { + args.clear(); + op = cImmed; + value = orig = v; + inverted = negated = false; + } + void NegateImmed() + { + negated = !negated; + UpdateValue(); + } + void InvertImmed() + { + inverted = !inverted; + UpdateValue(); + } + + bool IsOriginal() const { return !(IsInverted() || IsNegated()); } + bool IsInverted() const { return inverted; } + bool IsNegated() const { return negated; } + bool IsInvertedOriginal() const { return IsInverted() && !IsNegated(); } + bool IsNegatedOriginal() const { return !IsInverted() && IsNegated(); } + +private: + void UpdateValue() + { + value = orig; + if(IsInverted()) { value = 1.0 / value; + // FIXME: potential divide by zero. + } + if(IsNegated()) value = -value; + } + + double orig; + bool inverted; + bool negated; +protected: + // Ensure we don't accidentally copy this + void operator=(const CodeTreeData &b); +}; + + +class CodeTreeDataPtr +{ + typedef pair p_t; + typedef p_t* pp; + mutable pp p; + + void Alloc() const { ++p->second; } + void Dealloc() const { if(!--p->second) delete p; p = 0; } + + void PrepareForWrite() + { + // We're ready if we're the only owner. + if(p->second == 1) return; + + // Then make a clone. + p_t *newtree = new p_t(p->first, 1); + // Forget the old + Dealloc(); + // Keep the new + p = newtree; + } + +public: + CodeTreeDataPtr() : p(new p_t) { p->second = 1; } + CodeTreeDataPtr(const CodeTreeDataPtr &b): p(b.p) { Alloc(); } + ~CodeTreeDataPtr() { Dealloc(); } + const CodeTreeDataPtr &operator= (const CodeTreeDataPtr &b) + { + b.Alloc(); + Dealloc(); + p = b.p; + return *this; + } + const CodeTreeData *operator-> () const { return &p->first; } + const CodeTreeData &operator* () const { return p->first; } + CodeTreeData *operator-> () { PrepareForWrite(); return &p->first; } + CodeTreeData &operator* () { PrepareForWrite(); return p->first; } + + void Shock(); +}; + + +#define CHECKCONSTNEG(item, op) \ + ((op)==cMul) \ + ? (item).CheckConstInv() \ + : (item).CheckConstNeg() + +struct CodeTree +{ + CodeTreeDataPtr data; + +private: + typedef paramlist::iterator pit; + typedef paramlist::const_iterator pcit; + + /* + template inline void chk() const + { + } + */ + +public: + const pcit GetBegin() const { return data->args.begin(); } + const pcit GetEnd() const { return data->args.end(); } + const pit GetBegin() { return data->args.begin(); } + const pit GetEnd() { return data->args.end(); } + const SubTree& getp0() const { /*chk<1>();*/pcit tmp=GetBegin(); return *tmp; } + const SubTree& getp1() const { /*chk<2>();*/pcit tmp=GetBegin(); ++tmp; return *tmp; } + const SubTree& getp2() const { /*chk<3>();*/pcit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; } + unsigned GetArgCount() const { return data->args.size(); } + void Erase(const pit p) { data->args.erase(p); } + + SubTree& getp0() { /*chk<1>();*/pit tmp=GetBegin(); return *tmp; } + SubTree& getp1() { /*chk<2>();*/pit tmp=GetBegin(); ++tmp; return *tmp; } + SubTree& getp2() { /*chk<3>();*/pit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; } + + // set + void SetImmed(double v) { data->SetImmed(v); } + void SetOp(unsigned op) { data->SetOp(op); } + void SetVar(unsigned v) { data->SetVar(v); } + // get + double GetImmed() const { return data->GetImmed(); } + unsigned GetVar() const { return data->GetVar(); } + unsigned GetOp() const { return data->GetOp(); } + // test + bool IsImmed() const { return data->IsImmed(); } + bool IsVar() const { return data->IsVar(); } + // act + void AddParam(const SubTree &p) { data->AddParam(p); } + void NegateImmed() { data->NegateImmed(); } // don't use when op!=cImmed + void InvertImmed() { data->InvertImmed(); } // don't use when op!=cImmed + + compres NonZero() const { if(!IsImmed()) return maybe; + return GetImmed() != 0.0; } + compres IsPositive() const { if(!IsImmed()) return maybe; + return GetImmed() > 0.0; } + +private: + struct ConstList + { + double voidvalue; + list cp; + double value; + unsigned size() const { return cp.size(); } + }; + struct ConstList BuildConstList(); + void KillConst(const ConstList &cl) + { + for(list::const_iterator i=cl.cp.begin(); i!=cl.cp.end(); ++i) + Erase(*i); + } + void FinishConst(const ConstList &cl, bool has_voidvalue) + { + if(has_voidvalue) + { + if(cl.value != cl.voidvalue) + { + if(cl.size() == 1) return; // nothing to do + AddParam(cl.value); // add the resulting constant + } + } + else + { + if(cl.size() <= 0) return; // nothing to do + AddParam(cl.value); // add the resulting constant + } + KillConst(cl); // remove all listed constants + } + +public: + CodeTree() {} + CodeTree(double v) { SetImmed(v); } + + CodeTree(unsigned op, const SubTree &p) + { + SetOp(op); + AddParam(p); + } + CodeTree(unsigned op, const SubTree &p1, const SubTree &p2) + { + SetOp(op); + AddParam(p1); + AddParam(p2); + } + + bool operator== (const CodeTree& b) const; + bool operator< (const CodeTree& b) const; + +private: + bool IsSortable() const + { + switch(GetOp()) + { + case cAdd: case cMul: + case cEqual: + case cAnd: case cOr: + case cMax: case cMin: + return true; + default: + return false; + } + } + void SortIfPossible() + { + if(IsSortable()) + { + data->args.sort(); + } + } + + void ReplaceWithConst(double value) + { + SetImmed(value); + + /* REMEMBER TO CALL CheckConstInv / CheckConstNeg + * FOR PARENT SubTree, OR MAYHEM HAPPENS + */ + } + + void ReplaceWith(const CodeTree &b) + { + // If b is child of *this, mayhem + // happens. So we first make a clone + // and then proceed with copy. + CodeTreeDataPtr tmp = b.data; + tmp.Shock(); + data = tmp; + } + + void ReplaceWith(unsigned op, const SubTree &p) + { + ReplaceWith(CodeTree(op, p)); + } + + void ReplaceWith(unsigned op, const SubTree &p1, const SubTree &p2) + { + ReplaceWith(CodeTree(op, p1, p2)); + } + + void OptimizeConflict() + { + // This optimization does this: x-x = 0, x/x = 1, a+b-a = b. + + if(GetOp() == cAdd || GetOp() == cMul) + { + Redo: + pit a, b; + for(a=GetBegin(); a!=GetEnd(); ++a) + { + for(b=a; ++b != GetEnd(); ) + { + const SubTree &p1 = *a; + const SubTree &p2 = *b; + + if(GetOp() == cMul ? IsInverse(p1,p2) + : IsNegate(p1,p2)) + { + // These parameters complement each others out + Erase(b); + Erase(a); + goto Redo; + } + } + } + } + OptimizeRedundant(); + } + + void OptimizeRedundant() + { + // This optimization does this: min()=0, max()=0, add()=0, mul()=1 + + if(!GetArgCount()) + { + if(GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax) + ReplaceWithConst(0); + else if(GetOp() == cMul) + ReplaceWithConst(1); + return; + } + + // And this: mul(x) = x, min(x) = x, max(x) = x, add(x) = x + + if(GetArgCount() == 1) + { + if(GetOp() == cMul || GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax) + if(!getp0().getsign()) + ReplaceWith(*getp0()); + } + + OptimizeDoubleNegations(); + } + + void OptimizeDoubleNegations() + { + if(GetOp() == cAdd) + { + // Eschew double negations + + // If any of the elements is !cAdd of cMul + // and has a numeric constant, negate + // the constant and negate sign. + + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + SubTree &pa = *a; + if(pa.getsign() + && pa->GetOp() == cMul) + { + // Before doing the check, ensure that + // we are not changing a x*1 into x*-1. + pa->OptimizeConstantMath1(); + + // OptimizeConstantMath1() may have changed + // the type of the expression (but not its sign), + // so verify that we're still talking about a cMul group. + if(pa->GetOp() == cMul) + { + CodeTree &p = *pa; + for(pit b=p.GetBegin(); + b!=p.GetEnd(); ++b) + { + SubTree &pb = *b; + if(pb->IsImmed()) + { + pb.Negate(); + pa.Negate(); + break; + } + } + } + } + } + } + + if(GetOp() == cMul) + { + // If any of the elements is !cMul of cPow + // and has a numeric exponent, negate + // the exponent and negate sign. + + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + SubTree &pa = *a; + if(pa.getsign() && pa->GetOp() == cPow) + { + CodeTree &p = *pa; + if(p.getp1()->IsImmed()) + { + // negate ok for pow when op=cImmed + p.getp1().Negate(); + pa.Negate(); + } + } + } + } + } + + void OptimizeConstantMath1() + { + // This optimization does three things: + // - For adding groups: + // Constants are added together. + // - For multiplying groups: + // Constants are multiplied together. + // - For function calls: + // If all parameters are constants, + // the call is replaced with constant value. + + // First, do this: + OptimizeAddMulFlat(); + +#ifdef TREE_DEBUG + cout << "BEFORE MATH1 :" << *this << endl; +#endif + switch(GetOp()) + { + case cAdd: + { + ConstList cl = BuildConstList(); + FinishConst(cl, true); + break; + } + case cMul: + { + ConstList cl = BuildConstList(); + + if(cl.value == 0.0) ReplaceWithConst(0.0); + else FinishConst(cl, true); + + break; + } + case cMin: + case cMax: + { + ConstList cl = BuildConstList(); + FinishConst(cl, false); // No "default" value + break; + } + + #define ConstantUnaryFun(token, fun) \ + case token: { const SubTree &p0 = getp0(); \ + if(p0->IsImmed()) ReplaceWithConst(fun(p0->GetImmed())); \ + break; } + #define ConstantBinaryFun(token, fun) \ + case token: { const SubTree &p0 = getp0(); \ + const SubTree &p1 = getp1(); \ + if(p0->IsImmed() && \ + p1->IsImmed()) ReplaceWithConst(fun(p0->GetImmed(), p1->GetImmed())); \ + break; } + + // FIXME: potential invalid parameters for functions + // can cause exceptions here + + ConstantUnaryFun(cAbs, fabs); + ConstantUnaryFun(cAcos, acos); + ConstantUnaryFun(cAsin, asin); + ConstantUnaryFun(cAtan, atan); + ConstantUnaryFun(cCeil, ceil); + ConstantUnaryFun(cCos, cos); + ConstantUnaryFun(cCosh, cosh); + ConstantUnaryFun(cFloor, floor); + ConstantUnaryFun(cLog, log); + ConstantUnaryFun(cSin, sin); + ConstantUnaryFun(cSinh, sinh); + ConstantUnaryFun(cTan, tan); + ConstantUnaryFun(cTanh, tanh); + ConstantBinaryFun(cAtan2, atan2); + ConstantBinaryFun(cMod, fmod); // not a func, but belongs here too + ConstantBinaryFun(cPow, pow); + + case cNeg: + case cSub: + case cDiv: + /* Unreached (nonexistent operator) + * TODO: internal error here? + */ + break; + + case cCot: + case cCsc: + case cSec: + case cDeg: + case cRad: + case cLog10: + case cSqrt: + case cExp: + /* Unreached (nonexistent function) + * TODO: internal error here? + */ + break; + } +#ifdef TREE_DEBUG + cout << "AFTER MATH1 A :" << *this << endl; +#endif + + OptimizeConflict(); + +#ifdef TREE_DEBUG + cout << "AFTER MATH1 B :" << *this << endl; +#endif + } + + void OptimizeAddMulFlat() + { + // This optimization flattens the topography of the tree. + // Examples: + // x + (y+z) = x+y+z + // x * (y/z) = x*y/z + // x / (y/z) = x/y*z + // min(x, min(y,z)) = min(x,y,z) + // max(x, max(y,z)) = max(x,y,z) + + if(GetOp() == cAdd + || GetOp() == cMul + || GetOp() == cMin + || GetOp() == cMax) + { + // If children are same type as parent add them here + for(pit b, a=GetBegin(); a!=GetEnd(); a=b) + { + const SubTree &pa = *a; b=a; ++b; + if(pa->GetOp() != GetOp()) continue; + + // Child is same type + for(pcit c=pa->GetBegin(); + c!=pa->GetEnd(); + ++c) + { + const SubTree &pb = *c; + if(pa.getsign()) + { + // +a -(+b +c) + // means b and c will be negated + + SubTree tmp = pb; + if(GetOp() == cMul) + tmp.Invert(); + else + tmp.Negate(); + AddParam(tmp); + } + else + AddParam(pb); + } + Erase(a); + + // Note: OptimizeConstantMath1() would be a good thing to call next. + } + } + } + + void OptimizeLinearCombine() + { + // This optimization does the following: + // + // x*x*x*x -> x^4 + // x+x+x+x -> x*4 + // x*x -> x^2 + // x/z/z -> + // + // min(x,x,y) -> min(x,y) + // max(x,x,y) -> max(x,y) + // max(x,x) -> max(x) + + // Remove conflicts first, so we don't have to worry about signs. + OptimizeConflict(); + + bool didchanges = false; + bool FactorIsImportant = true; + unsigned ReplacementOpcode = GetOp(); + + switch(GetOp()) + { + case cAdd: + FactorIsImportant = true; + ReplacementOpcode = cMul; + goto Redo; + case cMul: + FactorIsImportant = true; + ReplacementOpcode = cPow; + goto Redo; + case cMin: + FactorIsImportant = false; + goto Redo; + case cMax: + FactorIsImportant = false; + goto Redo; + default: break; + Redo: + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + + list poslist; + + for(pit b=a; ++b!=GetEnd(); ) + { + const SubTree &pb = *b; + if(*pa == *pb) + poslist.push_back(b); + } + + unsigned min = 2; + if(poslist.size() >= min) + { + SubTree arvo = pa; + bool negate = arvo.getsign(); + + double factor = poslist.size() + 1; + + if(negate) + { + arvo.Negate(); + factor = -factor; + } + + if(FactorIsImportant) + { + CodeTree tmp; + tmp.SetOp(ReplacementOpcode); + tmp.AddParam(arvo); + tmp.AddParam(factor); + + list::const_iterator j; + for(j=poslist.begin(); j!=poslist.end(); ++j) + Erase(*j); + poslist.clear(); + *a = tmp; + } + else + { + list::const_iterator j; + bool first=true; + for(j=poslist.begin(); j!=poslist.end(); ++j, first=false) + if(!first) + Erase(*j); + } + didchanges = true; + goto Redo; + } + } + } + + if(didchanges) + { + // As a result, there might be need for this: + OptimizeAddMulFlat(); + // And this: + OptimizeRedundant(); + } + } + + void OptimizeLogarithm() + { + /* + This is basic logarithm math: + pow(X,Y)/log(Y) = X + log(X)/log(Y) = logY(X) + log(X^Y) = log(X)*Y + log(X*Y) = log(X)+log(Y) + exp(log(X)*Y) = X^Y + + This function does these optimizations: + pow(const_E, log(x)) = x + pow(const_E, log(x)*y) = x^y + pow(10, log(x)*const_L10I*y) = x^y + pow(z, log(x)/log(z)*y) = x^y + + And this: + log(x^z) = z * log(x) + Which automatically causes these too: + log(pow(const_E, x)) = x + log(pow(y, x)) = x * log(y) + log(pow(pow(const_E, y), x)) = x*y + + And it does this too: + log(x) + log(y) + log(z) = log(x * y * z) + log(x * exp(y)) = log(x) + y + + */ + + // Must be already in exponential form. + + // Optimize exponents before doing something. + OptimizeExponents(); + + if(GetOp() == cLog) + { + // We should have one parameter for log() function. + // If we don't, we're screwed. + + const SubTree &p = getp0(); + + if(p->GetOp() == cPow) + { + // Found log(x^y) + SubTree p0 = p->getp0(); // x + SubTree p1 = p->getp1(); // y + + // Build the new logarithm. + CodeTree tmp(GetOp(), p0); // log(x) + + // Become log(x) * y + ReplaceWith(cMul, tmp, p1); + } + else if(p->GetOp() == cMul) + { + // Redefine &p nonconst + SubTree &p = getp0(); + + p->OptimizeAddMulFlat(); + p->OptimizeExponents(); + + p.CheckConstInv(); + + list adds; + + for(pit b, a = p->GetBegin(); + a != p->GetEnd(); a=b) + { + SubTree &pa = *a; b=a; ++b; + if(pa->GetOp() == cPow + && pa->getp0()->IsImmed() + && pa->getp0()->GetImmed() == CONSTANT_E) + { + adds.push_back(pa->getp1()); + p->Erase(a); + continue; + } + } + if(adds.size()) + { + CodeTree tmp(cAdd, *this); + + list::const_iterator i; + for(i=adds.begin(); i!=adds.end(); ++i) + tmp.AddParam(*i); + + ReplaceWith(tmp); + } + } + } + if(GetOp() == cAdd) + { + // Check which ones are logs. + list poslist; + + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + if(pa->GetOp() == cLog) + poslist.push_back(a); + } + + if(poslist.size() >= 2) + { + CodeTree tmp(cMul, 1.0); // eek + + list::const_iterator j; + for(j=poslist.begin(); j!=poslist.end(); ++j) + { + const SubTree &pb = **j; + // Take all of its children + for(pcit b=pb->GetBegin(); + b!=pb->GetEnd(); + ++b) + { + SubTree tmp2 = *b; + if(pb.getsign()) tmp2.Negate(); + tmp.AddParam(tmp2); + } + Erase(*j); + } + poslist.clear(); + + AddParam(CodeTree(cLog, tmp)); + } + // Done, hopefully + } + if(GetOp() == cPow) + { + const SubTree &p0 = getp0(); + SubTree &p1 = getp1(); + + if(p0->IsImmed() && p0->GetImmed() == CONSTANT_E + && p1->GetOp() == cLog) + { + // pow(const_E, log(x)) = x + ReplaceWith(*(p1->getp0())); + } + else if(p1->GetOp() == cMul) + { + //bool didsomething = true; + + pit poslogpos; bool foundposlog = false; + pit neglogpos; bool foundneglog = false; + + ConstList cl = p1->BuildConstList(); + + for(pit a=p1->GetBegin(); a!=p1->GetEnd(); ++a) + { + const SubTree &pa = *a; + if(pa->GetOp() == cLog) + { + if(!pa.getsign()) + { + foundposlog = true; + poslogpos = a; + } + else if(*p0 == *(pa->getp0())) + { + foundneglog = true; + neglogpos = a; + } + } + } + + if(p0->IsImmed() + && p0->GetImmed() == 10.0 + && cl.value == CONSTANT_L10I + && foundposlog) + { + SubTree base = (*poslogpos)->getp0(); + p1->KillConst(cl); + p1->Erase(poslogpos); + p1->OptimizeRedundant(); + SubTree mul = p1; + + ReplaceWith(cPow, base, mul); + + // FIXME: what optimizations should be done now? + return; + } + + // Put back the constant + FinishConst(cl, true); + + if(p0->IsImmed() + && p0->GetImmed() == CONSTANT_E + && foundposlog) + { + SubTree base = (*poslogpos)->getp0(); + p1->Erase(poslogpos); + + p1->OptimizeRedundant(); + SubTree mul = p1; + + ReplaceWith(cPow, base, mul); + + // FIXME: what optimizations should be done now? + return; + } + + if(foundposlog + && foundneglog + && *((*neglogpos)->getp0()) == *p0) + { + SubTree base = (*poslogpos)->getp0(); + p1->Erase(poslogpos); + p1->Erase(neglogpos); + + p1->OptimizeRedundant(); + SubTree mul = p1; + + ReplaceWith(cPow, base, mul); + + // FIXME: what optimizations should be done now? + return; + } + } + } + } + + void OptimizeFunctionCalls() + { + /* Goals: sin(asin(x)) = x + * cos(acos(x)) = x + * tan(atan(x)) = x + * NOTE: + * Do NOT do these: + * asin(sin(x)) + * acos(cos(x)) + * atan(tan(x)) + * Because someone might want to wrap the angle. + */ + // FIXME: TODO + } + + void OptimizePowMulAdd() + { + // x^3 * x -> x^4 + // x*3 + x -> x*4 + // FIXME: Do those + + // x^1 -> x + if(GetOp() == cPow) + { + const SubTree &base = getp0(); + const SubTree &exponent = getp1(); + + if(exponent->IsImmed()) + { + if(exponent->GetImmed() == 1.0) + ReplaceWith(*base); + else if(exponent->GetImmed() == 0.0 + && base->NonZero()) + ReplaceWithConst(1.0); + } + } + } + + void OptimizeExponents() + { + /* Goals: + * (x^y)^z -> x^(y*z) + * x^y * x^z -> x^(y+z) + */ + // First move to exponential form. + OptimizeLinearCombine(); + + bool didchanges = false; + + Redo: + if(GetOp() == cPow) + { + // (x^y)^z -> x^(y*z) + + const SubTree &p0 = getp0(); + const SubTree &p1 = getp1(); + if(p0->GetOp() == cPow) + { + CodeTree tmp(cMul, p0->getp1(), p1); + tmp.Optimize(); + + ReplaceWith(cPow, p0->getp0(), tmp); + + didchanges = true; + goto Redo; + } + } + if(GetOp() == cMul) + { + // x^y * x^z -> x^(y+z) + + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + + if(pa->GetOp() != cPow) continue; + + list poslist; + + for(pit b=a; ++b != GetEnd(); ) + { + const SubTree &pb = *b; + if(pb->GetOp() == cPow + && *(pa->getp0()) + == *(pb->getp0())) + { + poslist.push_back(b); + } + } + + if(poslist.size() >= 1) + { + poslist.push_back(a); + + CodeTree base = *(pa->getp0()); + + CodeTree exponent(cAdd, 0.0); //eek + + // Collect all exponents to cAdd + list::const_iterator i; + for(i=poslist.begin(); i!=poslist.end(); ++i) + { + const SubTree &pb = **i; + + SubTree tmp2 = pb->getp1(); + if(pb.getsign()) tmp2.Invert(); + + exponent.AddParam(tmp2); + } + + exponent.Optimize(); + + CodeTree result(cPow, base, exponent); + + for(i=poslist.begin(); i!=poslist.end(); ++i) + Erase(*i); + poslist.clear(); + + AddParam(result); // We're cMul, remember + + didchanges = true; + goto Redo; + } + } + } + + OptimizePowMulAdd(); + + if(didchanges) + { + // As a result, there might be need for this: + OptimizeConflict(); + } + } + + void OptimizeLinearExplode() + { + // x^2 -> x*x + // But only if x is just a simple thing + + // Won't work on anything else. + if(GetOp() != cPow) return; + + // TODO TODO TODO + } + + void OptimizePascal() + { +#if 0 // Too big, too specific, etc + + // Won't work on anything else. + if(GetOp() != cAdd) return; + + // Must be done after OptimizeLinearCombine(); + + // Don't need pascal triangle + // Coefficient for x^a * y^b * z^c = 3! / (a! * b! * c!) + + // We are greedy and want other than just binomials + // FIXME + + // note: partial ones are also nice + // x*x + x*y + y*y + // = (x+y)^2 - x*y + // + // x x * x y * + y y * + + // -> x y + dup * x y * - +#endif + } + +public: + + void Optimize(); + + void Assemble(vector &byteCode, + vector &immed, + size_t& stacktop_cur, + size_t& stacktop_max) const; + + void FinalOptimize() + { + // First optimize each parameter. + for(pit a=GetBegin(); a!=GetEnd(); ++a) + (*a)->FinalOptimize(); + + /* These things are to be done: + * + * x * CONSTANT_DR -> cDeg(x) + * x * CONSTANT_RD -> cRad(x) + * pow(x, 0.5) -> sqrt(x) + * log(x) * CONSTANT_L10I -> log10(x) + * pow(CONSTANT_E, x) -> exp(x) + * inv(sin(x)) -> csc(x) + * inv(cos(x)) -> sec(x) + * inv(tan(x)) -> cot(x) + */ + + + if(GetOp() == cPow) + { + const SubTree &p0 = getp0(); + const SubTree &p1 = getp1(); + if(p0->GetOp() == cImmed + && p0->GetImmed() == CONSTANT_E) + { + ReplaceWith(cExp, p1); + } + else if(p1->GetOp() == cImmed + && p1->GetImmed() == 0.5) + { + ReplaceWith(cSqrt, p0); + } + } + if(GetOp() == cMul) + { + if(GetArgCount() == 1 && getp0().getsign()) + { + /***/if(getp0()->GetOp() == cSin)ReplaceWith(cCsc, getp0()->getp0()); + else if(getp0()->GetOp() == cCos)ReplaceWith(cSec, getp0()->getp0()); + else if(getp0()->GetOp() == cTan)ReplaceWith(cCot, getp0()->getp0()); + } + } + // Separate "if", because op may have just changed + if(GetOp() == cMul) + { + /* NOTE: Not gcov'd yet */ + + CodeTree *found_log = 0; + + ConstList cl = BuildConstList(); + + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + SubTree &pa = *a; + if(pa->GetOp() == cLog && !pa.getsign()) + found_log = &*pa; + } + if(cl.value == CONSTANT_L10I && found_log) + { + // Change the log() to log10() + found_log->SetOp(cLog10); + // And forget the constant + KillConst(cl); + } + else if(cl.value == CONSTANT_DR) + { + KillConst(cl); + OptimizeRedundant(); +#ifdef TREE_DEBUG + cout << "PRE_REP :" << (*this) << endl; +#endif + ReplaceWith(cDeg, *this); +#ifdef TREE_DEBUG + cout << "POST_REP :" << (*this) << endl; +#endif + } + else if(cl.value == CONSTANT_RD) + { + KillConst(cl); + OptimizeRedundant(); + ReplaceWith(cRad, *this); + } + else FinishConst(cl, true); + } + + SortIfPossible(); + } +}; + +#ifdef TREE_DEBUG +std::ostream& operator << (std::ostream& str, const CodeTree& tree) +{ + const CodeTreeData& data = *tree.data; + switch( (FUNCTIONPARSERTYPES::OPCODE) data.GetOp()) + { + case cImmed: str << data.GetImmed(); return str; + case cVar: str << "Var" << data.GetVar(); return str; + case cFCall: str << "FCall(Func" << data.GetFuncNo() << ")"; break; + case cPCall: str << "PCall(Func" << data.GetFuncNo() << ")"; break; + + case cAbs: str << "cAbs"; break; + case cAcos: str << "cAcos"; break; +#ifndef NO_ASINH + case cAcosh: str << "cAcosh"; break; +#endif + case cAsin: str << "cAsin"; break; +#ifndef NO_ASINH + case cAsinh: str << "cAsinh"; break; +#endif + case cAtan: str << "cAtan"; break; + case cAtan2: str << "cAtan2"; break; +#ifndef NO_ASINH + case cAtanh: str << "cAtanh"; break; +#endif + case cCeil: str << "cCeil"; break; + case cCos: str << "cCos"; break; + case cCosh: str << "cCosh"; break; + case cCot: str << "cCot"; break; + case cCsc: str << "cCsc"; break; +#ifndef DISABLE_EVAL + case cEval: str << "cEval"; break; +#endif + case cExp: str << "cExp"; break; + case cFloor: str << "cFloor"; break; + case cIf: str << "cIf"; break; + case cInt: str << "cInt"; break; + case cLog: str << "cLog"; break; + case cLog10: str << "cLog10"; break; + case cMax: str << "cMax"; break; + case cMin: str << "cMin"; break; + case cSec: str << "cSec"; break; + case cSin: str << "cSin"; break; + case cSinh: str << "cSinh"; break; + case cSqrt: str << "cSqrt"; break; + case cTan: str << "cTan"; break; + case cTanh: str << "cTanh"; break; + + // These do not need any ordering: + case cJump: str << "cJump"; break; + case cNeg: str << "cNeg"; break; + case cAdd: str << "cAdd"; break; + case cSub: str << "cSub"; break; + case cMul: str << "cMul"; break; + case cDiv: str << "cDiv"; break; + case cMod: str << "cMod"; break; + case cPow: str << "cPow"; break; + case cEqual: str << "cEqual"; break; + case cNEqual: str << "cNEqual"; break; + case cLess: str << "cLess"; break; + case cLessOrEq: str << "cLessOrEq"; break; + case cGreater: str << "cGreater"; break; + case cGreaterOrEq: str << "cGreaterOrEq"; break; + case cNot: str << "cNot"; break; + case cAnd: str << "cAnd"; break; + case cOr: str << "cOr"; break; + case cDeg: str << "cDeg"; break; + case cRad: str << "cRad"; break; + case cDup: str << "cDup"; break; + case cInv: str << "cInv"; break; + case VarBegin: str << "VarBegin"; break; + } + str << '('; + + bool first = true; + for(paramlist::const_iterator + i = data.args.begin(); i != data.args.end(); ++i) + { + if(first) first=false; else str << ", "; + const SubTree& sub = *i; + if(sub.getsign()) str << '!'; + str << *sub; + } + str << ')'; + + return str; +} +#endif + +void CodeTreeDataPtr::Shock() +{ + /* + PrepareForWrite(); + paramlist &p2 = (*this)->args; + for(paramlist::iterator i=p2.begin(); i!=p2.end(); ++i) + { + (*i)->data.Shock(); + } + */ +} + +CodeTree::ConstList CodeTree::BuildConstList() +{ + ConstList result; + result.value = + result.voidvalue = GetOp()==cMul ? 1.0 : 0.0; + + bool first = true; + + list &cp = result.cp; + for(pit b, a=GetBegin(); a!=GetEnd(); a=b) + { + SubTree &pa = *a; b=a; ++b; + if(!pa->IsImmed()) continue; + + double thisvalue = pa->GetImmed(); + if(thisvalue == result.voidvalue) + { + // This value is no good, forget it + Erase(a); + continue; + } + switch(GetOp()) + { + case cMul: + result.value *= thisvalue; + break; + case cAdd: + result.value += thisvalue; + break; + case cMin: + if(first || thisvalue < result.value) result.value = thisvalue; + break; + case cMax: + if(first || thisvalue > result.value) result.value = thisvalue; + break; + default: break; /* Unreached */ + } + cp.push_back(a); + first = false; + } + if(GetOp() == cMul) + { + /* + If one of the values is -1 and it's not the only value, + then some other value will be negated. + */ + for(bool done=false; cp.size() > 1 && !done; ) + { + done = true; + for(list::iterator b,a=cp.begin(); a!=cp.end(); a=b) + { + b=a; ++b; + if((**a)->GetImmed() == -1.0) + { + Erase(*a); + cp.erase(a); + + // take randomly something + (**cp.begin())->data->NegateImmed(); + if(cp.size() < 2)break; + done = false; + } + } + } + } + return result; +} + +void CodeTree::Assemble + (vector &byteCode, + vector &immed, + size_t& stacktop_cur, + size_t& stacktop_max) const +{ + #define AddCmd(op) byteCode.push_back((op)) + #define AddConst(v) do { \ + byteCode.push_back(cImmed); \ + immed.push_back((v)); \ + } while(0) + #define SimuPush(n) stacktop_cur += (n) + #define SimuPop(n) do { \ + if(stacktop_cur > stacktop_max) stacktop_max = stacktop_cur; \ + stacktop_cur -= (n); \ + } while(0) + + if(IsVar()) + { + SimuPush(1); + AddCmd(GetVar()); + return; + } + if(IsImmed()) + { + SimuPush(1); + AddConst(GetImmed()); + return; + } + + switch(GetOp()) + { + case cAdd: + case cMul: + case cMin: + case cMax: + { + unsigned opcount = 0; + for(pcit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + + if(opcount < 2) ++opcount; + + bool pnega = pa.getsign(); + + bool done = false; + if(pa->IsImmed()) + { + if(GetOp() == cMul + && pa->data->IsInverted() + && (pnega || opcount==2) + ) + { + CodeTree tmp = *pa; + tmp.data->InvertImmed(); + tmp.Assemble(byteCode, immed, + stacktop_cur, stacktop_max); + pnega = !pnega; + done = true; + } + else if(GetOp() == cAdd + && (pa->data->IsNegatedOriginal() + // || pa->GetImmed() < 0 + ) + && (pnega || opcount==2) + ) + { + CodeTree tmp = *pa; + tmp.data->NegateImmed(); + tmp.Assemble(byteCode, immed, + stacktop_cur, stacktop_max); + pnega = !pnega; + done = true; + } + } + if(!done) + pa->Assemble(byteCode, immed, stacktop_cur, stacktop_max); + + if(opcount == 2) + { + unsigned tmpop = GetOp(); + if(pnega) // negation in non-first operand + { + tmpop = (tmpop == cMul) ? cDiv : cSub; + } + SimuPop(1); + AddCmd(tmpop); + } + else if(pnega) // negation in the first operand + { + if(GetOp() == cMul) AddCmd(cInv); + else AddCmd(cNeg); + } + } + break; + } + case cIf: + { + // If the parameter amount is != 3, we're screwed. + getp0()->Assemble(byteCode, immed, stacktop_cur, stacktop_max); // expression + SimuPop(1); + + unsigned ofs = byteCode.size(); + AddCmd(cIf); + AddCmd(0); // code index + AddCmd(0); // immed index + + getp1()->Assemble(byteCode, immed, stacktop_cur, stacktop_max); // true branch + SimuPop(1); + + byteCode[ofs+1] = byteCode.size()+2; + byteCode[ofs+2] = immed.size(); + + ofs = byteCode.size(); + AddCmd(cJump); + AddCmd(0); // code index + AddCmd(0); // immed index + + getp2()->Assemble(byteCode, immed, stacktop_cur, stacktop_max); // false branch + SimuPop(1); + + byteCode[ofs+1] = byteCode.size()-1; + byteCode[ofs+2] = immed.size(); + + SimuPush(1); + + break; + } + case cFCall: + { + // If the parameter count is invalid, we're screwed. + size_t was_stacktop = stacktop_cur; + for(pcit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + pa->Assemble(byteCode, immed, stacktop_cur, stacktop_max); + } + AddCmd(GetOp()); + AddCmd(data->GetFuncNo()); + SimuPop(stacktop_cur - was_stacktop - 1); + break; + } + case cPCall: + { + // If the parameter count is invalid, we're screwed. + size_t was_stacktop = stacktop_cur; + for(pcit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + pa->Assemble(byteCode, immed, stacktop_cur, stacktop_max); + } + AddCmd(GetOp()); + AddCmd(data->GetFuncNo()); + SimuPop(stacktop_cur - was_stacktop - 1); + break; + } + default: + { + // If the parameter count is invalid, we're screwed. + size_t was_stacktop = stacktop_cur; + for(pcit a=GetBegin(); a!=GetEnd(); ++a) + { + const SubTree &pa = *a; + pa->Assemble(byteCode, immed, stacktop_cur, stacktop_max); + } + AddCmd(GetOp()); + SimuPop(stacktop_cur - was_stacktop - 1); + break; + } + } +} + +void CodeTree::Optimize() +{ + // Phase: + // Phase 0: Do local optimizations. + // Phase 1: Optimize each. + // Phase 2: Do local optimizations again. + + for(unsigned phase=0; phase<=2; ++phase) + { + if(phase == 1) + { + // Optimize each parameter. + for(pit a=GetBegin(); a!=GetEnd(); ++a) + { + (*a)->Optimize(); + CHECKCONSTNEG(*a, GetOp()); + } + continue; + } + if(phase == 0 || phase == 2) + { + // Do local optimizations. + + OptimizeConstantMath1(); + OptimizeLogarithm(); + OptimizeFunctionCalls(); + OptimizeExponents(); + OptimizeLinearExplode(); + OptimizePascal(); + + /* Optimization paths: + + doublenegations= + redundant= * doublenegations + conflict= * redundant + addmulflat= + constantmath1= addmulflat * conflict + linearcombine= conflict * addmulflat¹ redundant¹ + powmuladd= + exponents= linearcombine * powmuladd conflict¹ + logarithm= exponents * + functioncalls= IDLE + linearexplode= IDLE + pascal= IDLE + + * = actions here + ¹ = only if made changes + */ + } + } +} + + +bool CodeTree::operator== (const CodeTree& b) const +{ + if(GetOp() != b.GetOp()) return false; + if(IsImmed()) if(GetImmed() != b.GetImmed()) return false; + if(IsVar()) if(GetVar() != b.GetVar()) return false; + if(data->IsFunc()) + if(data->GetFuncNo() != b.data->GetFuncNo()) return false; + return data->args == b.data->args; +} + +bool CodeTree::operator< (const CodeTree& b) const +{ + if(GetArgCount() != b.GetArgCount()) + return GetArgCount() > b.GetArgCount(); + + if(GetOp() != b.GetOp()) + { + // sort immeds last + if(IsImmed() != b.IsImmed()) return IsImmed() < b.IsImmed(); + + return GetOp() < b.GetOp(); + } + + if(IsImmed()) + { + if(GetImmed() != b.GetImmed()) return GetImmed() < b.GetImmed(); + } + if(IsVar() && GetVar() != b.GetVar()) + { + return GetVar() < b.GetVar(); + } + if(data->IsFunc() && data->GetFuncNo() != b.data->GetFuncNo()) + { + return data->GetFuncNo() < b.data->GetFuncNo(); + } + + pcit i = GetBegin(), j = b.GetBegin(); + for(; i != GetEnd(); ++i, ++j) + { + const SubTree &pa = *i, &pb = *j; + + if(!(pa == pb)) + return pa < pb; + } + return false; +} + + +bool IsNegate(const SubTree &p1, const SubTree &p2) /*const */ +{ + if(p1->IsImmed() && p2->IsImmed()) + { + return p1->GetImmed() == -p2->GetImmed(); + } + if(p1.getsign() == p2.getsign()) return false; + return *p1 == *p2; +} +bool IsInverse(const SubTree &p1, const SubTree &p2) /*const*/ +{ + if(p1->IsImmed() && p2->IsImmed()) + { + // FIXME: potential divide by zero. + return p1->GetImmed() == 1.0 / p2->GetImmed(); + } + if(p1.getsign() == p2.getsign()) return false; + return *p1 == *p2; +} + +SubTree::SubTree() : tree(new CodeTree), sign(false) +{ +} + +SubTree::SubTree(const SubTree &b) : tree(new CodeTree(*b.tree)), sign(b.sign) +{ +} + +#define SubTreeDecl(p1, p2) \ + SubTree::SubTree p1 : tree(new CodeTree p2), sign(false) { } + +SubTreeDecl( (const CodeTree &b), (b) ) +SubTreeDecl( (double value), (value) ) + +#undef SubTreeDecl + +SubTree::~SubTree() +{ + delete tree; tree=0; +} + +const SubTree &SubTree::operator= (const SubTree &b) +{ + sign = b.sign; + CodeTree *oldtree = tree; + tree = new CodeTree(*b.tree); + delete oldtree; + return *this; +} +const SubTree &SubTree::operator= (const CodeTree &b) +{ + sign = false; + CodeTree *oldtree = tree; + tree = new CodeTree(b); + delete oldtree; + return *this; +} + +bool SubTree::operator< (const SubTree& b) const +{ + if(getsign() != b.getsign()) return getsign() < b.getsign(); + return *tree < *b.tree; +} +bool SubTree::operator== (const SubTree& b) const +{ + return sign == b.sign && *tree == *b.tree; +} +void SubTree::Negate() // Note: Parent must be cAdd +{ + flipsign(); + CheckConstNeg(); +} +void SubTree::CheckConstNeg() +{ + if(tree->IsImmed() && getsign()) + { + tree->NegateImmed(); + sign = false; + } +} +void SubTree::Invert() // Note: Parent must be cMul +{ + flipsign(); + CheckConstInv(); +} +void SubTree::CheckConstInv() +{ + if(tree->IsImmed() && getsign()) + { + tree->InvertImmed(); + sign = false; + } +} + +}//namespace + +void FunctionParser::MakeTree(void *r) const +{ + // Dirty hack. Should be fixed. + CodeTree* result = static_cast(r); + + vector stack(1); + + #define GROW(n) do { \ + stacktop += n; \ + if(stack.size() <= stacktop) stack.resize(stacktop+1); \ + } while(0) + + #define EAT(n, opcode) do { \ + unsigned newstacktop = stacktop-(n); \ + if((n) == 0) GROW(1); \ + stack[stacktop].SetOp((opcode)); \ + for(unsigned a=0, b=(n); a labels; + + const unsigned* const ByteCode = data->ByteCode; + const unsigned ByteCodeSize = data->ByteCodeSize; + const double* const Immed = data->Immed; + + for(unsigned IP=0, DP=0; ; ++IP) + { + while(labels.size() > 0 + && *labels.begin() == IP) + { + // The "else" of an "if" ends here + EAT(3, cIf); + labels.erase(labels.begin()); + } + + if(IP >= ByteCodeSize) + { + break; + } + + unsigned opcode = ByteCode[IP]; + + if(opcode == cIf) + { + IP += 2; + } + else if(opcode == cJump) + { + labels.push_front(ByteCode[IP+1]+1); + IP += 2; + } + else if(opcode == cImmed) + { + ADDCONST(Immed[DP++]); + } + else if(opcode < VarBegin) + { + switch(opcode) + { + // Unary operators + case cNeg: + { + EAT(1, cAdd); // Unary minus is negative adding. + stack[stacktop-1].getp0().Negate(); + break; + } + // Binary operators + case cSub: + { + EAT(2, cAdd); // Minus is negative adding + stack[stacktop-1].getp1().Negate(); + break; + } + case cDiv: + { + EAT(2, cMul); // Divide is inverse multiply + stack[stacktop-1].getp1().Invert(); + break; + } + + // ADD ALL TWO PARAMETER NON-FUNCTIONS HERE + case cAdd: case cMul: + case cMod: case cPow: + case cEqual: case cLess: case cGreater: + case cNEqual: case cLessOrEq: case cGreaterOrEq: + case cAnd: case cOr: + EAT(2, opcode); + break; + + // ADD ALL UNARY NON-FUNCTIONS HERE + case cNot: + EAT(1, opcode); + break; + + case cFCall: + { + unsigned index = ByteCode[++IP]; + unsigned params = data->FuncPtrs[index].params; + EAT(params, opcode); + stack[stacktop-1].data->SetFuncNo(index); + break; + } + case cPCall: + { + unsigned index = ByteCode[++IP]; + unsigned params = + data->FuncParsers[index]->data->varAmount; + EAT(params, opcode); + stack[stacktop-1].data->SetFuncNo(index); + break; + } + + // Converted to cMul on fly + case cDeg: + ADDCONST(CONSTANT_DR); + EAT(2, cMul); + break; + + // Converted to cMul on fly + case cRad: + ADDCONST(CONSTANT_RD); + EAT(2, cMul); + break; + + // Functions + default: + { + //assert(opcode >= cAbs); + unsigned funcno = opcode-cAbs; + assert(funcno < sizeof(Functions)/sizeof(Functions[0])); + const FuncDefinition& func = Functions[funcno]; + + //const FuncDefinition& func = Functions[opcode-cAbs]; + + unsigned paramcount = func.params; +#ifndef DISABLE_EVAL + if(opcode == cEval) paramcount = data->varAmount; +#endif + if(opcode == cSqrt) + { + // Converted on fly: sqrt(x) = x^0.5 + opcode = cPow; + paramcount = 2; + ADDCONST(0.5); + } + if(opcode == cExp) + { + // Converted on fly: exp(x) = CONSTANT_E^x + + opcode = cPow; + paramcount = 2; + // reverse the parameters... kludgey + stack[stacktop] = stack[stacktop-1]; + stack[stacktop-1].SetImmed(CONSTANT_E); + GROW(1); + } + bool do_inv = false; + if(opcode == cCot) { do_inv = true; opcode = cTan; } + if(opcode == cCsc) { do_inv = true; opcode = cSin; } + if(opcode == cSec) { do_inv = true; opcode = cCos; } + + bool do_log10 = false; + if(opcode == cLog10) + { + // Converted on fly: log10(x) = log(x) * CONSTANT_L10I + opcode = cLog; + do_log10 = true; + } + EAT(paramcount, opcode); + if(do_log10) + { + ADDCONST(CONSTANT_L10I); + EAT(2, cMul); + } + if(do_inv) + { + // Unary cMul, inverted. No need for "1.0" + EAT(1, cMul); + stack[stacktop-1].getp0().Invert(); + } + break; + } + } + } + else + { + stack[stacktop].SetVar(opcode); + GROW(1); + } + } + + if(!stacktop) + { + // ERROR: Stack does not have any values! + return; + } + + --stacktop; // Ignore the last element, it is always nop (cAdd). + + if(stacktop > 0) + { + // ERROR: Stack has too many values! + return; + } + + // Okay, the tree is now stack[0] + *result = stack[0]; +} + +void FunctionParser::Optimize() +{ + if(data->isOptimized) return; + CopyOnWrite(); + + CodeTree tree; + MakeTree(&tree); + +#ifdef TREE_DEBUG + cout << "BEFORE OPT :" << tree << endl; +#endif + // Do all sorts of optimizations + tree.Optimize(); + +#ifdef TREE_DEBUG + cout << "BEFORE FINAL :" << tree << endl; +#endif + + // Last changes before assembly + tree.FinalOptimize(); + +#ifdef TREE_DEBUG + cout << "AFTER :" << tree << endl; +#endif + // Now rebuild from the tree. + + vector byteCode; + vector immed; + +#if 0 + byteCode.resize(Comp.ByteCodeSize); + for(unsigned a=0; aStackSize < stacktop_max) + { + delete[] data->Stack; + data->Stack = new double[stacktop_max]; + data->StackSize = stacktop_max; + } +#endif + + delete[] data->ByteCode; data->ByteCode = 0; + if((data->ByteCodeSize = byteCode.size()) > 0) + { + data->ByteCode = new unsigned[data->ByteCodeSize]; + for(unsigned a=0; aByteCode[a] = byteCode[a]; + } + + delete[] data->Immed; data->Immed = 0; + if((data->ImmedSize = immed.size()) > 0) + { + data->Immed = new double[data->ImmedSize]; + for(unsigned a=0; aImmed[a] = immed[a]; + } + + data->isOptimized = true; +} + + +#endif // #ifdef SUPPORT_OPTIMIZER diff --git a/deal.II/contrib/functionparser/fptypes.h b/deal.II/contrib/functionparser/fptypes.h new file mode 100644 index 0000000000..a3534912cb --- /dev/null +++ b/deal.II/contrib/functionparser/fptypes.h @@ -0,0 +1,117 @@ +//=============================== +// Function parser v2.83 by Warp +//=============================== + +// NOTE: +// This file contains only internal types for the function parser library. +// You don't need to include this file in your code. Include "fparser.hh" +// only. + + + +namespace FUNCTIONPARSERTYPES +{ +// The functions must be in alphabetical order: + enum OPCODE + { + cAbs, cAcos, +#ifndef NO_ASINH + cAcosh, +#endif + cAsin, +#ifndef NO_ASINH + cAsinh, +#endif + cAtan, + cAtan2, +#ifndef NO_ASINH + cAtanh, +#endif + cCeil, cCos, cCosh, cCot, cCsc, +#ifndef DISABLE_EVAL + cEval, +#endif + cExp, cFloor, cIf, cInt, cLog, cLog10, cMax, cMin, + cSec, cSin, cSinh, cSqrt, cTan, cTanh, + +// These do not need any ordering: + cImmed, cJump, + cNeg, cAdd, cSub, cMul, cDiv, cMod, cPow, + cEqual, cNEqual, cLess, cLessOrEq, cGreater, cGreaterOrEq, + cNot, cAnd, cOr, + + cDeg, cRad, + + cFCall, cPCall, + +#ifdef SUPPORT_OPTIMIZER + cVar, cDup, cInv, +#endif + + VarBegin + }; + + struct FuncDefinition + { + const char* name; + unsigned nameLength; + unsigned opcode; + unsigned params; + + // This is basically strcmp(), but taking 'nameLength' as string + // length (not ending '\0'): + bool operator<(const FuncDefinition& rhs) const + { + for(unsigned i = 0; i < nameLength; ++i) + { + if(i == rhs.nameLength) return false; + const char c1 = name[i], c2 = rhs.name[i]; + if(c1 < c2) return true; + if(c2 < c1) return false; + } + return nameLength < rhs.nameLength; + } + }; + + +// This list must be in alphabetical order: + const FuncDefinition Functions[]= + { + { "abs", 3, cAbs, 1 }, + { "acos", 4, cAcos, 1 }, +#ifndef NO_ASINH + { "acosh", 5, cAcosh, 1 }, +#endif + { "asin", 4, cAsin, 1 }, +#ifndef NO_ASINH + { "asinh", 5, cAsinh, 1 }, +#endif + { "atan", 4, cAtan, 1 }, + { "atan2", 5, cAtan2, 2 }, +#ifndef NO_ASINH + { "atanh", 5, cAtanh, 1 }, +#endif + { "ceil", 4, cCeil, 1 }, + { "cos", 3, cCos, 1 }, + { "cosh", 4, cCosh, 1 }, + { "cot", 3, cCot, 1 }, + { "csc", 3, cCsc, 1 }, +#ifndef DISABLE_EVAL + { "eval", 4, cEval, 0 }, +#endif + { "exp", 3, cExp, 1 }, + { "floor", 5, cFloor, 1 }, + { "if", 2, cIf, 0 }, + { "int", 3, cInt, 1 }, + { "log", 3, cLog, 1 }, + { "log10", 5, cLog10, 1 }, + { "max", 3, cMax, 2 }, + { "min", 3, cMin, 2 }, + { "sec", 3, cSec, 1 }, + { "sin", 3, cSin, 1 }, + { "sinh", 4, cSinh, 1 }, + { "sqrt", 4, cSqrt, 1 }, + { "tan", 3, cTan, 1 }, + { "tanh", 4, cTanh, 1 } + }; +} diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 58b5117653..acea0135ce 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -526,6 +526,13 @@ inconvenience this causes.

deal.II

    +
  1. +

    + Upgraded: The FunctionParser classes now use version 2.83. +
    + (Luca Heltai 2008/12/08) +

    +
  2. Fixed: The GridGenerator::laplace_transform would only do at most 1000 -- 2.39.5