]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add the function parser.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 7 Mar 2005 02:33:55 +0000 (02:33 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 7 Mar 2005 02:33:55 +0000 (02:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@10018 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/contrib/functionparser/Makefile [new file with mode: 0644]
deal.II/contrib/functionparser/README [new file with mode: 0644]
deal.II/contrib/functionparser/fparser.cc [new file with mode: 0644]
deal.II/contrib/functionparser/fparser.h [new file with mode: 0644]
deal.II/contrib/functionparser/fparser.txt [new file with mode: 0644]

diff --git a/deal.II/contrib/functionparser/Makefile b/deal.II/contrib/functionparser/Makefile
new file mode 100644 (file)
index 0000000..7196804
--- /dev/null
@@ -0,0 +1,28 @@
+# $Id$
+# Copyright (C) 2005 by Wolfgang Bangerth
+
+
+D = ../..
+
+include $D/common/Make.global_options
+
+ifneq ($(with-multithreading),no)
+  MT = MT
+else
+  MT = ==
+endif
+
+
+
+# rules
+$(LIBDIR)/contrib/functionparser/fparser.$(OBJEXT): fparser.h fparser.cc
+       @echo =====function parser==optimized==$(MT)== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c fparser.cc -o $@
+
+
+
+clean:
+       -rm -f *~ 
+
+
+.PHONY: clean
diff --git a/deal.II/contrib/functionparser/README b/deal.II/contrib/functionparser/README
new file mode 100644 (file)
index 0000000..989892b
--- /dev/null
@@ -0,0 +1,7 @@
+The files in this directory are taken from the fparser distribution, version
+2.7. The code is subject to copyright by its author, Warp. The license is
+given in the file fparser.txt.
+
+In order to update the files to a newer version, just replace the files, but
+make sure you rename the fparser.hh file by fparser.h, to keep with the usual
+naming scheme in deal.II.
diff --git a/deal.II/contrib/functionparser/fparser.cc b/deal.II/contrib/functionparser/fparser.cc
new file mode 100644 (file)
index 0000000..44be603
--- /dev/null
@@ -0,0 +1,3179 @@
+//==============================
+// Function parser v2.7 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.
+#define SUPPORT_OPTIMIZER
+
+
+//============================================================================
+
+#include "fparser.h"
+
+#include <cstdlib>
+#include <cstring>
+#include <cctype>
+#include <cmath>
+
+using namespace std;
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626433832795
+#endif
+
+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]);
+
+
+    // BCB4 does not implement the standard lower_bound function.
+    // This is used instead:
+    const FuncDefinition* fp_lower_bound(const FuncDefinition* first,
+                                         const FuncDefinition* last,
+                                         const FuncDefinition& value)
+    {
+        while(first < last)
+        {
+            const FuncDefinition* middle = first+(last-first)/2;
+            if(*middle < value) first = middle+1;
+            else last = middle;
+        }
+        return last;
+    }
+
+
+    // Returns a pointer to the FuncDefinition instance which 'name' is
+    // the same as the one given by 'F'. If no such function name exists,
+    // returns 0.
+    inline const FuncDefinition* FindFunction(const char* F)
+    {
+        FuncDefinition func = { F, 0, 0, 0 };
+        while(isalnum(F[func.nameLength])) ++func.nameLength;
+        if(func.nameLength)
+        {
+            const FuncDefinition* found =
+                fp_lower_bound(Functions, Functions+FUNC_AMOUNT, func);
+            if(found == Functions+FUNC_AMOUNT || func < *found)
+                return 0;
+            return found;
+        }
+        return 0;
+    }
+};
+
+
+//---------------------------------------------------------------------------
+// Copy-on-write method
+//---------------------------------------------------------------------------
+inline void FunctionParser::copyOnWrite()
+{
+    if(data->referenceCounter > 1)
+    {
+        Data* oldData = data;
+        data = new Data(*oldData);
+        --(oldData->referenceCounter);
+        data->referenceCounter = 1;
+    }
+}
+
+
+//---------------------------------------------------------------------------
+// Constructors and destructors
+//---------------------------------------------------------------------------
+//===========================================================================
+FunctionParser::FunctionParser():
+    parseErrorType(FP_NO_ERROR), evalErrorType(0),
+    data(new Data)
+{
+    data->referenceCounter = 1;
+}
+
+FunctionParser::~FunctionParser()
+{
+    if(--(data->referenceCounter) == 0)
+    {
+        delete data;
+    }
+}
+
+FunctionParser::FunctionParser(const FunctionParser& cpy):
+    parseErrorType(cpy.parseErrorType),
+    evalErrorType(cpy.evalErrorType),
+    data(cpy.data)
+{
+    ++(data->referenceCounter);
+}
+
+FunctionParser& FunctionParser::operator=(const FunctionParser& cpy)
+{
+    if(data != cpy.data)
+    {
+        if(--(data->referenceCounter) == 0) delete data;
+
+        parseErrorType = cpy.parseErrorType;
+        evalErrorType = cpy.evalErrorType;
+        data = cpy.data;
+
+        ++(data->referenceCounter);
+    }
+
+    return *this;
+}
+
+
+FunctionParser::Data::Data():
+    useDegreeConversion(false),
+    ByteCode(0), ByteCodeSize(0),
+    Immed(0), ImmedSize(0),
+    Stack(0), StackSize(0)
+{}
+
+FunctionParser::Data::~Data()
+{
+    if(ByteCode) { delete[] ByteCode; ByteCode=0; }
+    if(Immed) { delete[] Immed; Immed=0; }
+    if(Stack) { delete[] Stack; Stack=0; }
+}
+
+// Makes a deep-copy of Data:
+FunctionParser::Data::Data(const Data& cpy):
+    varAmount(cpy.varAmount), useDegreeConversion(cpy.useDegreeConversion),
+    Variables(cpy.Variables), Constants(cpy.Constants),
+    FuncPtrNames(cpy.FuncPtrNames), FuncPtrs(cpy.FuncPtrs),
+    FuncParserNames(cpy.FuncParserNames), FuncParsers(cpy.FuncParsers),
+    ByteCode(0), ByteCodeSize(cpy.ByteCodeSize),
+    Immed(0), ImmedSize(cpy.ImmedSize),
+    Stack(0), StackSize(cpy.StackSize)
+{
+    if(ByteCodeSize) ByteCode = new unsigned[ByteCodeSize];
+    if(ImmedSize) Immed = new double[ImmedSize];
+    if(StackSize) Stack = new double[StackSize];
+
+    for(unsigned i=0; i<ByteCodeSize; ++i) ByteCode[i] = cpy.ByteCode[i];
+    for(unsigned i=0; i<ImmedSize; ++i) Immed[i] = cpy.Immed[i];
+
+    // No need to copy the stack contents because it's obsolete outside Eval()
+}
+
+
+//---------------------------------------------------------------------------
+// Function parsing
+//---------------------------------------------------------------------------
+//===========================================================================
+namespace
+{
+    // Error messages returned by ErrorMsg():
+    const char* ParseErrorMessage[]=
+    {
+        "Syntax error",                             // 0
+        "Mismatched parenthesis",                   // 1
+        "Missing ')'",                              // 2
+        "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
+        "Syntax error in parameter 'Vars' given to "
+        "FunctionParser::Parse()",                  // 7
+        "Illegal number of parameters to function", // 8
+        "Syntax error: Premature end of string",    // 9
+        "Syntax error: Expecting ( after function", // 10
+        ""
+    };
+
+
+    // Parse variables
+    bool ParseVars(const string& Vars, map<string, unsigned>& dest)
+    {
+        unsigned varNumber = VarBegin;
+        unsigned ind1 = 0, ind2;
+
+        while(ind1 < Vars.size())
+        {
+            if(!isalpha(Vars[ind1]) && Vars[ind1]!='_') return false;
+            for(ind2=ind1+1; ind2<Vars.size() && Vars[ind2]!=','; ++ind2)
+                if(!isalnum(Vars[ind2]) && Vars[ind2]!='_') return false;
+            const string varName = Vars.substr(ind1, ind2-ind1);
+
+            if(dest.insert(make_pair(varName, varNumber++)).second == false)
+                return false;
+
+            ind1 = ind2+1;
+        }
+        return true;
+    }
+};
+
+bool FunctionParser::isValidName(const std::string& name) const
+{
+    if(name.empty() || (!isalpha(name[0]) && name[0] != '_')) return false;
+    for(unsigned i=0; i<name.size(); ++i)
+        if(!isalnum(name[i]) && name[i] != '_') return false;
+
+    if(FindFunction(name.c_str())) return false;
+
+    return true;
+}
+
+
+// Constants:
+bool FunctionParser::AddConstant(const string& name, double value)
+{
+    if(isValidName(name))
+    {
+        const char* n = name.c_str();
+        if(FindVariable(n, data->FuncParserNames) !=
+           data->FuncParserNames.end() ||
+           FindVariable(n, data->FuncPtrNames) !=
+           data->FuncPtrNames.end())
+            return false;
+
+        copyOnWrite();
+
+        data->Constants[name] = value;
+        return true;
+    }
+    return false;
+}
+
+// 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))
+    {
+        const char* n = name.c_str();
+        if(FindVariable(n, data->FuncParserNames) !=
+           data->FuncParserNames.end() ||
+           FindConstant(n) != data->Constants.end())
+            return false;
+
+        copyOnWrite();
+
+        data->FuncPtrNames[name] = data->FuncPtrs.size();
+        data->FuncPtrs.push_back(Data::FuncPtrData(func, paramsAmount));
+        return true;
+    }
+    return false;
+}
+
+bool FunctionParser::checkRecursiveLinking(const FunctionParser* fp) const
+{
+    if(fp == this) return true;
+    for(unsigned i=0; i<fp->data->FuncParsers.size(); ++i)
+        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))
+    {
+        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;
+
+        copyOnWrite();
+
+        data->FuncParserNames[name] = data->FuncParsers.size();
+        data->FuncParsers.push_back(&parser);
+        return true;
+    }
+    return false;
+}
+
+
+
+// Main parsing function
+// ---------------------
+int FunctionParser::Parse(const std::string& Function,
+                          const std::string& Vars,
+                          bool useDegrees)
+{
+    copyOnWrite();
+
+    data->Variables.clear();
+
+    if(!ParseVars(Vars, data->Variables))
+    {
+        parseErrorType = INVALID_VARS;
+        return Function.size();
+    }
+    data->varAmount = data->Variables.size(); // this is for Eval()
+
+    const char* Func = Function.c_str();
+
+    parseErrorType = FP_NO_ERROR;
+
+    int Result = CheckSyntax(Func);
+    if(Result>=0) return Result;
+
+    data->useDegreeConversion = useDegrees;
+    if(!Compile(Func)) return Function.size();
+
+    data->Variables.clear();
+
+    parseErrorType = FP_NO_ERROR;
+    return -1;
+}
+
+namespace
+{
+    // Is given char an operator?
+    inline bool IsOperator(int c)
+    {
+        return strchr("+-*/%^=<>&|",c)!=NULL;
+    }
+
+    // skip whitespace
+    inline void sws(const char* F, int& Ind)
+    {
+        while(F[Ind] && isspace(F[Ind])) ++Ind;
+    }
+};
+
+// Returns an iterator to the variable with the same name as 'F', or to
+// Variables.end() if no such variable exists:
+inline FunctionParser::Data::VarMap_t::const_iterator
+FunctionParser::FindVariable(const char* F, const Data::VarMap_t& vars) const
+{
+    if(vars.size())
+    {
+        unsigned ind = 0;
+        while(isalnum(F[ind]) || F[ind] == '_') ++ind;
+        if(ind)
+        {
+            string name(F, ind);
+            return vars.find(name);
+        }
+    }
+    return vars.end();
+}
+
+inline FunctionParser::Data::ConstMap_t::const_iterator
+FunctionParser::FindConstant(const char* F) const
+{
+    if(data->Constants.size())
+    {
+        unsigned ind = 0;
+        while(isalnum(F[ind]) || F[ind] == '_') ++ind;
+        if(ind)
+        {
+            string name(F, ind);
+            return data->Constants.find(name);
+        }
+    }
+    return data->Constants.end();
+}
+
+//---------------------------------------------------------------------------
+// Check function string syntax
+// ----------------------------
+int FunctionParser::CheckSyntax(const char* Function)
+{
+    const Data::VarMap_t& Variables = data->Variables;
+    const Data::ConstMap_t& Constants = data->Constants;
+    const Data::VarMap_t& FuncPtrNames = data->FuncPtrNames;
+    const Data::VarMap_t& FuncParserNames = data->FuncParserNames;
+
+    vector<int> functionParenthDepth;
+
+    int Ind=0, ParenthCnt=0, c;
+    char* Ptr;
+
+    while(true)
+    {
+        sws(Function, Ind);
+        c=Function[Ind];
+
+// Check for valid operand (must appear)
+
+        // Check for leading -
+        if(c=='-') { sws(Function, ++Ind); c=Function[Ind]; }
+        if(c==0) { parseErrorType=PREMATURE_EOS; return Ind; }
+
+        // Check for math function
+        bool foundFunc = false;
+        const FuncDefinition* fptr = FindFunction(&Function[Ind]);
+        if(fptr)
+        {
+            Ind += fptr->nameLength;
+            foundFunc = true;
+        }
+        else
+        {
+            // Check for user-defined function
+            Data::VarMap_t::const_iterator fIter =
+                FindVariable(&Function[Ind], FuncPtrNames);
+            if(fIter != FuncPtrNames.end())
+            {
+                Ind += fIter->first.size();
+                foundFunc = true;
+            }
+            else
+            {
+                Data::VarMap_t::const_iterator pIter =
+                    FindVariable(&Function[Ind], FuncParserNames);
+                if(pIter != FuncParserNames.end())
+                {
+                    Ind += pIter->first.size();
+                    foundFunc = true;
+                }
+            }
+        }
+
+        if(foundFunc)
+        {
+            sws(Function, Ind);
+            c = Function[Ind];
+            if(c!='(') { parseErrorType=EXPECT_PARENTH_FUNC; return Ind; }
+            functionParenthDepth.push_back(ParenthCnt+1);
+        }
+
+        // Check for opening parenthesis
+        if(c=='(')
+        {
+            ++ParenthCnt;
+            sws(Function, ++Ind);
+            if(Function[Ind]==')') { parseErrorType=EMPTY_PARENTH; return Ind;}
+            continue;
+        }
+
+        // Check for number
+        if(isdigit(c) || (c=='.' && isdigit(Function[Ind+1])))
+        {
+            strtod(&Function[Ind], &Ptr);
+            Ind += int(Ptr-&Function[Ind]);
+            sws(Function, Ind);
+            c = Function[Ind];
+        }
+        else
+        { // Check for variable
+            Data::VarMap_t::const_iterator vIter =
+                FindVariable(&Function[Ind], Variables);
+            if(vIter != Variables.end())
+                Ind += vIter->first.size();
+            else
+            {
+                // Check for constant
+                Data::ConstMap_t::const_iterator cIter =
+                    FindConstant(&Function[Ind]);
+                if(cIter != Constants.end())
+                    Ind += cIter->first.size();
+                else
+                { parseErrorType=SYNTAX_ERROR; return Ind; }
+            }
+            sws(Function, Ind);
+            c = Function[Ind];
+        }
+
+        // Check for closing parenthesis
+        while(c==')')
+        {
+            if(functionParenthDepth.size() &&
+               functionParenthDepth.back() == ParenthCnt)
+                functionParenthDepth.pop_back();
+            if((--ParenthCnt)<0) { parseErrorType=MISM_PARENTH; return Ind; }
+            sws(Function, ++Ind);
+            c=Function[Ind];
+        }
+
+// If we get here, we have a legal operand and now a legal operator or
+// end of string must follow
+
+        // 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))
+        { 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;
+    } // while
+
+    // Check that all opened parentheses are also closed
+    if(ParenthCnt>0) { parseErrorType=MISSING_PARENTH; return Ind; }
+
+// The string is ok
+    parseErrorType=FP_NO_ERROR;
+    return -1;
+}
+
+
+// Compile function string to bytecode
+// -----------------------------------
+bool 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; }
+
+    vector<unsigned> byteCode; byteCode.reserve(1024);
+    tempByteCode = &byteCode;
+
+    vector<double> immed; immed.reserve(1024);
+    tempImmed = &immed;
+
+    data->StackSize = StackPtr = 0;
+
+    CompileExpression(Function, 0);
+    if(parseErrorType != FP_NO_ERROR) return false;
+
+    data->ByteCodeSize = byteCode.size();
+    data->ImmedSize = immed.size();
+
+    if(data->ByteCodeSize)
+    {
+        data->ByteCode = new unsigned[data->ByteCodeSize];
+        memcpy(data->ByteCode, &byteCode[0],
+               sizeof(unsigned)*data->ByteCodeSize);
+    }
+    if(data->ImmedSize)
+    {
+        data->Immed = new double[data->ImmedSize];
+        memcpy(data->Immed, &immed[0],
+               sizeof(double)*data->ImmedSize);
+    }
+    if(data->StackSize)
+        data->Stack = new double[data->StackSize];
+
+    return true;
+}
+
+
+inline void FunctionParser::AddCompiledByte(unsigned c)
+{
+    tempByteCode->push_back(c);
+}
+
+inline void FunctionParser::AddImmediate(double i)
+{
+    tempImmed->push_back(i);
+}
+
+inline void FunctionParser::AddFunctionOpcode(unsigned opcode)
+{
+    if(data->useDegreeConversion)
+        switch(opcode)
+        {
+          case cCos:
+          case cCosh:
+          case cCot:
+          case cCsc:
+          case cSec:
+          case cSin:
+          case cSinh:
+          case cTan:
+          case cTanh:
+              AddCompiledByte(cRad);
+        }
+
+    AddCompiledByte(opcode);
+
+    if(data->useDegreeConversion)
+        switch(opcode)
+        {
+          case cAcos:
+#ifndef NO_ASINH
+          case cAcosh:
+          case cAsinh:
+          case cAtanh:
+#endif
+          case cAsin:
+          case cAtan:
+          case cAtan2:
+              AddCompiledByte(cDeg);
+        }
+}
+
+inline void FunctionParser::incStackPtr()
+{
+    if(++StackPtr > data->StackSize) ++(data->StackSize);
+}
+
+
+// Compile if()
+int FunctionParser::CompileIf(const char* F, int ind)
+{
+    int ind2 = CompileExpression(F, ind, true); // condition
+    sws(F, ind2);
+    if(F[ind2] != ',') { parseErrorType=ILL_PARAMS_AMOUNT; return ind2; }
+    AddCompiledByte(cIf);
+    unsigned curByteCodeSize = tempByteCode->size();
+    AddCompiledByte(0); // Jump index; to be set later
+    AddCompiledByte(0); // Immed jump index; to be set later
+
+    --StackPtr;
+
+    ind2 = CompileExpression(F, ind2+1, true); // then
+    sws(F, ind2);
+    if(F[ind2] != ',') { parseErrorType=ILL_PARAMS_AMOUNT; return ind2; }
+    AddCompiledByte(cJump);
+    unsigned curByteCodeSize2 = tempByteCode->size();
+    unsigned curImmedSize2 = tempImmed->size();
+    AddCompiledByte(0); // Jump index; to be set later
+    AddCompiledByte(0); // Immed jump index; to be set later
+
+    --StackPtr;
+
+    ind2 = CompileExpression(F, ind2+1, true); // else
+    sws(F, ind2);
+    if(F[ind2] != ')') { parseErrorType=ILL_PARAMS_AMOUNT; return ind2; }
+
+    // Set jump indices
+    (*tempByteCode)[curByteCodeSize] = curByteCodeSize2+1;
+    (*tempByteCode)[curByteCodeSize+1] = curImmedSize2;
+    (*tempByteCode)[curByteCodeSize2] = tempByteCode->size()-1;
+    (*tempByteCode)[curByteCodeSize2+1] = tempImmed->size();
+
+    return ind2+1;
+}
+
+int FunctionParser::CompileFunctionParams(const char* F, int ind,
+                                          unsigned requiredParams)
+{
+    unsigned curStackPtr = StackPtr;
+    int ind2 = CompileExpression(F, ind);
+
+    if(StackPtr != curStackPtr+requiredParams)
+    { parseErrorType=ILL_PARAMS_AMOUNT; return ind; }
+
+    StackPtr -= requiredParams - 1;
+    sws(F, ind2);
+    return ind2+1; // F[ind2] is ')'
+}
+
+// Compiles element
+int FunctionParser::CompileElement(const char* F, int ind)
+{
+    sws(F, ind);
+    char c = F[ind];
+
+    if(c == '(')
+    {
+        ind = CompileExpression(F, ind+1);
+        sws(F, ind);
+        return ind+1; // F[ind] is ')'
+    }
+
+    if(isdigit(c) || c=='.' /*|| c=='-'*/) // Number
+    {
+        const char* startPtr = &F[ind];
+        char* endPtr;
+        double val = strtod(startPtr, &endPtr);
+        AddImmediate(val);
+        AddCompiledByte(cImmed);
+        incStackPtr();
+        return ind+(endPtr-startPtr);
+    }
+
+    if(isalpha(c) || c == '_') // Function, variable or constant
+    {
+        const FuncDefinition* func = FindFunction(F+ind);
+        if(func) // is function
+        {
+            int ind2 = ind + func->nameLength;
+            sws(F, ind2); // F[ind2] is '('
+            if(strcmp(func->name, "if") == 0) // "if" is a special case
+            {
+                return CompileIf(F, ind2+1);
+            }
+
+#ifndef DISABLE_EVAL
+            unsigned requiredParams =
+                strcmp(func->name, "eval") == 0 ?
+                data->Variables.size() : func->params;
+#else
+            unsigned requiredParams = func->params;
+#endif
+            ind2 = CompileFunctionParams(F, ind2+1, requiredParams);
+            AddFunctionOpcode(func->opcode);
+            return ind2; // F[ind2-1] is ')'
+        }
+
+        Data::VarMap_t::const_iterator vIter =
+            FindVariable(F+ind, data->Variables);
+        if(vIter != data->Variables.end()) // is variable
+        {
+            AddCompiledByte(vIter->second);
+            incStackPtr();
+            return ind + vIter->first.size();
+        }
+
+        Data::ConstMap_t::const_iterator cIter = FindConstant(F+ind);
+        if(cIter != data->Constants.end()) // is constant
+        {
+            AddImmediate(cIter->second);
+            AddCompiledByte(cImmed);
+            incStackPtr();
+            return ind + cIter->first.size();
+        }
+
+        Data::VarMap_t::const_iterator fIter =
+            FindVariable(F+ind, data->FuncPtrNames);
+        if(fIter != data->FuncPtrNames.end()) // is user-defined func pointer
+        {
+            unsigned index = fIter->second;
+
+            int ind2 = ind + fIter->first.length();
+            sws(F, ind2); // F[ind2] is '('
+
+            ind2 = CompileFunctionParams(F, ind2+1,
+                                         data->FuncPtrs[index].params);
+
+            AddCompiledByte(cFCall);
+            AddCompiledByte(index);
+            return ind2;
+        }
+
+        Data::VarMap_t::const_iterator pIter =
+            FindVariable(F+ind, data->FuncParserNames);
+        if(pIter != data->FuncParserNames.end()) // is user-defined func parser
+        {
+            unsigned index = pIter->second;
+
+            int ind2 = ind + pIter->first.length();
+            sws(F, ind2); // F[ind2] is '('
+
+            ind2 = CompileFunctionParams
+                (F, ind2+1, data->FuncParsers[index]->data->varAmount);
+
+            AddCompiledByte(cPCall);
+            AddCompiledByte(index);
+            return ind2;
+        }
+    }
+
+    parseErrorType = UNEXPECTED_ERROR;
+    return ind;
+}
+
+// Compiles '^'
+int FunctionParser::CompilePow(const char* F, int ind)
+{
+    int ind2 = CompileElement(F, ind);
+    sws(F, ind2);
+
+    while(F[ind2] == '^')
+    {
+        ind2 = CompileUnaryMinus(F, ind2+1);
+        sws(F, ind2);
+        AddCompiledByte(cPow);
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles unary '-'
+int FunctionParser::CompileUnaryMinus(const char* F, int ind)
+{
+    sws(F, ind);
+    if(F[ind] == '-')
+    {
+        int ind2 = ind+1;
+        sws(F, ind2);
+        ind2 = CompilePow(F, ind2);
+        sws(F, ind2);
+
+        // if we are negating a constant, negate the constant itself:
+        if(tempByteCode->back() == cImmed)
+            tempImmed->back() = -tempImmed->back();
+
+        // if we are negating a negation, we can remove both:
+        else if(tempByteCode->back() == cNeg)
+            tempByteCode->pop_back();
+
+        else
+            AddCompiledByte(cNeg);
+
+        return ind2;
+    }
+
+    int ind2 = CompilePow(F, ind);
+    sws(F, ind2);
+    return ind2;
+}
+
+// Compiles '*', '/' and '%'
+int FunctionParser::CompileMult(const char* F, int ind)
+{
+    int ind2 = CompileUnaryMinus(F, ind);
+    sws(F, ind2);
+    char op;
+
+    while((op = F[ind2]) == '*' || op == '/' || op == '%')
+    {
+        ind2 = CompileUnaryMinus(F, ind2+1);
+        sws(F, ind2);
+        switch(op)
+        {
+          case '*': AddCompiledByte(cMul); break;
+          case '/': AddCompiledByte(cDiv); break;
+          case '%': AddCompiledByte(cMod); break;
+        }
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles '+' and '-'
+int FunctionParser::CompileAddition(const char* F, int ind)
+{
+    int ind2 = CompileMult(F, ind);
+    sws(F, ind2);
+    char op;
+
+    while((op = F[ind2]) == '+' || op == '-')
+    {
+        ind2 = CompileMult(F, ind2+1);
+        sws(F, ind2);
+        AddCompiledByte(op=='+' ? cAdd : cSub);
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles '=', '<' and '>'
+int FunctionParser::CompileComparison(const char* F, int ind)
+{
+    int ind2 = CompileAddition(F, ind);
+    sws(F, ind2);
+    char op;
+
+    while((op = F[ind2]) == '=' || op == '<' || op == '>')
+    {
+        ind2 = CompileAddition(F, ind2+1);
+        sws(F, ind2);
+        switch(op)
+        {
+          case '=': AddCompiledByte(cEqual); break;
+          case '<': AddCompiledByte(cLess); break;
+          case '>': AddCompiledByte(cGreater); break;
+        }
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles '&'
+int FunctionParser::CompileAnd(const char* F, int ind)
+{
+    int ind2 = CompileComparison(F, ind);
+    sws(F, ind2);
+
+    while(F[ind2] == '&')
+    {
+        ind2 = CompileComparison(F, ind2+1);
+        sws(F, ind2);
+        AddCompiledByte(cAnd);
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles '|'
+int FunctionParser::CompileOr(const char* F, int ind)
+{
+    int ind2 = CompileAnd(F, ind);
+    sws(F, ind2);
+
+    while(F[ind2] == '|')
+    {
+        ind2 = CompileAnd(F, ind2+1);
+        sws(F, ind2);
+        AddCompiledByte(cOr);
+        --StackPtr;
+    }
+
+    return ind2;
+}
+
+// Compiles ','
+int FunctionParser::CompileExpression(const char* F, int ind, bool stopAtComma)
+{
+    int ind2 = CompileOr(F, ind);
+    sws(F, ind2);
+
+    if(stopAtComma) return ind2;
+
+    while(F[ind2] == ',')
+    {
+        ind2 = CompileOr(F, ind2+1);
+        sws(F, ind2);
+    }
+
+    return ind2;
+}
+
+
+// Return parse error message
+// --------------------------
+const char* FunctionParser::ErrorMsg() const
+{
+    if(parseErrorType != FP_NO_ERROR) return ParseErrorMessage[parseErrorType];
+    return 0;
+}
+
+//---------------------------------------------------------------------------
+// Function evaluation
+//---------------------------------------------------------------------------
+//===========================================================================
+namespace
+{
+    inline int doubleToInt(double d)
+    {
+        return d<0 ? -int((-d)+.5) : int(d+.5);
+    }
+
+    inline double Min(double d1, double d2)
+    {
+        return d1<d2 ? d1 : d2;
+    }
+    inline double Max(double d1, double d2)
+    {
+        return d1>d2 ? d1 : d2;
+    }
+
+
+    inline double DegreesToRadians(double degrees)
+    {
+        return degrees*(M_PI/180.0);
+    }
+    inline double RadiansToDegrees(double radians)
+    {
+        return radians*(180.0/M_PI);
+    }
+}
+
+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;
+
+    for(IP=0; IP<ByteCodeSize; ++IP)
+    {
+        switch(ByteCode[IP])
+        {
+// Functions:
+          case   cAbs: Stack[SP] = fabs(Stack[SP]); break;
+          case  cAcos: if(Stack[SP] < -1 || Stack[SP] > 1)
+                       { evalErrorType=4; return 0; }
+                       Stack[SP] = acos(Stack[SP]); break;
+#ifndef NO_ASINH
+          case cAcosh: Stack[SP] = acosh(Stack[SP]); break;
+#endif
+          case  cAsin: if(Stack[SP] < -1 || Stack[SP] > 1)
+                       { evalErrorType=4; return 0; }
+                       Stack[SP] = asin(Stack[SP]); break;
+#ifndef NO_ASINH
+          case cAsinh: Stack[SP] = asinh(Stack[SP]); break;
+#endif
+          case  cAtan: Stack[SP] = atan(Stack[SP]); break;
+          case cAtan2: Stack[SP-1] = atan2(Stack[SP-1], Stack[SP]);
+                       --SP; break;
+#ifndef NO_ASINH
+          case cAtanh: Stack[SP] = atanh(Stack[SP]); break;
+#endif
+          case  cCeil: Stack[SP] = ceil(Stack[SP]); break;
+          case   cCos: Stack[SP] = cos(Stack[SP]); break;
+          case  cCosh: Stack[SP] = cosh(Stack[SP]); break;
+
+          case   cCot:
+              {
+                  double t = tan(Stack[SP]);
+                  if(t == 0) { evalErrorType=1; return 0; }
+                  Stack[SP] = 1/t; break;
+              }
+          case   cCsc:
+              {
+                  double s = sin(Stack[SP]);
+                  if(s == 0) { evalErrorType=1; return 0; }
+                  Stack[SP] = 1/s; break;
+              }
+
+
+#ifndef DISABLE_EVAL
+          case  cEval:
+              {
+                  data->Stack = new double[data->StackSize];
+                  double retVal = Eval(&Stack[SP-data->varAmount+1]);
+                  delete[] data->Stack;
+                  data->Stack = Stack;
+                  SP -= data->varAmount-1;
+                  Stack[SP] = retVal;
+                  break;
+              }
+#endif
+
+          case   cExp: Stack[SP] = exp(Stack[SP]); break;
+          case cFloor: Stack[SP] = floor(Stack[SP]); break;
+
+          case    cIf:
+              {
+                  unsigned jumpAddr = ByteCode[++IP];
+                  unsigned immedAddr = ByteCode[++IP];
+                  if(doubleToInt(Stack[SP]) == 0)
+                  {
+                      IP = jumpAddr;
+                      DP = immedAddr;
+                  }
+                  --SP; break;
+              }
+
+          case   cInt: Stack[SP] = floor(Stack[SP]+.5); break;
+          case   cLog: if(Stack[SP] <= 0) { evalErrorType=3; return 0; }
+                       Stack[SP] = log(Stack[SP]); break;
+          case cLog10: if(Stack[SP] <= 0) { evalErrorType=3; return 0; }
+                       Stack[SP] = log10(Stack[SP]); break;
+          case   cMax: Stack[SP-1] = Max(Stack[SP-1], Stack[SP]);
+                       --SP; break;
+          case   cMin: Stack[SP-1] = Min(Stack[SP-1], Stack[SP]);
+                       --SP; break;
+          case   cSec:
+              {
+                  double c = cos(Stack[SP]);
+                  if(c == 0) { evalErrorType=1; return 0; }
+                  Stack[SP] = 1/c; break;
+              }
+          case   cSin: Stack[SP] = sin(Stack[SP]); break;
+          case  cSinh: Stack[SP] = sinh(Stack[SP]); break;
+          case  cSqrt: if(Stack[SP] < 0) { evalErrorType=2; return 0; }
+                       Stack[SP] = sqrt(Stack[SP]); break;
+          case   cTan: Stack[SP] = tan(Stack[SP]); break;
+          case  cTanh: Stack[SP] = tanh(Stack[SP]); break;
+
+
+// Misc:
+          case cImmed: Stack[++SP] = Immed[DP++]; break;
+          case  cJump: DP = ByteCode[IP+2];
+                       IP = ByteCode[IP+1];
+                       break;
+
+// Operators:
+          case   cNeg: Stack[SP] = -Stack[SP]; break;
+          case   cAdd: Stack[SP-1] += Stack[SP]; --SP; break;
+          case   cSub: Stack[SP-1] -= Stack[SP]; --SP; break;
+          case   cMul: Stack[SP-1] *= Stack[SP]; --SP; break;
+          case   cDiv: if(Stack[SP] == 0) { evalErrorType=1; return 0; }
+                       Stack[SP-1] /= Stack[SP]; --SP; break;
+          case   cMod: if(Stack[SP] == 0) { evalErrorType=1; return 0; }
+                       Stack[SP-1] = fmod(Stack[SP-1], Stack[SP]);
+                       --SP; break;
+          case   cPow: Stack[SP-1] = pow(Stack[SP-1], Stack[SP]);
+                       --SP; break;
+
+          case cEqual: Stack[SP-1] = (Stack[SP-1] == Stack[SP]);
+                       --SP; break;
+          case  cLess: Stack[SP-1] = (Stack[SP-1] < Stack[SP]);
+                       --SP; break;
+          case cGreater: Stack[SP-1] = (Stack[SP-1] > Stack[SP]);
+                         --SP; break;
+          case   cAnd: Stack[SP-1] =
+                           (doubleToInt(Stack[SP-1]) &&
+                            doubleToInt(Stack[SP]));
+                       --SP; break;
+          case    cOr: Stack[SP-1] =
+                           (doubleToInt(Stack[SP-1]) ||
+                            doubleToInt(Stack[SP]));
+                       --SP; break;
+
+// Degrees-radians conversion:
+          case   cDeg: Stack[SP] = RadiansToDegrees(Stack[SP]); break;
+          case   cRad: Stack[SP] = DegreesToRadians(Stack[SP]); break;
+
+// User-defined function calls:
+          case cFCall:
+              {
+                  unsigned index = ByteCode[++IP];
+                  unsigned params = data->FuncPtrs[index].params;
+                  double retVal =
+                      data->FuncPtrs[index].ptr(&Stack[SP-params+1]);
+                  SP -= params-1;
+                  Stack[SP] = retVal;
+                  break;
+              }
+
+          case cPCall:
+              {
+                  unsigned index = ByteCode[++IP];
+                  unsigned params = data->FuncParsers[index]->data->varAmount;
+                  double retVal =
+                      data->FuncParsers[index]->Eval(&Stack[SP-params+1]);
+                  SP -= params-1;
+                  Stack[SP] = retVal;
+                  break;
+              }
+
+
+#ifdef SUPPORT_OPTIMIZER
+          case   cVar: break; // Paranoia. These should never exist
+          case   cDup: Stack[SP+1] = Stack[SP]; ++SP; break;
+          case   cInv:
+              if(Stack[SP] == 0.0) { evalErrorType=1; return 0; }
+              Stack[SP] = 1.0/Stack[SP];
+              break;
+#endif
+
+// Variables:
+          default:
+              Stack[++SP] = Vars[ByteCode[IP]-VarBegin];
+        }
+    }
+
+    evalErrorType=0;
+    return Stack[SP];
+}
+
+
+#ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT
+namespace
+{
+    inline void printHex(std::ostream& dest, unsigned n)
+    {
+        dest.width(8); dest.fill('0'); hex(dest); //uppercase(dest);
+        dest << n;
+    }
+}
+
+void FunctionParser::PrintByteCode(std::ostream& dest) const
+{
+    const unsigned* const ByteCode = data->ByteCode;
+    const double* const Immed = data->Immed;
+
+    for(unsigned IP=0, DP=0; IP<data->ByteCodeSize; ++IP)
+    {
+        printHex(dest, IP);
+        dest << ": ";
+
+        unsigned opcode = ByteCode[IP];
+
+        switch(opcode)
+        {
+          case cIf:
+              dest << "jz\t";
+              printHex(dest, ByteCode[IP+1]+1);
+              dest << endl;
+              IP += 2;
+              break;
+
+          case cJump:
+              dest << "jump\t";
+              printHex(dest, ByteCode[IP+1]+1);
+              dest << endl;
+              IP += 2;
+              break;
+          case cImmed:
+              dest.precision(10);
+              dest << "push\t" << Immed[DP++] << endl;
+              break;
+
+          case cFCall:
+              {
+                  unsigned index = ByteCode[++IP];
+                  Data::VarMap_t::const_iterator iter =
+                      data->FuncPtrNames.begin();
+                  while(iter->second != index) ++iter;
+                  dest << "call\t" << iter->first << endl;
+                  break;
+              }
+
+          case cPCall:
+              {
+                  unsigned index = ByteCode[++IP];
+                  Data::VarMap_t::const_iterator iter =
+                      data->FuncParserNames.begin();
+                  while(iter->second != index) ++iter;
+                  dest << "call\t" << iter->first << endl;
+                  break;
+              }
+
+          default:
+              if(opcode < VarBegin)
+              {
+                  string n;
+                  switch(opcode)
+                  {
+                    case cNeg: n = "neg"; break;
+                    case cAdd: n = "add"; break;
+                    case cSub: n = "sub"; break;
+                    case cMul: n = "mul"; break;
+                    case cDiv: n = "div"; break;
+                    case cMod: n = "mod"; break;
+                    case cPow: n = "pow"; break;
+                    case cEqual: n = "eq"; break;
+                    case cLess: n = "lt"; break;
+                    case cGreater: n = "gt"; break;
+                    case cAnd: n = "and"; break;
+                    case cOr: n = "or"; break;
+                    case cDeg: n = "deg"; break;
+                    case cRad: n = "rad"; break;
+
+#ifndef DISABLE_EVAL
+                    case cEval: n = "call\t0"; break;
+#endif
+
+#ifdef SUPPORT_OPTIMIZER
+                    case cVar: n = "(var)"; break;
+                    case cDup: n = "dup"; break;
+                    case cInv: n = "inv"; break;
+#endif
+
+                    default: n = Functions[opcode-cAbs].name;
+                  }
+                  dest << n << endl;
+              }
+              else
+              {
+                  dest << "push\tVar" << opcode-VarBegin << endl;
+              }
+        }
+    }
+}
+#endif
+
+
+//========================================================================
+// Optimization code was contributed by Bisqwit (http://iki.fi/bisqwit/)
+//========================================================================
+#ifdef SUPPORT_OPTIMIZER
+
+#include <list>
+#include <utility>
+
+#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; }
+    // 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<SubTree> 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 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<CodeTreeData, unsigned> 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<unsigned v> 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<pit> cp;
+        double value;
+        unsigned size() const { return cp.size(); }
+    };
+    struct ConstList BuildConstList();
+    void KillConst(const ConstList &cl)
+    {
+        for(list<pit>::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<pit> 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<pit>::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<SubTree> 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<SubTree>::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<pit> 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<pit>::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<pit> 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<pit>::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<unsigned> &byteCode,
+                  vector<double>   &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<pit> &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<pit>::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<unsigned> &byteCode,
+    vector<double>   &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<CodeTree*>(r);
+
+    vector<CodeTree> 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<b; ++a) \
+            stack[stacktop].AddParam(stack[newstacktop+a]); \
+        stack.erase(stack.begin() + newstacktop, \
+                    stack.begin() + stacktop); \
+        stacktop = newstacktop; GROW(1); \
+    } while(0)
+
+    #define ADDCONST(n) do { \
+        stack[stacktop].SetImmed((n)); \
+        GROW(1); \
+    } while(0)
+
+    unsigned stacktop=0;
+
+    list<unsigned> 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<unsigned> byteCode;
+    vector<double> immed;
+
+#if 0
+    byteCode.resize(Comp.ByteCodeSize);
+    for(unsigned a=0; a<Comp.ByteCodeSize; ++a)byteCode[a] = Comp.ByteCode[a];
+
+    immed.resize(Comp.ImmedSize);
+    for(unsigned a=0; a<Comp.ImmedSize; ++a)immed[a] = Comp.Immed[a];
+#else
+    byteCode.clear(); immed.clear();
+    tree.Assemble(byteCode, immed);
+#endif
+
+    delete[] data->ByteCode; data->ByteCode = 0;
+    if((data->ByteCodeSize = byteCode.size()) > 0)
+    {
+        data->ByteCode = new unsigned[data->ByteCodeSize];
+        for(unsigned a=0; a<byteCode.size(); ++a)
+            data->ByteCode[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; a<immed.size(); ++a)
+            data->Immed[a] = immed[a];
+    }
+}
+
+
+#else /* !SUPPORT_OPTIMIZER */
+
+/* keep the linker happy */
+void FunctionParser::MakeTree(CodeTree *) const {}
+void FunctionParser::Optimize()
+{
+    // Do nothing if no optimizations are supported.
+}
+#endif
diff --git a/deal.II/contrib/functionparser/fparser.h b/deal.II/contrib/functionparser/fparser.h
new file mode 100644 (file)
index 0000000..aee49a2
--- /dev/null
@@ -0,0 +1,155 @@
+/***************************************************************************\
+|* Function parser v2.7 by Warp                                            *|
+|* ----------------------------                                            *|
+|* Parses and evaluates the given function with the given variable values. *|
+|*                                                                         *|
+\***************************************************************************/
+
+#ifndef ONCE_FPARSER_H_
+#define ONCE_FPARSER_H_
+
+#include <string>
+#include <map>
+#include <vector>
+
+#ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT
+#include <iostream>
+#endif
+
+class FunctionParser
+{
+public:
+    enum ParseErrorType
+    {
+        SYNTAX_ERROR=0, MISM_PARENTH, MISSING_PARENTH, EMPTY_PARENTH,
+        EXPECT_OPERATOR, OUT_OF_MEMORY, UNEXPECTED_ERROR, INVALID_VARS,
+        ILL_PARAMS_AMOUNT, PREMATURE_EOS, EXPECT_PARENTH_FUNC,
+        FP_NO_ERROR
+    };
+
+
+    int Parse(const std::string& Function, const std::string& Vars,
+              bool useDegrees = false);
+    const char* ErrorMsg() const;
+    inline ParseErrorType GetParseErrorType() const { return parseErrorType; }
+
+    double Eval(const double* Vars);
+    inline int EvalError() const { return evalErrorType; }
+
+    bool AddConstant(const std::string& name, double value);
+
+    typedef double (*FunctionPtr)(const double*);
+
+    bool AddFunction(const std::string& name,
+                     FunctionPtr, unsigned paramsAmount);
+    bool AddFunction(const std::string& name, FunctionParser&);
+
+    void Optimize();
+
+
+    FunctionParser();
+    ~FunctionParser();
+
+    // Copy constructor and assignment operator (implemented using the
+    // copy-on-write technique for efficiency):
+    FunctionParser(const FunctionParser&);
+    FunctionParser& operator=(const FunctionParser&);
+
+
+#ifdef FUNCTIONPARSER_SUPPORT_DEBUG_OUTPUT
+    // For debugging purposes only:
+    void PrintByteCode(std::ostream& dest) const;
+#endif
+
+
+
+//========================================================================
+private:
+//========================================================================
+
+// Private data:
+// ------------
+    ParseErrorType parseErrorType;
+    int evalErrorType;
+
+    struct Data
+    {
+        unsigned referenceCounter;
+
+        int varAmount;
+        bool useDegreeConversion;
+
+        typedef std::map<std::string, unsigned> VarMap_t;
+        VarMap_t Variables;
+
+        typedef std::map<std::string, double> ConstMap_t;
+        ConstMap_t Constants;
+
+        VarMap_t FuncPtrNames;
+        struct FuncPtrData
+        {
+            FunctionPtr ptr; unsigned params;
+            FuncPtrData(FunctionPtr p, unsigned par): ptr(p), params(par) {}
+        };
+        std::vector<FuncPtrData> FuncPtrs;
+
+        VarMap_t FuncParserNames;
+        std::vector<FunctionParser*> FuncParsers;
+
+        unsigned* ByteCode;
+        unsigned ByteCodeSize;
+        double* Immed;
+        unsigned ImmedSize;
+        double* Stack;
+        unsigned StackSize;
+
+        Data();
+        ~Data();
+        Data(const Data&);
+
+        Data& operator=(const Data&); // not implemented on purpose
+    };
+
+    Data* data;
+
+    // Temp data needed in Compile():
+    unsigned StackPtr;
+    std::vector<unsigned>* tempByteCode;
+    std::vector<double>* tempImmed;
+
+
+// Private methods:
+// ---------------
+    inline void copyOnWrite();
+
+
+    bool checkRecursiveLinking(const FunctionParser*) const;
+
+    bool isValidName(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;
+    int CheckSyntax(const char*);
+    bool Compile(const char*);
+    bool IsVariable(int);
+    void AddCompiledByte(unsigned);
+    void AddImmediate(double);
+    void AddFunctionOpcode(unsigned);
+    inline void incStackPtr();
+    int CompileIf(const char*, int);
+    int CompileFunctionParams(const char*, int, unsigned);
+    int CompileElement(const char*, int);
+    int CompilePow(const char*, int);
+    int CompileUnaryMinus(const char*, int);
+    int CompileMult(const char*, int);
+    int CompileAddition(const char*, int);
+    int CompileComparison(const char*, int);
+    int CompileAnd(const char*, int);
+    int CompileOr(const char*, int);
+    int CompileExpression(const char*, int, bool=false);
+
+
+    void MakeTree(void*) const;
+};
+
+#endif
diff --git a/deal.II/contrib/functionparser/fparser.txt b/deal.II/contrib/functionparser/fparser.txt
new file mode 100644 (file)
index 0000000..cd2263a
--- /dev/null
@@ -0,0 +1,777 @@
+  Function parser for C++  v2.7 by Warp.
+  =====================================
+
+  Optimization code contributed by Bisqwit (http://iki.fi/bisqwit/)
+
+
+  The usage license of this library is located at the end of this text file.
+
+
+
+  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.
+
+
+
+=============================================================================
+  - Preface
+=============================================================================
+
+  Often people need to ask some mathematical expression from the user and
+then evaluate values for that expression. The simplest example is a program
+which draws the graphic of a user-defined function on screen.
+
+  This library adds C-style function string parsing to the program. This
+means that you can evaluate the string "sqrt(1-x^2+y^2)" with given values
+of 'x' and 'y'.
+
+  The library is intended to be very fast. It byte-compiles the function
+string at parse time and interpretes this byte-code at evaluation time.
+The evaluation is straightforward and no recursions are done (uses stack
+arithmetic).
+  Empirical tests show that it indeed is very fast (specially compared to
+libraries which evaluate functions by just interpreting the raw function
+string).
+
+  The library is made in ISO C++ and requires a standard-conforming C++
+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).
+
+
+  * Conditional compiling:
+    ---------------------
+
+    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:
+
+  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.
+
+  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.
+
+  SUPPORT_OPTIMIZER : (Default on)
+       If you are not going to use the Optimize() method, you can comment
+       this line out to speed-up the compilation of fparser.cc a bit, as
+       well as making the binary a bit smaller. (Optimize() can still be
+       called, but it will not do anything.)
+
+
+  * Copying and assignment:
+    ----------------------
+
+    The class implements a safe copy constructor and assignment operator.
+
+    It uses the copy-on-write technique for efficiency. This means that
+  when copying or assigning a FunctionParser instance, the internal data
+  (which in some cases can be quite lengthy) is not immediately copied
+  but only when the contents of the copy (or the original) are changed.
+    This means that copying/assigning is a very fast operation, and if
+  the copies are never modified then actual data copying never happens
+  either.
+
+    The Eval() and EvalError() methods of the copy can be called without
+  the internal data being copied.
+    Calling Parse(), Optimize() or the user-defined constant/function adding
+  methods will cause a deep-copy.
+
+    (C++ basics: The copy constructor is called when a new FunctionParser
+     instance is initialized with another, ie. like:
+
+       FunctionParser fp2 = fp1; // or: FunctionParser fp2(fp1);
+
+     or when a function takes a FunctionParser instance as parameter, eg:
+
+       void foo(FunctionParser p) // takes an instance of FunctionParser
+       { ... }
+
+     The assignment operator is called when a FunctionParser instance is
+     assigned to another, like "fp2 = fp1;".)
+
+
+  * Short descriptions of FunctionParser methods:
+    --------------------------------------------
+
+int Parse(const std::string& Function, const std::string& Vars,
+          bool useDegrees = false);
+
+    Parses the given function and compiles it to internal format.
+    Return value is -1 if successful, else the index value to the location
+    of the error.
+
+
+const char* ErrorMsg(void) const;
+
+    Returns an error message corresponding to the error in Parse(), or 0 if
+    no such error occurred.
+
+
+ParseErrorType GetParseErrorType() const;
+
+    Returns the type of parsing error which occurred. Possible return types
+    are described in the long description.
+
+
+double Eval(const double* Vars);
+
+    Evaluates the function given to Parse().
+
+
+int EvalError(void) const;
+
+    Returns 0 if no error happened in the previous call to Eval(), else an
+    error code >0.
+
+
+void Optimize();
+
+    Tries to optimize the bytecode for faster evaluation.
+
+
+bool AddConstant(const std::string& name, double value);
+
+    Add a constant to the parser. Returns false if the name of the constant
+    is invalid, else true.
+
+
+bool AddFunction(const std::string& name,
+                 double (*functionPtr)(const double*),
+                 unsigned paramsAmount);
+
+    Add a user-defined function to the parser (as a function pointer).
+    Returns false if the name of the function is invalid, else true.
+
+
+bool AddFunction(const std::string& name, FunctionParser&);
+
+    Add a user-defined function to the parser (as a FunctionParser instance).
+    Returns false if the name of the function is invalid, else true.
+
+
+
+  * Long descriptions of FunctionParser methods:
+    -------------------------------------------
+
+---------------------------------------------------------------------------
+int Parse(const std::string& Function, const std::string& Vars,
+          bool useDegrees = false);
+---------------------------------------------------------------------------
+
+      Parses the given function (and compiles it to internal format).
+    Destroys previous function. Following calls to Eval() will evaluate
+    the given function.
+      The strings given as parameters are not needed anymore after parsing.
+
+    Parameters:
+      Function  : String containing the function to parse.
+      Vars      : String containing the variable names, separated by commas.
+                  Eg. "x,y", "VarX,VarY,VarZ,n" or "x1,x2,x3,x4,__VAR__".
+      useDegrees: (Optional.) Whether to use degrees or radians in
+                  trigonometric functions. (Default: radians)
+
+    Variables can have any size and they are case sensitive (ie. "var",
+    "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.
+
+    Using longer variable names causes no overhead whatsoever to the Eval()
+    method, so it's completely safe to use variable names of any size.
+
+    The third, optional parameter specifies whether angles should be
+    interpreted as radians or degrees in trigonometrical functions.
+    If not specified, the default value is radians.
+
+    Return values:
+    -On success the function returns -1.
+    -On error the function returns an index to where the error was found
+     (0 is the first character, 1 the second, etc). If the error was not
+     a parsing error returns an index to the end of the string + 1.
+
+    Example: parser.Parse("3*x+y", "x,y");
+
+
+---------------------------------------------------------------------------
+const char* ErrorMsg(void) const;
+---------------------------------------------------------------------------
+
+    Returns a pointer to an error message string corresponding to the error
+    caused by Parse() (you can use this to print the proper error message to
+    the user). If no such error has occurred, returns 0.
+
+
+---------------------------------------------------------------------------
+ParseErrorType GetParseErrorType() const;
+---------------------------------------------------------------------------
+
+    Returns the type of parse error which occurred.
+
+    This method can be used to get the error type if ErrorMsg() is not
+    enough for printing the error message. In other words, this can be
+    used for printing customized error messages (eg. in another language).
+    If the default error messages suffice, then this method doesn't need
+    to be called.
+
+    FunctionParser::ParseErrorType is an enumerated type inside the class
+    (ie. its values are accessed like "FunctionParser::SYNTAX_ERROR").
+
+    The possible values for FunctionParser::ParseErrorType are listed below,
+    along with their equivalent error message returned by the ErrorMsg()
+    method:
+
+FP_NO_ERROR        : If no error occurred in the previous call to Parse().
+SYNTAX_ERROR       : "Syntax error"
+MISM_PARENTH       : "Mismatched parenthesis"
+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"
+INVALID_VARS       : "Syntax error in parameter 'Vars' given to "
+                     "FunctionParser::Parse()"
+ILL_PARAMS_AMOUNT  : "Illegal number of parameters to function"
+PREMATURE_EOS      : "Syntax error: Premature end of string"
+EXPECT_PARENTH_FUNC: "Syntax error: Expecting ( after function"
+
+
+---------------------------------------------------------------------------
+double Eval(const double* Vars);
+---------------------------------------------------------------------------
+
+    Evaluates the function given to Parse().
+    The array given as parameter must contain the same amount of values as
+    the amount of variables given to Parse(). Each value corresponds to each
+    variable, in the same order.
+
+    Return values:
+    -On success returns the evaluated value of the function given to
+     Parse().
+    -On error (such as division by 0) the return value is unspecified,
+     probably 0.
+
+    Example:
+
+      double Vars[] = {1, -2.5};
+      double result = parser.Eval(Vars);
+
+
+---------------------------------------------------------------------------
+int EvalError(void) const;
+---------------------------------------------------------------------------
+
+    Used to test if the call to Eval() succeeded.
+
+    Return values:
+      If there was no error in the previous call to Eval(), returns 0,
+      else returns a positive value as follows:
+        1: division by zero
+        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)
+
+
+---------------------------------------------------------------------------
+void Optimize();
+---------------------------------------------------------------------------
+
+    This method can be called after calling the Parse() method. It will try
+    to simplify the internal bytecode so that it will evaluate faster (it
+    tries to reduce the amount of opcodes in the bytecode).
+
+      For example, the bytecode for the function "5+x*y-25*4/8" will be
+    reduced to a bytecode equivalent to the function "x*y-7.5" (the original
+    11 opcodes will be reduced to 5). Besides calculating constant expressions
+    (like in the example), it also performs other types of simplifications
+    with variable and function expressions.
+
+      This method is quite slow and the decision of whether to use it or
+    not should depend on the type of application. If a function is parsed
+    once and evaluated millions of times, then calling Optimize() may speed-up
+    noticeably. However, if there are tons of functions to parse and each one
+    is evaluated once or just a few times, then calling Optimize() will only
+    slow down the program.
+      Also, if the original function is expected to be optimal, then calling
+    Optimize() would be useless.
+
+      Note: Currently this method does not make any checks (like Eval() does)
+    and thus things like "1/0" will cause undefined behaviour. (On the other
+    hand, if such expression is given to the parser, Eval() will always give
+    an error code, no matter what the parameters.) If caching this type of
+    errors is important, a work-around is to call Eval() once before calling
+    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).
+
+    (If you are interested in seeing how this method optimizes the opcode,
+    you can call the PrintByteCode() method before and after the call to
+    Optimize() to see the difference.)
+
+
+---------------------------------------------------------------------------
+bool AddConstant(const std::string& name, double value);
+---------------------------------------------------------------------------
+
+    This method can be used to add constants to the parser. Syntactically
+    constants are identical to variables (ie. they follow the same naming
+    rules and they can be used in the function string in the same way as
+    variables), but internally constants are directly replaced with their
+    value at parse time.
+
+      Constants used by a function must be added before calling Parse()
+    for that function. Constants are preserved between Parse() calls in
+    the current FunctionParser instance, so they don't need to be added
+    but once. (If you use the same constant in several instances of
+    FunctionParser, you will need to add it to all the instances separately.)
+
+      Constants can be added at any time and the value of old constants can
+    be changed, but new additions and changes will only have effect the next
+    time Parse() is called. (That is, changing the value of a constant
+    after calling Parse() and before calling Eval() will have no effect.)
+
+      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);
+
+    Now for example parser.Parse("x*pi", "x"); will be identical to the
+    call parser.Parse("x*3.14159265", "x");
+
+
+---------------------------------------------------------------------------
+bool AddFunction(const std::string& name,
+                 double (*functionPtr)(const double*),
+                 unsigned paramsAmount);
+---------------------------------------------------------------------------
+
+    This method can be used to add new functions to the parser. For example,
+    if you would like to add a function "sqr(A)" which squares the value
+    of A, you can do it with this method (so that you don't need to touch
+    the source code of the parser).
+
+      The method takes three parameters:
+
+    - The name of the function. The name follows the same naming conventions
+      as variable names.
+
+    - A C++ function, which will be called when evaluating the function
+      string (if the user-given function is called there). The C++ function
+      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 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
+    reserved), else true. If the return value is false, nothing is added.
+
+    Example:
+    Suppose we have a C++ function like this:
+
+    double Square(const double* p)
+    {
+        return p[0]*p[0];
+    }
+
+    Now we can add this function to the parser like this:
+
+    parser.AddFunction("sqr", Square, 1);
+
+    parser.Parse("2*sqr(x)", "x");
+
+
+    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.
+
+
+---------------------------------------------------------------------------
+bool AddFunction(const std::string& name, FunctionParser&);
+---------------------------------------------------------------------------
+
+    This method is almost identical to the previous AddFunction(), but
+    instead of taking a C++ function, it takes another FunctionParser
+    instance.
+
+    There are some important restrictions on making a FunctionParser instance
+    call another:
+
+    - The FunctionParser instance given as parameter must be initialized
+      with a Parse() call before giving it as parameter. That is, if you
+      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
+      parameter must not change after it has been given to the AddFunction()
+      of another instance. Changing the number of parameters 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);
+
+
+---------------------------------------------------------------------------
+
+  Example program:
+
+#include "fparser.hh"
+#include <iostream>
+
+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;
+}
+
+
+
+=============================================================================
+  - The function string
+=============================================================================
+
+  The function string understood by the class is very similar to the C-syntax.
+  Arithmetic float expressions can be created from float literals, variables
+or functions using the following operators in this order of precedence:
+
+   ()             expressions in parentheses first
+   -A             unary minus
+   A^B            exponentiation (A raised to the power B)
+   A*B  A/B  A%B  multiplication, division and modulo
+   A+B  A-B       addition and subtraction
+   A=B  A<B  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 class supports these functions:
+
+  abs(A)    : Absolute value of A. If A is negative, returns -A otherwise
+              returns A.
+  acos(A)   : Arc-cosine of A. Returns the angle, measured in radians,
+              whose cosine is A.
+  acosh(A)  : Same as acos() but for hyperbolic cosine.
+  asin(A)   : Arc-sine of A. Returns the angle, measured in radians, whose
+              sine is A.
+  asinh(A)  : Same as asin() but for hyperbolic sine.
+  atan(A)   : Arc-tangent of (A). Returns the angle, measured in radians,
+              whose tangent is (A).
+  atan2(A,B): Arc-tangent of A/B. The two main differences to atan() is
+              that it will return the right angle depending on the signs of
+              A and B (atan() can only return values betwen -pi/2 and pi/2),
+              and that the return value of pi/2 and -pi/2 are possible.
+  atanh(A)  : Same as atan() but for hyperbolic tangent.
+  ceil(A)   : Ceiling of A. Returns the smallest integer greater than A.
+              Rounds up to the next higher integer.
+  cos(A)    : Cosine of A. Returns the cosine of the angle A, where A is
+              measured in radians.
+  cosh(A)   : Same as cos() but for hyperbolic cosine.
+  cot(A)    : Cotangent of A (equivalent to 1/tan(A)).
+  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
+              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
+              non-repeating value approximately equal to 2.71828182846.
+  floor(A)  : Floor of A. Returns the largest integer less than A. Rounds
+              down to the next lower integer.
+  if(A,B,C) : If int(A) differs from 0, the return value of this function is B,
+              else C. Only the parameter which needs to be evaluated is
+              evaluated, the other parameter is skipped; this makes it safe to
+              use eval() in them.
+  int(A)    : Rounds A to the closest integer. 0.5 is rounded to 1.
+  log(A)    : Natural (base e) logarithm of A.
+  log10(A)  : Base 10 logarithm of A.
+  max(A,B)  : If A>B, the result is A, else B.
+  min(A,B)  : If A<B, the result is A, else B.
+  sec(A)    : Secant of A (equivalent to 1/cos(A)).
+  sin(A)    : Sine of A. Returns the sine of the angle A, where A is
+              measured in radians.
+  sinh(A)   : Same as sin() but for hyperbolic sine.
+  sqrt(A)   : Square root of A. Returns the value whose square is A.
+  tan(A)    : Tangent of A. Returns the tangent of the angle A, where A
+              is measured in radians.
+  tanh(A)   : Same as tan() but for hyperbolic tangent.
+
+
+  Examples of function string understood by the class:
+
+  "1+2"
+  "x-1"
+  "-sin(sqrt(x^2+y^2))"
+  "sqrt(XCoord*XCoord + YCoord*YCoord)"
+
+  An example of a recursive function is the factorial function:
+
+  "if(n>1, n*eval(n-1), 1)"
+
+  Note that a recursive call has some overhead, which makes it a bit slower
+  than any other operation. It may be a good idea to avoid recursive functions
+  in very time-critical applications. Recursion also takes some memory, so
+  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).
+
+
+=============================================================================
+  - Contacting the author
+=============================================================================
+
+  Any comments, bug reports, etc. should be sent to warp@iki.fi
+
+
+=============================================================================
+  - The algorithm used in the library
+=============================================================================
+
+  The whole idea behind the algorithm is to convert the regular infix
+format (the regular syntax for mathematical operations in most languages,
+like C and the input of the library) to postfix format. The postfix format
+is also called stack arithmetic since an expression in postfix format
+can be evaluated using a stack and operating with the top of the stack.
+
+  For example:
+
+  infix    postfix
+  2+3      2 3 +
+  1+2+3    1 2 + 3 +
+  5*2+8/2  5 2 * 8 2 / +
+  (5+9)*3  5 9 + 3 *
+
+  The postfix notation should be read in this way:
+
+  Let's take for example the expression: 5 2 * 8 2 / +
+  - Put 5 on the stack
+  - Put 2 on the stack
+  - Multiply the two values on the top of the stack and put the result on
+    the stack (removing the two old values)
+  - Put 8 on the stack
+  - Put 2 on the stack
+  - Divide the two values on the top of the stack
+  - Add the two values on the top of the stack (which are in this case
+    the result of 5*2 and 8/2, that is, 10 and 4).
+
+  At the end there's only one value in the stack, and that value is the
+result of the expression.
+
+  Why stack arithmetic?
+
+  The last example above can give you a hint.
+  In infix format operators have precedence and we have to use parentheses to
+group operations with lower precedence to be calculated before operations
+with higher precedence.
+  This causes a problem when evaluating an infix expression, specially
+when converting it to byte code. For example in this kind of expression:
+    (x+1)/(y+2)
+we have to calculate first the two additions before we can calculate the
+division. We have to also keep counting parentheses, since there can be
+a countless amount of nested parentheses. This usually means that you
+have to do some type of recursion.
+
+  The most simple and efficient way of calculating this is to convert it
+to postfix notation.
+  The postfix notation has the advantage that you can make all operations
+in a straightforward way. You just evaluate the expression from left to
+right, applying each operation directly and that's it. There are no
+parentheses to worry about. You don't need recursion anywhere.
+  You have to keep a stack, of course, but that's extremely easily done.
+Also you just operate with the top of the stack, which makes it very easy.
+You never have to go deeper than 2 items in the stack.
+  And even better: Evaluating an expression in postfix format is never
+slower than in infix format. All the contrary, in many cases it's a lot
+faster (eg. because all parentheses are optimized away).
+  The above example could be expressed in postfix format:
+    x 1 + y 2 + /
+
+  The good thing about the postfix notation is also the fact that it can
+be extremely easily expressed in bytecode form.
+  You only need a byte value for each operation, for each variable and
+to push a constant to the stack.
+  Then you can interpret this bytecode straightforwardly. You just interpret
+it byte by byte, from the beginning to the end. You never have to go back,
+make loops or anything.
+
+  This is what makes byte-coded stack arithmetic so fast.
+
+
+
+=============================================================================
+  Usage license:
+=============================================================================
+
+Copyright © 2003 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
+("Library" below).
+  The reason for having two distinct usage licenses is to make the library
+compatible with the GPL license while still being usable in other non-GPL
+(even commercial) software.
+
+A) If the Software using the Library is distributed under the GPL license,
+   then the Library can be used under the GPL license as well.
+
+   The Library will be under the GPL license only when used with the
+   Software. If the Library is separated from the Software and used in
+   another different software under a different license, then the Library
+   will have the B) license below.
+
+   Exception to the above: If the Library is modified for the GPL Software,
+   then the Library cannot be used with the B) license without the express
+   permission of the author of the modifications. A modified library will
+   be under the GPL license by default. That is, only the original,
+   unmodified version of the Library can be taken to another software
+   with the B) license below.
+
+   The author of the Software should provide an URL to the original
+   version of the Library if the one used in the Software has been
+   modified. (http://iki.fi/warp/FunctionParser/)
+
+   This text file must be distributed in its original intact form along
+   with the sources of the Library. (Documentation about possible
+   modifications to the library should be put in a different text file.)
+
+B) If the Software using the Library is not distributed under the GPL
+   license but under any other license, then the following usage license
+   applies to the Library:
+
+  1. This library is free for non-commercial usage. You can do whatever you
+     like with it as long as you don't claim you made it yourself.
+
+  2. It is possible to use this library in a commercial program, but in this
+     case you MUST contact me first (warp@iki.fi) and ask express permission
+     for this. (Read explanation at the end of the file.)
+       If you are making a free program or a shareware program with just a
+     nominal price (5 US dollars or less), you don't have to ask for
+     permission.
+       In any case, I DON'T WANT MONEY for the usage of this library. It is
+     free, period.
+
+  3. You can make any modifications you want to it so that it conforms your
+     needs. If you make modifications to it, you have, of course, credits for
+     the modified parts.
+
+  4. If you use this library in your own program, you don't have to provide
+     the source code if you don't want to (ie. the source code of your program
+     or this library).
+       If you DO include the source code for this library, this text file
+     must be included in its original intact form.
+
+  5. If you distribute a program which uses this library, and specially if you
+     provide the source code, proper credits MUST be included. Trying to
+     obfuscate the fact that this library is not made by you or that it is
+     free is expressly prohibited. When crediting the usage of this library,
+     it's enough to include my name and email address, that is:
+     "Juha Nieminen (warp@iki.fi)". Also a URL to the library download page
+     would be nice, although not required. The official URL is:
+       http://iki.fi/warp/FunctionParser/
+
+  6. And the necessary "lawyer stuff":
+
+     The above copyright notice and this permission notice shall be
+     included in all copies or substantial portions of the Software.
+
+     The software is provided "as is", without warranty of any kind,
+     express or implied, including but not limited to the warranties of
+     merchantability, fitness for a particular purpose and noninfringement.
+     In no event shall the authors or copyright holders be liable for any
+     claim, damages or other liability, whether in an action of contract,
+     tort or otherwise, arising from, out of or in connection with the
+     software or the use or other dealings in the software.
+
+
+---  Explanation of the section 2 of the B) license above:
+
+  The section 2 tries to define "fair use" of the library in commercial
+programs.
+  "Fair use" of the library means that the program is not heavily dependent
+on the library, but the library only provides a minor secondary feature
+to the program.
+  "Heavily dependent" means that the program depends so much on the library
+that without it the functionality of the program would be seriously
+degraded or the program would even become completely non-functional.
+
+  In other words: If the program does not depend heavily on the library,
+that is, the library only provides a minor secondary feature which could
+be removed without the program being degraded in any considerable way,
+then it's OK to use the library in the commercial program.
+  If, however, the program depends so heavily on the library that
+removing it would make the program non-functional or degrade its
+functionality considerably, then it's NOT OK to use the library.
+
+  The ideology behind this is that it's not fair to use a free library
+as a base for a commercial program, but it's fair if the library is
+just a minor, unimportant extra.
+
+  If you are going to ask me for permission to use the library in a
+commercial program, please describe the feature which the library will
+be providing and how important it is to the program.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.