]> Shamusworld >> Repos - architektonas/blob - fparser/src/fparser.cpp
Refactored CAD tool bars to use predefined actions.
[architektonas] / fparser / src / fparser.cpp
1 //===============================
2 // Function parser v2.51 by Warp
3 //===============================
4
5 // Comment out the following line if your compiler supports the (non-standard)
6 // asinh, acosh and atanh functions and you want them to be supported. If
7 // you are not sure, just leave it (those function will then not be supported).
8 #define NO_ASINH
9
10
11 // Uncomment the following line to disable the eval() function if it could
12 // be too dangerous in the target application:
13 //#define DISABLE_EVAL
14
15
16 // Comment this line out if you are not going to use the optimizer and want
17 // a slightly smaller library. The Optimize() method can still be called,
18 // but it will not do anything.
19 // If you are unsure, just leave it. It won't slow down the other parts of
20 // the library.
21 #if !defined(_MSC_VER)
22 #define SUPPORT_OPTIMIZER
23 #endif
24
25
26 //============================================================================
27
28 #include "fparser.h"
29
30 #include <cstdlib>
31 #include <cstring>
32 #include <cctype>
33 #include <cmath>
34 #include <new>
35 #include <algorithm>
36
37 using namespace std;
38
39 #ifndef M_PI
40 #define M_PI 3.1415926535897932384626433832795
41 #endif
42
43 namespace
44 {
45 // The functions must be in alphabetical order:
46     enum OPCODE
47     {
48         cAbs, cAcos,
49 #ifndef NO_ASINH
50         cAcosh,
51 #endif
52         cAsin,
53 #ifndef NO_ASINH
54         cAsinh,
55 #endif
56         cAtan,
57         cAtan2,
58 #ifndef NO_ASINH
59         cAtanh,
60 #endif
61         cCeil, cCos, cCosh, cCot, cCsc,
62 #ifndef DISABLE_EVAL
63         cEval,
64 #endif
65         cExp, cFloor, cIf, cInt, cLog, cLog10, cMax, cMin,
66         cSec, cSin, cSinh, cSqrt, cTan, cTanh,
67
68 // These do not need any ordering:
69         cImmed, cJump,
70         cNeg, cAdd, cSub, cMul, cDiv, cMod, cPow,
71         cEqual, cLess, cGreater, cAnd, cOr,
72
73         cDeg, cRad,
74
75         cFCall, cPCall,
76
77 #ifdef SUPPORT_OPTIMIZER
78         cVar, cDup, cInv,
79 #endif
80
81         VarBegin
82     };
83
84     struct FuncDefinition
85     {
86         const char* name;
87         unsigned nameLength;
88         unsigned opcode;
89         unsigned params;
90
91         // This is basically strcmp(), but taking 'nameLength' as string
92         // length (not ending '\0'):
93         bool operator<(const FuncDefinition& rhs) const
94         {
95             for(unsigned i = 0; i < nameLength; ++i)
96             {
97                 if(i == rhs.nameLength) return false;
98                 const char c1 = name[i], c2 = rhs.name[i];
99                 if(c1 < c2) return true;
100                 if(c2 < c1) return false;
101             }
102             return nameLength < rhs.nameLength;
103         }
104     };
105
106
107 // This list must be in alphabetical order:
108     const FuncDefinition Functions[]=
109     {
110         { "abs", 3, cAbs, 1 },
111         { "acos", 4, cAcos, 1 },
112 #ifndef NO_ASINH
113         { "acosh", 5, cAcosh, 1 },
114 #endif
115         { "asin", 4, cAsin, 1 },
116 #ifndef NO_ASINH
117         { "asinh", 5, cAsinh, 1 },
118 #endif
119         { "atan", 4, cAtan, 1 },
120         { "atan2", 5, cAtan2, 2 },
121 #ifndef NO_ASINH
122         { "atanh", 5, cAtanh, 1 },
123 #endif
124         { "ceil", 4, cCeil, 1 },
125         { "cos", 3, cCos, 1 },
126         { "cosh", 4, cCosh, 1 },
127         { "cot", 3, cCot, 1 },
128         { "csc", 3, cCsc, 1 },
129 #ifndef DISABLE_EVAL
130         { "eval", 4, cEval, 0 },
131 #endif
132         { "exp", 3, cExp, 1 },
133         { "floor", 5, cFloor, 1 },
134         { "if", 2, cIf, 0 },
135         { "int", 3, cInt, 1 },
136         { "log", 3, cLog, 1 },
137         { "log10", 5, cLog10, 1 },
138         { "max", 3, cMax, 2 },
139         { "min", 3, cMin, 2 },
140         { "sec", 3, cSec, 1 },
141         { "sin", 3, cSin, 1 },
142         { "sinh", 4, cSinh, 1 },
143         { "sqrt", 4, cSqrt, 1 },
144         { "tan", 3, cTan, 1 },
145         { "tanh", 4, cTanh, 1 }
146     };
147
148     const unsigned FUNC_AMOUNT = sizeof(Functions)/sizeof(Functions[0]);
149
150
151     // Returns a pointer to the FuncDefinition instance which 'name' is
152     // the same as the one given by 'F'. If no such function name exists,
153     // returns 0.
154     inline const FuncDefinition* FindFunction(const char* F)
155     {
156         FuncDefinition func = { F, 0, 0, 0 };
157         while(isalnum(F[func.nameLength])) ++func.nameLength;
158         if(func.nameLength)
159         {
160             const FuncDefinition* found =
161                 lower_bound(Functions, Functions+FUNC_AMOUNT, func);
162             if(found == Functions+FUNC_AMOUNT || func < *found)
163                 return 0;
164             return found;
165         }
166         return 0;
167     }
168 };
169
170 //---------------------------------------------------------------------------
171 // Constructors and destructors
172 //---------------------------------------------------------------------------
173 //===========================================================================
174 FunctionParser::FunctionParser():
175     ParseErrorType(-1), EvalErrorType(0)
176 {}
177
178 FunctionParser::~FunctionParser()
179 {}
180
181 FunctionParser::CompiledCode::CompiledCode():
182     ByteCode(0), ByteCodeSize(0),
183     Immed(0), ImmedSize(0),
184     Stack(0), StackSize(0)
185 {}
186
187 FunctionParser::CompiledCode::~CompiledCode()
188 {
189     if(ByteCode) { delete[] ByteCode; ByteCode=0; }
190     if(Immed) { delete[] Immed; Immed=0; }
191     if(Stack) { delete[] Stack; Stack=0; }
192 }
193
194
195 //---------------------------------------------------------------------------
196 // Function parsing
197 //---------------------------------------------------------------------------
198 //===========================================================================
199 namespace
200 {
201     // Error messages returned by ErrorMsg():
202     const char* ParseErrorMessage[]=
203     {
204         "Syntax error",                             // 0
205         "Mismatched parenthesis",                   // 1
206         "Missing ')'",                              // 2
207         "Empty parentheses",                        // 3
208         "Syntax error: Operator expected",          // 4
209         "Not enough memory",                        // 5
210         "An unexpected error ocurred. Please make a full bug report "
211         "to warp@iki.fi",                           // 6
212         "Syntax error in parameter 'Vars' given to "
213         "FunctionParser::Parse()",                  // 7
214         "Illegal number of parameters to function", // 8
215         "Syntax error: Premature end of string",    // 9
216         "Syntax error: Expecting ( after function", // 10
217         ""
218     };
219
220
221     // Parse variables
222     bool ParseVars(const string& Vars, map<string, unsigned>& dest)
223     {
224         unsigned varNumber = VarBegin;
225         unsigned ind1 = 0, ind2;
226
227         while(ind1 < Vars.size())
228         {
229             if(!isalpha(Vars[ind1]) && Vars[ind1]!='_') return false;
230             for(ind2=ind1+1; ind2<Vars.size() && Vars[ind2]!=','; ++ind2)
231                 if(!isalnum(Vars[ind2]) && Vars[ind2]!='_') return false;
232             const string varName = Vars.substr(ind1, ind2-ind1);
233
234             if(dest.insert(make_pair(varName, varNumber++)).second == false)
235                 return false;
236
237             ind1 = ind2+1;
238         }
239         return true;
240     }
241 };
242
243 bool FunctionParser::isValidName(const std::string& name)
244 {
245     if(name.empty() || (!isalpha(name[0]) && name[0] != '_')) return false;
246     for(unsigned i=0; i<name.size(); ++i)
247         if(!isalnum(name[i]) && name[i] != '_') return false;
248
249     if(FindFunction(name.c_str())) return false;
250
251     return true;
252 }
253
254
255 // Constants:
256 bool FunctionParser::AddConstant(const string& name, double value)
257 {
258     if(isValidName(name))
259     {
260         const char* n = name.c_str();
261         if(FindVariable(n, FuncParserNames) != FuncParserNames.end() ||
262            FindVariable(n, FuncPtrNames) != FuncPtrNames.end())
263             return false;
264
265         Constants[name] = value;
266         return true;
267     }
268     return false;
269 }
270
271 // Function pointers
272 bool FunctionParser::AddFunction(const std::string& name,
273                                  FunctionPtr func, unsigned paramsAmount)
274 {
275     if(paramsAmount == 0) return false; // Currently must be at least one
276
277     if(isValidName(name))
278     {
279         const char* n = name.c_str();
280         if(FindVariable(n, FuncParserNames) != FuncParserNames.end() ||
281            FindConstant(n) != Constants.end())
282             return false;
283
284         FuncPtrNames[name] = FuncPtrs.size();
285         FuncPtrs.push_back(FuncPtrData(func, paramsAmount));
286         return true;
287     }
288     return false;
289 }
290
291 bool FunctionParser::checkRecursiveLinking(const FunctionParser* fp)
292 {
293     if(fp == this) return true;
294     for(unsigned i=0; i<fp->FuncParsers.size(); ++i)
295         if(checkRecursiveLinking(fp->FuncParsers[i])) return true;
296     return false;
297 }
298
299 bool FunctionParser::AddFunction(const std::string& name,
300                                  FunctionParser& parser)
301 {
302     if(parser.varAmount == 0) return false; // Currently must be at least one
303
304     if(isValidName(name))
305     {
306         const char* n = name.c_str();
307         if(FindVariable(n, FuncPtrNames) != FuncPtrNames.end() ||
308            FindConstant(n) != Constants.end())
309             return false;
310
311         if(checkRecursiveLinking(&parser)) return false;
312
313         FuncParserNames[name] = FuncParsers.size();
314         FuncParsers.push_back(&parser);
315         return true;
316     }
317     return false;
318 }
319
320
321
322 // Main parsing function
323 // ---------------------
324 int FunctionParser::Parse(const std::string& Function,
325                           const std::string& Vars,
326                           bool useDegrees)
327 {
328     Variables.clear();
329
330     if(!ParseVars(Vars, Variables))
331     {
332         ParseErrorType = 7;
333         return Function.size();
334     }
335     varAmount = Variables.size(); // this is for Eval()
336
337     const char* Func = Function.c_str();
338
339     ParseErrorType = -1;
340
341     int Result = CheckSyntax(Func);
342     if(Result>=0) return Result;
343
344     useDegreeConversion = useDegrees;
345     if(!Compile(Func)) return Function.size();
346
347     Variables.clear();
348
349     ParseErrorType = -1;
350     return -1;
351 }
352
353 namespace
354 {
355     // Is given char an operator?
356     inline bool IsOperator(int c)
357     {
358         return strchr("+-*/%^=<>&|,",c)!=NULL;
359     }
360
361     // skip whitespace
362     inline void sws(const char* F, int& Ind)
363     {
364         while(F[Ind] && F[Ind] == ' ') ++Ind;
365     }
366 };
367
368 // Returns an iterator to the variable with the same name as 'F', or to
369 // Variables.end() if no such variable exists:
370 inline FunctionParser::VarMap_t::const_iterator
371 FunctionParser::FindVariable(const char* F, const VarMap_t& vars)
372 {
373     if(vars.size())
374     {
375         unsigned ind = 0;
376         while(isalnum(F[ind]) || F[ind] == '_') ++ind;
377         if(ind)
378         {
379             string name(F, ind);
380             return vars.find(name);
381         }
382     }
383     return vars.end();
384 }
385
386 inline FunctionParser::ConstMap_t::const_iterator
387 FunctionParser::FindConstant(const char* F)
388 {
389     if(Constants.size())
390     {
391         unsigned ind = 0;
392         while(isalnum(F[ind]) || F[ind] == '_') ++ind;
393         if(ind)
394         {
395             string name(F, ind);
396             return Constants.find(name);
397         }
398     }
399     return Constants.end();
400 }
401
402 //---------------------------------------------------------------------------
403 // Check function string syntax
404 // ----------------------------
405 int FunctionParser::CheckSyntax(const char* Function)
406 {
407     int Ind=0, ParenthCnt=0, c;
408     char* Ptr;
409
410     while(true)
411     {
412         sws(Function, Ind);
413         c=Function[Ind];
414
415 // Check for valid operand (must appear)
416
417         // Check for leading -
418         if(c=='-') { sws(Function, ++Ind); c=Function[Ind]; }
419         if(c==0) { ParseErrorType=9; return Ind; }
420
421         // Check for math function
422         bool foundFunc = false;
423         const FuncDefinition* fptr = FindFunction(&Function[Ind]);
424         if(fptr)
425         {
426             Ind += fptr->nameLength;
427             foundFunc = true;
428         }
429         else
430         {
431             // Check for user-defined function
432             VarMap_t::const_iterator fIter =
433                 FindVariable(&Function[Ind], FuncPtrNames);
434             if(fIter != FuncPtrNames.end())
435             {
436                 Ind += fIter->first.size();
437                 foundFunc = true;
438             }
439             else
440             {
441                 VarMap_t::const_iterator pIter =
442                     FindVariable(&Function[Ind], FuncParserNames);
443                 if(pIter != FuncParserNames.end())
444                 {
445                     Ind += pIter->first.size();
446                     foundFunc = true;
447                 }
448             }
449         }
450
451         if(foundFunc)
452         {
453             sws(Function, Ind);
454             c = Function[Ind];
455             if(c!='(') { ParseErrorType=10; return Ind; }
456         }
457
458         // Check for opening parenthesis
459         if(c=='(')
460         {
461             ++ParenthCnt;
462             sws(Function, ++Ind);
463             if(Function[Ind]==')') { ParseErrorType=3; return Ind; }
464             continue;
465         }
466
467         // Check for number
468         if(isdigit(c) || (c=='.' && isdigit(Function[Ind+1])))
469         {
470             strtod(&Function[Ind], &Ptr);
471             Ind += int(Ptr-&Function[Ind]);
472             sws(Function, Ind);
473             c = Function[Ind];
474         }
475         else
476         { // Check for variable
477             VarMap_t::const_iterator vIter =
478                 FindVariable(&Function[Ind], Variables);
479             if(vIter != Variables.end())
480                 Ind += vIter->first.size();
481             else
482             {
483                 // Check for constant
484                 ConstMap_t::const_iterator cIter =
485                     FindConstant(&Function[Ind]);
486                 if(cIter != Constants.end())
487                     Ind += cIter->first.size();
488                 else
489                 { ParseErrorType=0; return Ind; }
490             }
491             sws(Function, Ind);
492             c = Function[Ind];
493         }
494
495         // Check for closing parenthesis
496         while(c==')')
497         {
498             if((--ParenthCnt)<0) { ParseErrorType=1; return Ind; }
499             sws(Function, ++Ind);
500             c=Function[Ind];
501         }
502
503 // If we get here, we have a legal operand and now a legal operator or
504 // end of string must follow
505
506         // Check for EOS
507         if(c==0) break; // The only way to end the checking loop without error
508         // Check for operator
509         if(!IsOperator(c)) { ParseErrorType=4; return Ind; }
510
511 // If we get here, we have an operand and an operator; the next loop will
512 // check for another operand (must appear)
513         ++Ind;
514     } // while
515
516     // Check that all opened parentheses are also closed
517     if(ParenthCnt>0) { ParseErrorType=2; return Ind; }
518
519 // The string is ok
520     ParseErrorType=-1;
521     return -1;
522 }
523
524
525 // Compile function string to bytecode
526 // -----------------------------------
527 bool FunctionParser::Compile(const char* Function)
528 {
529     if(Comp.ByteCode) { delete[] Comp.ByteCode; Comp.ByteCode=0; }
530     if(Comp.Immed) { delete[] Comp.Immed; Comp.Immed=0; }
531     if(Comp.Stack) { delete[] Comp.Stack; Comp.Stack=0; }
532
533     vector<unsigned> byteCode; byteCode.reserve(1024);
534     tempByteCode = &byteCode;
535
536     vector<double> immed; immed.reserve(1024);
537     tempImmed = &immed;
538
539     Comp.StackSize = Comp.StackPtr = 0;
540
541     CompileExpression(Function, 0);
542     if(ParseErrorType >= 0) return false;
543
544     Comp.ByteCodeSize = byteCode.size();
545     Comp.ImmedSize = immed.size();
546
547     if(Comp.ByteCodeSize)
548     {
549         Comp.ByteCode = new unsigned[Comp.ByteCodeSize];
550         memcpy(Comp.ByteCode, &byteCode[0],
551                sizeof(unsigned)*Comp.ByteCodeSize);
552     }
553     if(Comp.ImmedSize)
554     {
555         Comp.Immed = new double[Comp.ImmedSize];
556         memcpy(Comp.Immed, &immed[0],
557                sizeof(double)*Comp.ImmedSize);
558     }
559     if(Comp.StackSize)
560         Comp.Stack = new double[Comp.StackSize];
561
562     return true;
563 }
564
565
566 inline void FunctionParser::AddCompiledByte(unsigned c)
567 {
568     tempByteCode->push_back(c);
569 }
570
571 inline void FunctionParser::AddImmediate(double i)
572 {
573     tempImmed->push_back(i);
574 }
575
576 inline void FunctionParser::AddFunctionOpcode(unsigned opcode)
577 {
578     if(useDegreeConversion)
579         switch(opcode)
580         {
581           case cCos:
582           case cCosh:
583           case cCot:
584           case cCsc:
585           case cSec:
586           case cSin:
587           case cSinh:
588           case cTan:
589           case cTanh:
590               AddCompiledByte(cRad);
591         }
592
593     AddCompiledByte(opcode);
594
595     if(useDegreeConversion)
596         switch(opcode)
597         {
598           case cAcos:
599 #ifndef NO_ASINH
600           case cAcosh:
601           case cAsinh:
602           case cAtanh:
603 #endif
604           case cAsin:
605           case cAtan:
606           case cAtan2:
607               AddCompiledByte(cDeg);
608         }
609 }
610
611 // Compile if()
612 int FunctionParser::CompileIf(const char* F, int ind)
613 {
614     int ind2 = CompileExpression(F, ind, true); // condition
615     sws(F, ind2);
616     if(F[ind2] != ',') { ParseErrorType=8; return ind2; }
617     AddCompiledByte(cIf);
618     unsigned curByteCodeSize = tempByteCode->size();
619     AddCompiledByte(0); // Jump index; to be set later
620     AddCompiledByte(0); // Immed jump index; to be set later
621
622     --Comp.StackPtr;
623
624     ind2 = CompileExpression(F, ind2+1, true); // then
625     sws(F, ind2);
626     if(F[ind2] != ',') { ParseErrorType=8; return ind2; }
627     AddCompiledByte(cJump);
628     unsigned curByteCodeSize2 = tempByteCode->size();
629     unsigned curImmedSize2 = tempImmed->size();
630     AddCompiledByte(0); // Jump index; to be set later
631     AddCompiledByte(0); // Immed jump index; to be set later
632
633     --Comp.StackPtr;
634
635     ind2 = CompileExpression(F, ind2+1, true); // else
636     sws(F, ind2);
637     if(F[ind2] != ')') { ParseErrorType=8; return ind2; }
638
639     // Set jump indices
640     (*tempByteCode)[curByteCodeSize] = curByteCodeSize2+1;
641     (*tempByteCode)[curByteCodeSize+1] = curImmedSize2;
642     (*tempByteCode)[curByteCodeSize2] = tempByteCode->size()-1;
643     (*tempByteCode)[curByteCodeSize2+1] = tempImmed->size();
644
645     return ind2+1;
646 }
647
648 int FunctionParser::CompileFunctionParams(const char* F, int ind,
649                                           unsigned requiredParams)
650 {
651     unsigned curStackPtr = Comp.StackPtr;
652     int ind2 = CompileExpression(F, ind);
653
654     if(Comp.StackPtr != curStackPtr+requiredParams)
655     { ParseErrorType=8; return ind; }
656
657     Comp.StackPtr -= requiredParams - 1;
658     sws(F, ind2);
659     return ind2+1; // F[ind2] is ')'
660 }
661
662 // Compiles element
663 int FunctionParser::CompileElement(const char* F, int ind)
664 {
665     sws(F, ind);
666     char c = F[ind];
667
668     if(c == '(')
669     {
670         ind = CompileExpression(F, ind+1);
671         sws(F, ind);
672         return ind+1; // F[ind] is ')'
673     }
674     else if(c == '-')
675     {
676         char c2 = F[ind+1];
677         if(!isdigit(c2) && c2!='.')
678         {
679             int ind2 = CompileElement(F, ind+1);
680             AddCompiledByte(cNeg);
681             return ind2;
682         }
683     }
684
685     if(isdigit(c) || c=='.' || c=='-') // Number
686     {
687         const char* startPtr = &F[ind];
688         char* endPtr;
689         double val = strtod(startPtr, &endPtr);
690         AddImmediate(val);
691         AddCompiledByte(cImmed);
692         ++Comp.StackPtr; if(Comp.StackPtr>Comp.StackSize) Comp.StackSize++;
693         return ind+(endPtr-startPtr);
694     }
695
696     if(isalpha(c) || c == '_') // Function, variable or constant
697     {
698         const FuncDefinition* func = FindFunction(F+ind);
699         if(func) // is function
700         {
701             int ind2 = ind + func->nameLength;
702             sws(F, ind2); // F[ind2] is '('
703             if(strcmp(func->name, "if") == 0) // "if" is a special case
704             {
705                 return CompileIf(F, ind2+1);
706             }
707
708 #ifndef DISABLE_EVAL
709             unsigned requiredParams =
710                 strcmp(func->name, "eval") == 0 ?
711                 Variables.size() : func->params;
712 #else
713             unsigned requiredParams = func->params;
714 #endif
715             ind2 = CompileFunctionParams(F, ind2+1, requiredParams);
716             AddFunctionOpcode(func->opcode);
717             return ind2; // F[ind2-1] is ')'
718         }
719
720         VarMap_t::const_iterator vIter = FindVariable(F+ind, Variables);
721         if(vIter != Variables.end()) // is variable
722         {
723             AddCompiledByte(vIter->second);
724             ++Comp.StackPtr; if(Comp.StackPtr>Comp.StackSize) Comp.StackSize++;
725             return ind + vIter->first.size();
726         }
727
728         ConstMap_t::const_iterator cIter = FindConstant(F+ind);
729         if(cIter != Constants.end()) // is constant
730         {
731             AddImmediate(cIter->second);
732             AddCompiledByte(cImmed);
733             ++Comp.StackPtr; if(Comp.StackPtr>Comp.StackSize) Comp.StackSize++;
734             return ind + cIter->first.size();
735         }
736
737         VarMap_t::const_iterator fIter = FindVariable(F+ind, FuncPtrNames);
738         if(fIter != FuncPtrNames.end()) // is user-defined function pointer
739         {
740             unsigned index = fIter->second;
741
742             int ind2 = ind + fIter->first.length();
743             sws(F, ind2); // F[ind2] is '('
744
745             ind2 = CompileFunctionParams(F, ind2+1, FuncPtrs[index].params);
746
747             AddCompiledByte(cFCall);
748             AddCompiledByte(index);
749             return ind2;
750         }
751
752         VarMap_t::const_iterator pIter = FindVariable(F+ind, FuncParserNames);
753         if(pIter != FuncParserNames.end()) // is user-defined function parser
754         {
755             unsigned index = pIter->second;
756
757             int ind2 = ind + pIter->first.length();
758             sws(F, ind2); // F[ind2] is '('
759
760             ind2 = CompileFunctionParams(F, ind2+1,
761                                          FuncParsers[index]->varAmount);
762
763             AddCompiledByte(cPCall);
764             AddCompiledByte(index);
765             return ind2;
766         }
767     }
768
769     ParseErrorType = 6;
770     return ind;
771 }
772
773 // Compiles '^'
774 int FunctionParser::CompilePow(const char* F, int ind)
775 {
776     int ind2 = CompileElement(F, ind);
777     sws(F, ind2);
778
779     while(F[ind2] == '^')
780     {
781         ind2 = CompileElement(F, ind2+1);
782         sws(F, ind2);
783         AddCompiledByte(cPow);
784         --Comp.StackPtr;
785     }
786
787     return ind2;
788 }
789
790 // Compiles '*', '/' and '%'
791 int FunctionParser::CompileMult(const char* F, int ind)
792 {
793     int ind2 = CompilePow(F, ind);
794     sws(F, ind2);
795     char op;
796
797     while((op = F[ind2]) == '*' || op == '/' || op == '%')
798     {
799         ind2 = CompilePow(F, ind2+1);
800         sws(F, ind2);
801         switch(op)
802         {
803           case '*': AddCompiledByte(cMul); break;
804           case '/': AddCompiledByte(cDiv); break;
805           case '%': AddCompiledByte(cMod); break;
806         }
807         --Comp.StackPtr;
808     }
809
810     return ind2;
811 }
812
813 // Compiles '+' and '-'
814 int FunctionParser::CompileAddition(const char* F, int ind)
815 {
816     int ind2 = CompileMult(F, ind);
817     sws(F, ind2);
818     char op;
819
820     while((op = F[ind2]) == '+' || op == '-')
821     {
822         ind2 = CompileMult(F, ind2+1);
823         sws(F, ind2);
824         AddCompiledByte(op=='+' ? cAdd : cSub);
825         --Comp.StackPtr;
826     }
827
828     return ind2;
829 }
830
831 // Compiles '=', '<' and '>'
832 int FunctionParser::CompileComparison(const char* F, int ind)
833 {
834     int ind2 = CompileAddition(F, ind);
835     sws(F, ind2);
836     char op;
837
838     while((op = F[ind2]) == '=' || op == '<' || op == '>')
839     {
840         ind2 = CompileAddition(F, ind2+1);
841         sws(F, ind2);
842         switch(op)
843         {
844           case '=': AddCompiledByte(cEqual); break;
845           case '<': AddCompiledByte(cLess); break;
846           case '>': AddCompiledByte(cGreater); break;
847         }
848         --Comp.StackPtr;
849     }
850
851     return ind2;
852 }
853
854 // Compiles '&'
855 int FunctionParser::CompileAnd(const char* F, int ind)
856 {
857     int ind2 = CompileComparison(F, ind);
858     sws(F, ind2);
859
860     while(F[ind2] == '&')
861     {
862         ind2 = CompileComparison(F, ind2+1);
863         sws(F, ind2);
864         AddCompiledByte(cAnd);
865         --Comp.StackPtr;
866     }
867
868     return ind2;
869 }
870
871 // Compiles '|'
872 int FunctionParser::CompileOr(const char* F, int ind)
873 {
874     int ind2 = CompileAnd(F, ind);
875     sws(F, ind2);
876
877     while(F[ind2] == '|')
878     {
879         ind2 = CompileAnd(F, ind2+1);
880         sws(F, ind2);
881         AddCompiledByte(cOr);
882         --Comp.StackPtr;
883     }
884
885     return ind2;
886 }
887
888 // Compiles ','
889 int FunctionParser::CompileExpression(const char* F, int ind, bool stopAtComma)
890 {
891     int ind2 = CompileOr(F, ind);
892     sws(F, ind2);
893
894     if(stopAtComma) return ind2;
895
896     while(F[ind2] == ',')
897     {
898         ind2 = CompileOr(F, ind2+1);
899         sws(F, ind2);
900     }
901
902     return ind2;
903 }
904
905
906 // Return parse error message
907 // --------------------------
908 const char* FunctionParser::ErrorMsg(void) const
909 {
910     if(ParseErrorType>=0) return ParseErrorMessage[ParseErrorType];
911     return 0;
912 }
913
914 //---------------------------------------------------------------------------
915 // Function evaluation
916 //---------------------------------------------------------------------------
917 //===========================================================================
918 namespace
919 {
920     inline int doubleToInt(double d)
921     {
922         return d<0 ? -int((-d)+.5) : int(d+.5);
923     }
924
925     inline double Min(double d1, double d2)
926     {
927         return d1<d2 ? d1 : d2;
928     }
929     inline double Max(double d1, double d2)
930     {
931         return d1>d2 ? d1 : d2;
932     }
933
934
935     inline double DegreesToRadians(double degrees)
936     {
937         return degrees*(M_PI/180.0);
938     }
939     inline double RadiansToDegrees(double radians)
940     {
941         return radians*(180.0/M_PI);
942     }
943 }
944
945 double FunctionParser::Eval(const double* Vars)
946 {
947     unsigned IP, DP=0;
948     int SP=-1;
949
950     for(IP=0; IP<Comp.ByteCodeSize; IP++)
951     {
952         switch(Comp.ByteCode[IP])
953         {
954 // Functions:
955           case   cAbs: Comp.Stack[SP]=fabs(Comp.Stack[SP]); break;
956           case  cAcos: if(Comp.Stack[SP]<-1 || Comp.Stack[SP]>1)
957                        { EvalErrorType=4; return 0; }
958                        Comp.Stack[SP]=acos(Comp.Stack[SP]); break;
959 #ifndef NO_ASINH
960           case cAcosh: Comp.Stack[SP]=acosh(Comp.Stack[SP]); break;
961 #endif
962           case  cAsin: if(Comp.Stack[SP]<-1 || Comp.Stack[SP]>1)
963                        { EvalErrorType=4; return 0; }
964                        Comp.Stack[SP]=asin(Comp.Stack[SP]); break;
965 #ifndef NO_ASINH
966           case cAsinh: Comp.Stack[SP]=asinh(Comp.Stack[SP]); break;
967 #endif
968           case  cAtan: Comp.Stack[SP]=atan(Comp.Stack[SP]); break;
969           case cAtan2: Comp.Stack[SP-1]=atan2(Comp.Stack[SP-1],Comp.Stack[SP]);
970                        SP--; break;
971 #ifndef NO_ASINH
972           case cAtanh: Comp.Stack[SP]=atanh(Comp.Stack[SP]); break;
973 #endif
974           case  cCeil: Comp.Stack[SP]=ceil(Comp.Stack[SP]); break;
975           case   cCos: Comp.Stack[SP]=cos(Comp.Stack[SP]); break;
976           case  cCosh: Comp.Stack[SP]=cosh(Comp.Stack[SP]); break;
977
978           case   cCot:
979               {
980                   double t = tan(Comp.Stack[SP]);
981                   if(t == 0) { EvalErrorType=1; return 0; }
982                   Comp.Stack[SP] = 1/t; break;
983               }
984           case   cCsc:
985               {
986                   double s = sin(Comp.Stack[SP]);
987                   if(s == 0) { EvalErrorType=1; return 0; }
988                   Comp.Stack[SP] = 1/s; break;
989               }
990
991
992 #ifndef DISABLE_EVAL
993           case  cEval:
994               {
995                   double* tmpStack = Comp.Stack;
996                   Comp.Stack = new double[Comp.StackSize];
997                   double retVal = Eval(&tmpStack[SP-varAmount+1]);
998                   delete[] Comp.Stack;
999                   Comp.Stack = tmpStack;
1000                   SP -= varAmount-1;
1001                   Comp.Stack[SP] = retVal;
1002                   break;
1003               }
1004 #endif
1005
1006           case   cExp: Comp.Stack[SP]=exp(Comp.Stack[SP]); break;
1007           case cFloor: Comp.Stack[SP]=floor(Comp.Stack[SP]); break;
1008
1009           case    cIf:
1010               {
1011                   unsigned jumpAddr = Comp.ByteCode[++IP];
1012                   unsigned immedAddr = Comp.ByteCode[++IP];
1013                   if(doubleToInt(Comp.Stack[SP]) == 0)
1014                   {
1015                       IP = jumpAddr;
1016                       DP = immedAddr;
1017                   }
1018                   SP--; break;
1019               }
1020
1021           case   cInt: Comp.Stack[SP]=floor(Comp.Stack[SP]+.5); break;
1022           case   cLog: if(Comp.Stack[SP]<=0) { EvalErrorType=3; return 0; }
1023                        Comp.Stack[SP]=log(Comp.Stack[SP]); break;
1024           case cLog10: if(Comp.Stack[SP]<=0) { EvalErrorType=3; return 0; }
1025                        Comp.Stack[SP]=log10(Comp.Stack[SP]); break;
1026           case   cMax: Comp.Stack[SP-1]=Max(Comp.Stack[SP-1],Comp.Stack[SP]);
1027                        SP--; break;
1028           case   cMin: Comp.Stack[SP-1]=Min(Comp.Stack[SP-1],Comp.Stack[SP]);
1029                        SP--; break;
1030           case   cSec:
1031               {
1032                   double c = cos(Comp.Stack[SP]);
1033                   if(c == 0) { EvalErrorType=1; return 0; }
1034                   Comp.Stack[SP] = 1/c; break;
1035               }
1036           case   cSin: Comp.Stack[SP]=sin(Comp.Stack[SP]); break;
1037           case  cSinh: Comp.Stack[SP]=sinh(Comp.Stack[SP]); break;
1038           case  cSqrt: if(Comp.Stack[SP]<0) { EvalErrorType=2; return 0; }
1039                        Comp.Stack[SP]=sqrt(Comp.Stack[SP]); break;
1040           case   cTan: Comp.Stack[SP]=tan(Comp.Stack[SP]); break;
1041           case  cTanh: Comp.Stack[SP]=tanh(Comp.Stack[SP]); break;
1042
1043
1044 // Misc:
1045           case cImmed: Comp.Stack[++SP]=Comp.Immed[DP++]; break;
1046           case  cJump: DP = Comp.ByteCode[IP+2];
1047                        IP = Comp.ByteCode[IP+1];
1048                        break;
1049
1050 // Operators:
1051           case   cNeg: Comp.Stack[SP]=-Comp.Stack[SP]; break;
1052           case   cAdd: Comp.Stack[SP-1]+=Comp.Stack[SP]; SP--; break;
1053           case   cSub: Comp.Stack[SP-1]-=Comp.Stack[SP]; SP--; break;
1054           case   cMul: Comp.Stack[SP-1]*=Comp.Stack[SP]; SP--; break;
1055           case   cDiv: if(Comp.Stack[SP]==0) { EvalErrorType=1; return 0; }
1056                        Comp.Stack[SP-1]/=Comp.Stack[SP]; SP--; break;
1057           case   cMod: if(Comp.Stack[SP]==0) { EvalErrorType=1; return 0; }
1058                        Comp.Stack[SP-1]=fmod(Comp.Stack[SP-1],Comp.Stack[SP]);
1059                        SP--; break;
1060           case   cPow: Comp.Stack[SP-1]=pow(Comp.Stack[SP-1],Comp.Stack[SP]);
1061                        SP--; break;
1062
1063           case cEqual: Comp.Stack[SP-1] = (Comp.Stack[SP-1]==Comp.Stack[SP]);
1064                        SP--; break;
1065           case  cLess: Comp.Stack[SP-1] = (Comp.Stack[SP-1]<Comp.Stack[SP]);
1066                        SP--; break;
1067           case cGreater: Comp.Stack[SP-1] = (Comp.Stack[SP-1]>Comp.Stack[SP]);
1068                          SP--; break;
1069           case   cAnd: Comp.Stack[SP-1] =
1070                            (doubleToInt(Comp.Stack[SP-1]) &&
1071                             doubleToInt(Comp.Stack[SP]));
1072                        SP--; break;
1073           case    cOr: Comp.Stack[SP-1] =
1074                            (doubleToInt(Comp.Stack[SP-1]) ||
1075                             doubleToInt(Comp.Stack[SP]));
1076                        SP--; break;
1077
1078 // Degrees-radians conversion:
1079           case   cDeg: Comp.Stack[SP]=RadiansToDegrees(Comp.Stack[SP]); break;
1080           case   cRad: Comp.Stack[SP]=DegreesToRadians(Comp.Stack[SP]); break;
1081
1082 // User-defined function calls:
1083           case cFCall:
1084               {
1085                   unsigned index = Comp.ByteCode[++IP];
1086                   unsigned params = FuncPtrs[index].params;
1087                   double retVal =
1088                       FuncPtrs[index].ptr(&Comp.Stack[SP-params+1]);
1089                   SP -= params-1;
1090                   Comp.Stack[SP] = retVal;
1091                   break;
1092               }
1093
1094           case cPCall:
1095               {
1096                   unsigned index = Comp.ByteCode[++IP];
1097                   unsigned params = FuncParsers[index]->varAmount;
1098                   double retVal =
1099                       FuncParsers[index]->Eval(&Comp.Stack[SP-params+1]);
1100                   SP -= params-1;
1101                   Comp.Stack[SP] = retVal;
1102                   break;
1103               }
1104
1105
1106 #ifdef SUPPORT_OPTIMIZER
1107           case   cVar: break; // Paranoia. These should never exist
1108           case   cDup: Comp.Stack[SP+1]=Comp.Stack[SP]; ++SP; break;
1109           case   cInv:
1110               if(Comp.Stack[SP]==0.0) { EvalErrorType=1; return 0; }
1111               Comp.Stack[SP]=1.0/Comp.Stack[SP];
1112               break;
1113 #endif
1114
1115 // Variables:
1116           default:
1117               Comp.Stack[++SP]=Vars[Comp.ByteCode[IP]-VarBegin];
1118         }
1119     }
1120
1121     EvalErrorType=0;
1122     return Comp.Stack[SP];
1123 }
1124
1125
1126 void FunctionParser::PrintByteCode(std::ostream& dest) const
1127 {
1128     for(unsigned IP=0, DP=0; IP<Comp.ByteCodeSize; IP++)
1129     {
1130         dest.width(8); dest.fill('0'); hex(dest); //uppercase(dest);
1131         dest << IP << ": ";
1132
1133         unsigned opcode = Comp.ByteCode[IP];
1134
1135         switch(opcode)
1136         {
1137           case cIf:
1138               dest << "jz\t";
1139               dest.width(8); dest.fill('0'); hex(dest); //uppercase(dest);
1140               dest << Comp.ByteCode[IP+1]+1 << endl;
1141               IP += 2;
1142               break;
1143
1144           case cJump:
1145               dest << "jump\t";
1146               dest.width(8); dest.fill('0'); hex(dest); //uppercase(dest);
1147               dest << Comp.ByteCode[IP+1]+1 << endl;
1148               IP += 2;
1149               break;
1150           case cImmed:
1151               dest.precision(10);
1152               dest << "push\t" << Comp.Immed[DP++] << endl;
1153               break;
1154
1155           case cFCall:
1156               {
1157                   unsigned index = Comp.ByteCode[++IP];
1158                   VarMap_t::const_iterator iter = FuncPtrNames.begin();
1159                   while(iter->second != index) ++iter;
1160                   dest << "call\t" << iter->first << endl;
1161                   break;
1162               }
1163
1164           case cPCall:
1165               {
1166                   unsigned index = Comp.ByteCode[++IP];
1167                   VarMap_t::const_iterator iter = FuncParserNames.begin();
1168                   while(iter->second != index) ++iter;
1169                   dest << "call\t" << iter->first << endl;
1170                   break;
1171               }
1172
1173           default:
1174               if(opcode < VarBegin)
1175               {
1176                   string n;
1177                   switch(opcode)
1178                   {
1179                     case cNeg: n = "neg"; break;
1180                     case cAdd: n = "add"; break;
1181                     case cSub: n = "sub"; break;
1182                     case cMul: n = "mul"; break;
1183                     case cDiv: n = "div"; break;
1184                     case cMod: n = "mod"; break;
1185                     case cPow: n = "pow"; break;
1186                     case cEqual: n = "eq"; break;
1187                     case cLess: n = "lt"; break;
1188                     case cGreater: n = "gt"; break;
1189                     case cAnd: n = "and"; break;
1190                     case cOr: n = "or"; break;
1191                     case cDeg: n = "deg"; break;
1192                     case cRad: n = "rad"; break;
1193
1194 #ifndef DISABLE_EVAL
1195                     case cEval: n = "call\t0"; break;
1196 #endif
1197
1198 #ifdef SUPPORT_OPTIMIZER
1199                     case cVar: n = "(var)"; break;
1200                     case cDup: n = "dup"; break;
1201                     case cInv: n = "inv"; break;
1202 #endif
1203
1204                     default: n = Functions[opcode-cAbs].name;
1205                   }
1206                   dest << n << endl;
1207               }
1208               else
1209               {
1210                   dest << "push\tVar" << opcode-VarBegin << endl;
1211               }
1212         }
1213     }
1214 }
1215
1216
1217
1218 //========================================================================
1219 // Optimization code was contributed by Bisqwit (http://iki.fi/bisqwit/)
1220 //========================================================================
1221 #ifdef SUPPORT_OPTIMIZER
1222
1223 #include <list>
1224 #include <utility>
1225
1226 #define CONSTANT_E     2.71828182845904509080  // exp(1)
1227 #define CONSTANT_PI    M_PI                    // atan2(0,-1)
1228 #define CONSTANT_L10   2.30258509299404590109  // log(10)
1229 #define CONSTANT_L10I  0.43429448190325176116  // 1/log(10)
1230 #define CONSTANT_L10E  CONSTANT_L10I           // log10(e)
1231 #define CONSTANT_L10EI CONSTANT_L10            // 1/log10(e)
1232 #define CONSTANT_DR    (180.0 / M_PI)          // 180/pi
1233 #define CONSTANT_RD    (M_PI / 180.0)          // pi/180
1234
1235 class compres
1236 {
1237     // states: 0=false, 1=true, 2=unknown
1238 public:
1239     compres(bool b) : state(b) {}
1240     compres(char v) : state(v) {}
1241     // is it?
1242     operator bool() const { return state != 0; }
1243     // is it not?
1244     bool operator! () const { return state != 1; }
1245     bool operator==(bool b) const { return state != !b; }
1246     bool operator!=(bool b) const { return state != b; }
1247 private:
1248     char state;
1249 };
1250
1251 namespace {
1252 const compres maybe = (char)2;
1253 }
1254
1255 class SubTree
1256 {
1257     struct CodeTree *tree;
1258     bool sign;  // Only possible when parent is cAdd or cMul
1259
1260     inline void flipsign() { sign = !sign; }
1261 public:
1262     SubTree();
1263     SubTree(double value);
1264     SubTree(const SubTree &b);
1265     SubTree(const struct CodeTree &b);
1266
1267     ~SubTree();
1268     const SubTree &operator= (const SubTree &b);
1269     const SubTree &operator= (const CodeTree &b);
1270
1271     bool getsign() const { return sign; }
1272
1273     const struct CodeTree* operator-> () const { return tree; }
1274     const struct CodeTree& operator* () const { return *tree; }
1275     struct CodeTree* operator-> () { return tree; }
1276     struct CodeTree& operator* () { return *tree; }
1277
1278     bool operator< (const SubTree& b) const;
1279     bool operator== (const SubTree& b) const;
1280     void Negate(); // Note: Parent must be cAdd
1281     void Invert(); // Note: Parent must be cMul
1282
1283     void CheckConstNeg();
1284     void CheckConstInv();
1285 };
1286
1287 namespace {
1288 bool IsNegate(const SubTree &p1, const SubTree &p2);
1289 bool IsInverse(const SubTree &p1, const SubTree &p2);
1290 }
1291
1292 typedef list<SubTree> paramlist;
1293
1294 struct CodeTreeData
1295 {
1296     paramlist args;
1297
1298 private:
1299     unsigned op;       // Operation
1300     double value;      // In case of cImmed
1301     unsigned var;      // In case of cVar
1302     unsigned funcno;   // In case of cFCall, cPCall
1303
1304 public:
1305     CodeTreeData() : op(cAdd) {}
1306     ~CodeTreeData() {}
1307
1308     void SetOp(unsigned newop)     { op=newop; }
1309     void SetFuncNo(unsigned newno) { funcno=newno; }
1310     unsigned GetFuncNo() const { return funcno; }
1311
1312     bool IsFunc() const  { return op == cFCall || op == cPCall; }
1313     bool IsImmed() const { return op == cImmed; }
1314     bool IsVar() const   { return op == cVar; }
1315     inline unsigned GetOp() const { return op; }
1316     inline double GetImmed() const
1317     {
1318         return value;
1319     }
1320     inline unsigned GetVar() const
1321     {
1322         return var;
1323     }
1324
1325     void AddParam(const SubTree &p)
1326     {
1327         args.push_back(p);
1328     }
1329     void SetVar(unsigned v)
1330     {
1331         args.clear();
1332         op  = cVar;
1333         var = v;
1334     }
1335     void SetImmed(double v)
1336     {
1337         args.clear();
1338         op       = cImmed;
1339         value    = orig = v;
1340         inverted = negated = false;
1341     }
1342     void NegateImmed()
1343     {
1344         negated = !negated;
1345         UpdateValue();
1346     }
1347     void InvertImmed()
1348     {
1349         inverted = !inverted;
1350         UpdateValue();
1351     }
1352
1353     bool IsOriginal() const { return !(IsInverted() || IsNegated()); }
1354     bool IsInverted() const { return inverted; }
1355     bool IsNegated() const { return negated; }
1356     bool IsInvertedOriginal() const { return IsInverted() && !IsNegated(); }
1357     bool IsNegatedOriginal() const { return !IsInverted() && IsNegated(); }
1358
1359 private:
1360     void UpdateValue()
1361     {
1362         value = orig;
1363         if(IsInverted()) { value = 1.0 / value;
1364                            // FIXME: potential divide by zero.
1365                          }
1366         if(IsNegated()) value = -value;
1367     }
1368
1369     double orig;
1370     bool inverted;
1371     bool negated;
1372 protected:
1373     // Ensure we don't accidentally copy this
1374     void operator=(const CodeTreeData &b);
1375 };
1376
1377 class CodeTreeDataPtr
1378 {
1379     typedef pair<CodeTreeData, unsigned> p_t;
1380     typedef p_t* pp;
1381     mutable pp p;
1382
1383     void Alloc()   const { ++p->second; }
1384     void Dealloc() const { if(!--p->second) delete p; p = 0; }
1385
1386     void PrepareForWrite()
1387     {
1388         // We're ready if we're the only owner.
1389         if(p->second == 1) return;
1390
1391         // Then make a clone.
1392         p_t *newtree = new p_t(p->first, 1);
1393         // Forget the old
1394         Dealloc();
1395         // Keep the new
1396         p = newtree;
1397     }
1398
1399 public:
1400     CodeTreeDataPtr() : p(new p_t) { p->second = 1; }
1401     CodeTreeDataPtr(const CodeTreeDataPtr &b): p(b.p) { Alloc(); }
1402     ~CodeTreeDataPtr() { Dealloc(); }
1403     const CodeTreeDataPtr &operator= (const CodeTreeDataPtr &b)
1404     {
1405         b.Alloc();
1406         Dealloc();
1407         p = b.p;
1408         return *this;
1409     }
1410     const CodeTreeData *operator-> () const { return &p->first; }
1411     const CodeTreeData &operator*  () const { return p->first; }
1412     CodeTreeData *operator-> () { PrepareForWrite(); return &p->first; }
1413     CodeTreeData &operator*  () { PrepareForWrite(); return p->first; }
1414
1415     void Shock();
1416 };
1417
1418 #define CHECKCONSTNEG(item, op) \
1419     ((op)==cMul) \
1420        ? (item).CheckConstInv() \
1421        : (item).CheckConstNeg()
1422
1423 struct CodeTree
1424 {
1425     CodeTreeDataPtr data;
1426
1427 private:
1428     typedef paramlist::iterator pit;
1429     typedef paramlist::const_iterator pcit;
1430
1431     template<unsigned v> inline void chk() const
1432     {
1433     }
1434
1435 public:
1436     const pcit GetBegin() const { return data->args.begin(); }
1437     const pcit GetEnd()   const { return data->args.end(); }
1438     const pit GetBegin() { return data->args.begin(); }
1439     const pit GetEnd()   { return data->args.end(); }
1440     const SubTree& getp0() const { chk<1>();pcit tmp=GetBegin();               return *tmp; }
1441     const SubTree& getp1() const { chk<2>();pcit tmp=GetBegin(); ++tmp;        return *tmp; }
1442     const SubTree& getp2() const { chk<3>();pcit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; }
1443     unsigned GetArgCount() const { return data->args.size(); }
1444     void Erase(const pit p)      { data->args.erase(p); }
1445
1446     SubTree& getp0() { chk<1>();pit tmp=GetBegin();               return *tmp; }
1447     SubTree& getp1() { chk<2>();pit tmp=GetBegin(); ++tmp;        return *tmp; }
1448     SubTree& getp2() { chk<3>();pit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; }
1449
1450     // set
1451     void SetImmed(double v) { data->SetImmed(v); }
1452     void SetOp(unsigned op) { data->SetOp(op); }
1453     void SetVar(unsigned v) { data->SetVar(v); }
1454     // get
1455     double GetImmed() const { return data->GetImmed(); }
1456     unsigned GetVar() const { return data->GetVar(); }
1457     unsigned GetOp() const  { return data->GetOp(); }
1458     // test
1459     bool IsImmed() const { return data->IsImmed(); }
1460     bool IsVar()   const { return data->IsVar(); }
1461     // act
1462     void AddParam(const SubTree &p) { data->AddParam(p); }
1463     void NegateImmed() { data->NegateImmed(); } // don't use when op!=cImmed
1464     void InvertImmed() { data->InvertImmed(); } // don't use when op!=cImmed
1465
1466     compres NonZero() const { if(!IsImmed()) return maybe;
1467                               return GetImmed() != 0.0; }
1468     compres IsPositive() const { if(!IsImmed()) return maybe;
1469                                  return GetImmed() > 0.0; }
1470
1471 private:
1472     struct ConstList
1473     {
1474         double voidvalue;
1475         list<pit> cp;
1476         double value;
1477         unsigned size() const { return cp.size(); }
1478     };
1479     struct ConstList BuildConstList();
1480     void KillConst(const ConstList &cl)
1481     {
1482         for(list<pit>::const_iterator i=cl.cp.begin(); i!=cl.cp.end(); ++i)
1483             Erase(*i);
1484     }
1485     void FinishConst(const ConstList &cl)
1486     {
1487         if(cl.value != cl.voidvalue && cl.size() > 1) AddParam(cl.value);
1488         if(cl.value == cl.voidvalue || cl.size() > 1) KillConst(cl);
1489     }
1490
1491 public:
1492     CodeTree() {}
1493     CodeTree(double v) { SetImmed(v); }
1494
1495     CodeTree(unsigned op, const SubTree &p)
1496     {
1497         SetOp(op);
1498         AddParam(p);
1499     }
1500     CodeTree(unsigned op, const SubTree &p1, const SubTree &p2)
1501     {
1502         SetOp(op);
1503         AddParam(p1);
1504         AddParam(p2);
1505     }
1506
1507     bool operator== (const CodeTree& b) const;
1508     bool operator< (const CodeTree& b) const;
1509
1510 private:
1511     bool IsSortable() const
1512     {
1513         switch(GetOp())
1514         {
1515             case cAdd:  case cMul:
1516             case cEqual:
1517             case cAnd: case cOr:
1518             case cMax: case cMin:
1519                 return true;
1520             default:
1521                 return false;
1522         }
1523     }
1524     void SortIfPossible()
1525     {
1526         if(IsSortable())
1527         {
1528             data->args.sort();
1529         }
1530     }
1531
1532     void ReplaceWithConst(double value)
1533     {
1534         SetImmed(value);
1535
1536         /* REMEMBER TO CALL CheckConstInv / CheckConstNeg
1537          * FOR PARENT SubTree, OR MAYHEM HAPPENS
1538          */
1539     }
1540
1541     void ReplaceWith(const CodeTree &b)
1542     {
1543         // If b is child of *this, mayhem
1544         // happens. So we first make a clone
1545         // and then proceed with copy.
1546         CodeTreeDataPtr tmp = b.data;
1547         tmp.Shock();
1548         data = tmp;
1549     }
1550
1551     void ReplaceWith(unsigned op, const SubTree &p)
1552     {
1553         ReplaceWith(CodeTree(op, p));
1554     }
1555
1556     void ReplaceWith(unsigned op, const SubTree &p1, const SubTree &p2)
1557     {
1558         ReplaceWith(CodeTree(op, p1, p2));
1559     }
1560
1561     void OptimizeConflict()
1562     {
1563         // This optimization does this: x-x = 0, x/x = 1, a+b-a = b.
1564
1565         if(GetOp() == cAdd || GetOp() == cMul)
1566         {
1567         Redo:
1568             pit a, b;
1569             for(a=GetBegin(); a!=GetEnd(); ++a)
1570             {
1571                 for(b=GetBegin(); ++b != GetEnd(); )
1572                 {
1573                     const SubTree &p1 = *a;
1574                     const SubTree &p2 = *b;
1575
1576                     if(GetOp() == cMul ? IsInverse(p1,p2)
1577                                        : IsNegate(p1,p2))
1578                     {
1579                         // These parameters complement each others out
1580                         Erase(b);
1581                         Erase(a);
1582                         goto Redo;
1583                     }
1584                 }
1585             }
1586         }
1587         OptimizeRedundant();
1588     }
1589
1590     void OptimizeRedundant()
1591     {
1592         // This optimization does this: min()=0, max()=0, add()=0, mul()=1
1593
1594         if(!GetArgCount())
1595         {
1596             if(GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax)
1597                 ReplaceWithConst(0);
1598             else if(GetOp() == cMul)
1599                 ReplaceWithConst(1);
1600             return;
1601         }
1602
1603         // And this: mul(x) = x, min(x) = x, max(x) = x, add(x) = x
1604
1605         if(GetArgCount() == 1)
1606         {
1607             if(GetOp() == cMul || GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax)
1608                 if(!getp0().getsign())
1609                     ReplaceWith(*getp0());
1610         }
1611
1612         OptimizeDoubleNegations();
1613     }
1614
1615     void OptimizeDoubleNegations()
1616     {
1617         if(GetOp() == cAdd)
1618         {
1619             // Eschew double negations
1620
1621             // If any of the elements is cMul
1622             // and has a numeric constant, negate
1623             // the constant and negate sign.
1624
1625             for(pit a=GetBegin(); a!=GetEnd(); ++a)
1626             {
1627                 SubTree &pa = *a;
1628                 if(pa.getsign()
1629                 && pa->GetOp() == cMul)
1630                 {
1631                     CodeTree &p = *pa;
1632                     for(pit b=p.GetBegin();
1633                             b!=p.GetEnd(); ++b)
1634                     {
1635                         SubTree &pb = *b;
1636                         if(pb->IsImmed())
1637                         {
1638                             pb.Negate();
1639                             pa.Negate();
1640                             break;
1641                         }
1642                     }
1643                 }
1644             }
1645         }
1646
1647         if(GetOp() == cMul)
1648         {
1649             // If any of the elements is cPow
1650             // and has a numeric exponent, negate
1651             // the exponent and negate sign.
1652
1653             for(pit a=GetBegin(); a!=GetEnd(); ++a)
1654             {
1655                 SubTree &pa = *a;
1656                 if(pa.getsign() && pa->GetOp() == cPow)
1657                 {
1658                     CodeTree &p = *pa;
1659                     if(p.getp1()->IsImmed())
1660                     {
1661                         // negate ok for pow when op=cImmed
1662                         p.getp1().Negate();
1663                         pa.Negate();
1664                     }
1665                 }
1666             }
1667         }
1668     }
1669
1670     void OptimizeConstantMath1()
1671     {
1672         // This optimization does three things:
1673         //      - For adding groups:
1674         //          Constants are added together.
1675         //      - For multiplying groups:
1676         //          Constants are multiplied together.
1677         //      - For function calls:
1678         //          If all parameters are constants,
1679         //          the call is replaced with constant value.
1680
1681         // First, do this:
1682         OptimizeAddMulFlat();
1683
1684         switch(GetOp())
1685         {
1686             case cAdd:
1687             {
1688                 ConstList cl = BuildConstList();
1689                 FinishConst(cl);
1690                 break;
1691             }
1692             case cMul:
1693             {
1694                 ConstList cl = BuildConstList();
1695
1696                 if(cl.value == 0.0) ReplaceWithConst(0.0);
1697                 else FinishConst(cl);
1698
1699                 break;
1700             }
1701             #define ConstantUnaryFun(token, fun) \
1702                 case token: { const SubTree &p0 = getp0(); \
1703                     if(p0->IsImmed()) ReplaceWithConst(fun(p0->GetImmed())); \
1704                     break; }
1705             #define ConstantBinaryFun(token, fun) \
1706                 case token: { const SubTree &p0 = getp0(); \
1707                               const SubTree &p1 = getp1(); \
1708                     if(p0->IsImmed() && \
1709                        p1->IsImmed()) ReplaceWithConst(fun(p0->GetImmed(), p1->GetImmed())); \
1710                     break; }
1711
1712             // FIXME: potential invalid parameters for functions
1713             //        can cause exceptions here
1714
1715             ConstantUnaryFun(cAbs,   fabs);
1716             ConstantUnaryFun(cAcos,  acos);
1717             ConstantUnaryFun(cAsin,  asin);
1718             ConstantUnaryFun(cAtan,  atan);
1719             ConstantUnaryFun(cCeil,  ceil);
1720             ConstantUnaryFun(cCos,   cos);
1721             ConstantUnaryFun(cCosh,  cosh);
1722             ConstantUnaryFun(cFloor, floor);
1723             ConstantUnaryFun(cLog,   log);
1724             ConstantUnaryFun(cSin,   sin);
1725             ConstantUnaryFun(cSinh,  sinh);
1726             ConstantUnaryFun(cTan,   tan);
1727             ConstantUnaryFun(cTanh,  tanh);
1728             ConstantBinaryFun(cAtan2, atan2);
1729             ConstantBinaryFun(cMax,   Max);
1730             ConstantBinaryFun(cMin,   Min);
1731             ConstantBinaryFun(cMod,   fmod); // not a func, but belongs here too
1732             ConstantBinaryFun(cPow,   pow);
1733
1734             case cNeg:
1735             case cSub:
1736             case cDiv:
1737                 /* Unreached (nonexistent operator)
1738                  * TODO: internal error here?
1739                  */
1740                 break;
1741
1742             case cCot:
1743             case cCsc:
1744             case cSec:
1745             case cDeg:
1746             case cRad:
1747             case cLog10:
1748             case cSqrt:
1749             case cExp:
1750                 /* Unreached (nonexistent function)
1751                  * TODO: internal error here?
1752                  */
1753                  break;
1754         }
1755
1756         OptimizeConflict();
1757     }
1758
1759     void OptimizeAddMulFlat()
1760     {
1761         // This optimization flattens the topography of the tree.
1762         //   Examples:
1763         //       x + (y+z) = x+y+z
1764         //       x * (y/z) = x*y/z
1765         //       x / (y/z) = x/y*z
1766
1767         if(GetOp() == cAdd || GetOp() == cMul)
1768         {
1769             // If children are same type as parent add them here
1770             for(pit b, a=GetBegin(); a!=GetEnd(); a=b)
1771             {
1772                 const SubTree &pa = *a;  b=a; ++b;
1773                 if(pa->GetOp() != GetOp()) continue;
1774
1775                 // Child is same type
1776                 for(pcit c=pa->GetBegin();
1777                          c!=pa->GetEnd();
1778                          ++c)
1779                 {
1780                     const SubTree &pb = *c;
1781                     if(pa.getsign())
1782                     {
1783                         // +a -(+b +c)
1784                         // means b and c will be negated
1785
1786                         SubTree tmp = pb;
1787                         if(GetOp() == cMul)
1788                             tmp.Invert();
1789                         else
1790                             tmp.Negate();
1791                         AddParam(tmp);
1792                     }
1793                     else
1794                         AddParam(pb);
1795                 }
1796                 Erase(a);
1797
1798                 // Note: OptimizeConstantMath1() would be a good thing to call next.
1799             }
1800         }
1801     }
1802
1803     void OptimizeLinearCombine()
1804     {
1805         // This optimization does the following:
1806         //
1807         //   x*x*x*x -> x^4
1808         //   x+x+x+x -> x*4
1809         //   x*x     -> x^2
1810         //   x/z/z   ->
1811         //
1812
1813         // Remove conflicts first, so we don't have to worry about signs.
1814         OptimizeConflict();
1815
1816         bool didchanges = false;
1817         if(GetOp() == cAdd || GetOp() == cMul)
1818         {
1819         Redo:
1820             for(pit a=GetBegin(); a!=GetEnd(); ++a)
1821             {
1822                 const SubTree &pa = *a;
1823
1824                 list<pit> poslist;
1825
1826                 for(pit b=a; ++b!=GetEnd(); )
1827                 {
1828                     const SubTree &pb = *b;
1829                     if(*pa == *pb)
1830                         poslist.push_back(b);
1831                 }
1832
1833                 unsigned min = 2;
1834                 if(poslist.size() >= min)
1835                 {
1836                     SubTree arvo = pa;
1837                     bool negate = arvo.getsign();
1838
1839                     double factor = poslist.size() + 1;
1840
1841                     if(negate)
1842                     {
1843                         arvo.Negate();
1844                         factor = -factor;
1845                     }
1846
1847                     CodeTree tmp(GetOp()==cAdd ? cMul : cPow,
1848                                  arvo,
1849                                  factor);
1850
1851                     list<pit>::const_iterator j;
1852                     for(j=poslist.begin(); j!=poslist.end(); ++j)
1853                         Erase(*j);
1854                     poslist.clear();
1855
1856                     *a = tmp;
1857                     didchanges = true;
1858                     goto Redo;
1859                 }
1860             }
1861         }
1862         if(didchanges)
1863         {
1864             // As a result, there might be need for this:
1865             OptimizeAddMulFlat();
1866             // And this:
1867             OptimizeRedundant();
1868         }
1869     }
1870
1871     void OptimizeLogarithm()
1872     {
1873         /*
1874             This is basic logarithm math:
1875               pow(X,Y)/log(Y) = X
1876               log(X)/log(Y) = logY(X)
1877               log(X^Y)      = log(X)*Y
1878               log(X*Y)      = log(X)+log(Y)
1879               exp(log(X)*Y) = X^Y
1880
1881             This function does these optimizations:
1882                pow(const_E, log(x))   = x
1883                pow(const_E, log(x)*y) = x^y
1884                pow(10,      log(x)*const_L10I*y) = x^y
1885                pow(z,       log(x)/log(z)*y)     = x^y
1886
1887             And this:
1888                log(x^z)             = z * log(x)
1889             Which automatically causes these too:
1890                log(pow(const_E, x))         = x
1891                log(pow(y,       x))         = x * log(y)
1892                log(pow(pow(const_E, y), x)) = x*y
1893
1894             And it does this too:
1895                log(x) + log(y) + log(z) = log(x * y * z)
1896                log(x * exp(y)) = log(x) + y
1897
1898         */
1899
1900         // Must be already in exponential form.
1901
1902         // Optimize exponents before doing something.
1903         OptimizeExponents();
1904
1905         if(GetOp() == cLog)
1906         {
1907             // We should have one parameter for log() function.
1908             // If we don't, we're screwed.
1909
1910             const SubTree &p = getp0();
1911
1912             if(p->GetOp() == cPow)
1913             {
1914                 // Found log(x^y)
1915                 SubTree p0 = p->getp0(); // x
1916                 SubTree p1 = p->getp1(); // y
1917
1918                 // Build the new logarithm.
1919                 CodeTree tmp(GetOp(), p0);  // log(x)
1920
1921                 // Become log(x) * y
1922                 ReplaceWith(cMul, tmp, p1);
1923             }
1924             else if(p->GetOp() == cMul)
1925             {
1926                 // Redefine &p nonconst
1927                 SubTree &p = getp0();
1928
1929                 p->OptimizeAddMulFlat();
1930                 p->OptimizeExponents();
1931                 CHECKCONSTNEG(p, p->GetOp());
1932
1933                 list<SubTree> adds;
1934
1935                 for(pit b, a = p->GetBegin();
1936                            a != p->GetEnd(); a=b)
1937                 {
1938                     SubTree &pa = *a;  b=a; ++b;
1939                     if(pa->GetOp() == cPow
1940                     && pa->getp0()->IsImmed()
1941                     && pa->getp0()->GetImmed() == CONSTANT_E)
1942                     {
1943                         adds.push_back(pa->getp1());
1944                         p->Erase(a);
1945                         continue;
1946                     }
1947                 }
1948                 if(adds.size())
1949                 {
1950                     CodeTree tmp(cAdd, *this);
1951
1952                     list<SubTree>::const_iterator i;
1953                     for(i=adds.begin(); i!=adds.end(); ++i)
1954                         tmp.AddParam(*i);
1955
1956                     ReplaceWith(tmp);
1957                 }
1958             }
1959         }
1960         if(GetOp() == cAdd)
1961         {
1962             // Check which ones are logs.
1963             list<pit> poslist;
1964
1965             for(pit a=GetBegin(); a!=GetEnd(); ++a)
1966             {
1967                 const SubTree &pa = *a;
1968                 if(pa->GetOp() == cLog)
1969                     poslist.push_back(a);
1970             }
1971
1972             if(poslist.size() >= 2)
1973             {
1974                 CodeTree tmp(cMul, 1.0); // eek
1975
1976                 list<pit>::const_iterator j;
1977                 for(j=poslist.begin(); j!=poslist.end(); ++j)
1978                 {
1979                     const SubTree &pb = **j;
1980                     // Take all of its children
1981                     for(pcit b=pb->GetBegin();
1982                              b!=pb->GetEnd();
1983                              ++b)
1984                     {
1985                         SubTree tmp2 = *b;
1986                         if(pb.getsign()) tmp2.Negate();
1987                         tmp.AddParam(tmp2);
1988                     }
1989                     Erase(*j);
1990                 }
1991                 poslist.clear();
1992
1993                 AddParam(CodeTree(cLog, tmp));
1994             }
1995             // Done, hopefully
1996         }
1997         if(GetOp() == cPow)
1998         {
1999             const SubTree &p0 = getp0();
2000             SubTree &p1 = getp1();
2001
2002             if(p0->IsImmed() && p0->GetImmed() == CONSTANT_E
2003             && p1->GetOp() == cLog)
2004             {
2005                 // pow(const_E, log(x)) = x
2006                 ReplaceWith(*(p1->getp0()));
2007             }
2008             else if(p1->GetOp() == cMul)
2009             {
2010                 //bool didsomething = true;
2011
2012                 pit poslogpos; bool foundposlog = false;
2013                 pit neglogpos; bool foundneglog = false;
2014
2015                 ConstList cl = p1->BuildConstList();
2016
2017                 for(pit a=p1->GetBegin(); a!=p1->GetEnd(); ++a)
2018                 {
2019                     const SubTree &pa = *a;
2020                     if(pa->GetOp() == cLog)
2021                     {
2022                         if(!pa.getsign())
2023                         {
2024                             foundposlog = true;
2025                             poslogpos   = a;
2026                         }
2027                         else if(*p0 == *(pa->getp0()))
2028                         {
2029                             foundneglog = true;
2030                             neglogpos   = a;
2031                         }
2032                     }
2033                 }
2034
2035                 if(p0->IsImmed()
2036                 && p0->GetImmed() == 10.0
2037                 && cl.value == CONSTANT_L10I
2038                 && foundposlog)
2039                 {
2040                     SubTree base = (*poslogpos)->getp0();
2041                     p1->KillConst(cl);
2042                     p1->Erase(poslogpos);
2043                     p1->OptimizeRedundant();
2044                     SubTree mul = p1;
2045
2046                     ReplaceWith(cPow, base, mul);
2047
2048                     // FIXME: what optimizations should be done now?
2049                     return;
2050                 }
2051
2052                 // Put back the constant
2053                 FinishConst(cl);
2054
2055                 if(p0->IsImmed()
2056                 && p0->GetImmed() == CONSTANT_E
2057                 && foundposlog)
2058                 {
2059                     SubTree base = (*poslogpos)->getp0();
2060                     p1->Erase(poslogpos);
2061
2062                     p1->OptimizeRedundant();
2063                     SubTree mul = p1;
2064
2065                     ReplaceWith(cPow, base, mul);
2066
2067                     // FIXME: what optimizations should be done now?
2068                     return;
2069                 }
2070
2071                 if(foundposlog
2072                 && foundneglog
2073                 && *((*neglogpos)->getp0()) == *p0)
2074                 {
2075                     SubTree base = (*poslogpos)->getp0();
2076                     p1->Erase(poslogpos);
2077                     p1->Erase(neglogpos);
2078
2079                     p1->OptimizeRedundant();
2080                     SubTree mul = p1;
2081
2082                     ReplaceWith(cPow, base, mul);
2083
2084                     // FIXME: what optimizations should be done now?
2085                     return;
2086                 }
2087             }
2088         }
2089     }
2090
2091     void OptimizeFunctionCalls()
2092     {
2093         /* Goals: sin(asin(x)) = x
2094          *        cos(acos(x)) = x
2095          *        tan(atan(x)) = x
2096          * NOTE:
2097          *   Do NOT do these:
2098          *     asin(sin(x))
2099          *     acos(cos(x))
2100          *     atan(tan(x))
2101          *   Because someone might want to wrap the angle.
2102          */
2103         // FIXME: TODO
2104     }
2105
2106     void OptimizePowMulAdd()
2107     {
2108         // x^3 * x -> x^4
2109         // x*3 + x -> x*4
2110         // FIXME: Do those
2111
2112         // x^1 -> x
2113         if(GetOp() == cPow)
2114         {
2115             const SubTree &base     = getp0();
2116             const SubTree &exponent = getp1();
2117
2118             if(exponent->IsImmed())
2119             {
2120                 if(exponent->GetImmed() == 1.0)
2121                     ReplaceWith(*base);
2122                 else if(exponent->GetImmed() == 0.0
2123                      && base->NonZero())
2124                     ReplaceWithConst(1.0);
2125             }
2126         }
2127     }
2128
2129     void OptimizeExponents()
2130     {
2131         /* Goals:
2132          *     (x^y)^z   -> x^(y*z)
2133          *     x^y * x^z -> x^(y+z)
2134          */
2135         // First move to exponential form.
2136         OptimizeLinearCombine();
2137
2138         bool didchanges = false;
2139
2140     Redo:
2141         if(GetOp() == cPow)
2142         {
2143             // (x^y)^z   -> x^(y*z)
2144
2145             const SubTree &p0 = getp0();
2146             const SubTree &p1 = getp1();
2147             if(p0->GetOp() == cPow)
2148             {
2149                 CodeTree tmp(cMul, p0->getp1(), p1);
2150                 tmp.Optimize();
2151
2152                 ReplaceWith(cPow, p0->getp0(), tmp);
2153
2154                 didchanges = true;
2155                 goto Redo;
2156             }
2157         }
2158         if(GetOp() == cMul)
2159         {
2160             // x^y * x^z -> x^(y+z)
2161
2162             for(pit a=GetBegin(); a!=GetEnd(); ++a)
2163             {
2164                 const SubTree &pa = *a;
2165
2166                 if(pa->GetOp() != cPow) continue;
2167
2168                 list<pit> poslist;
2169
2170                 for(pit b=a; ++b != GetEnd(); )
2171                 {
2172                     const SubTree &pb = *b;
2173                     if(pb->GetOp() == cPow
2174                     && *(pa->getp0())
2175                     == *(pb->getp0()))
2176                     {
2177                         poslist.push_back(b);
2178                     }
2179                 }
2180
2181                 if(poslist.size() >= 1)
2182                 {
2183                     poslist.push_back(a);
2184
2185                     CodeTree base = *(pa->getp0());
2186
2187                     CodeTree exponent(cAdd, 0.0); //eek
2188
2189                     // Collect all exponents to cAdd
2190                     list<pit>::const_iterator i;
2191                     for(i=poslist.begin(); i!=poslist.end(); ++i)
2192                     {
2193                         const SubTree &pb = **i;
2194
2195                         SubTree tmp2 = pb->getp1();
2196                         if(pb.getsign()) tmp2.Invert();
2197
2198                         exponent.AddParam(tmp2);
2199                     }
2200
2201                     exponent.Optimize();
2202
2203                     CodeTree result(cPow, base, exponent);
2204
2205                     for(i=poslist.begin(); i!=poslist.end(); ++i)
2206                         Erase(*i);
2207                     poslist.clear();
2208
2209                     AddParam(result); // We're cMul, remember
2210
2211                     didchanges = true;
2212                     goto Redo;
2213                 }
2214             }
2215         }
2216
2217         OptimizePowMulAdd();
2218
2219         if(didchanges)
2220         {
2221             // As a result, there might be need for this:
2222             OptimizeConflict();
2223         }
2224     }
2225
2226     void OptimizeLinearExplode()
2227     {
2228         // x^2 -> x*x
2229         // But only if x is just a simple thing
2230
2231         // Won't work on anything else.
2232         if(GetOp() != cPow) return;
2233
2234         // TODO TODO TODO
2235     }
2236
2237     void OptimizePascal()
2238     {
2239 #if 0    // Too big, too specific, etc
2240
2241         // Won't work on anything else.
2242         if(GetOp() != cAdd) return;
2243
2244         // Must be done after OptimizeLinearCombine();
2245
2246         // Don't need pascal triangle
2247         // Coefficient for x^a * y^b * z^c = 3! / (a! * b! * c!)
2248
2249         // We are greedy and want other than just binomials
2250         // FIXME
2251
2252         // note: partial ones are also nice
2253         //     x*x + x*y + y*y
2254         //   = (x+y)^2 - x*y
2255         //
2256         //     x x * x y * + y y * +
2257         // ->  x y + dup * x y * -
2258 #endif
2259     }
2260
2261 public:
2262
2263     void Optimize();
2264
2265     void Assemble(vector<unsigned> &byteCode,
2266                   vector<double>   &immed) const;
2267
2268     void FinalOptimize()
2269     {
2270         // First optimize each parameter.
2271         for(pit a=GetBegin(); a!=GetEnd(); ++a)
2272             (*a)->FinalOptimize();
2273
2274         /* These things are to be done:
2275          *
2276          * x * CONSTANT_DR        -> cDeg(x)
2277          * x * CONSTANT_RD        -> cRad(x)
2278          * pow(x, 0.5)            -> sqrt(x)
2279          * log(x) * CONSTANT_L10I -> log10(x)
2280          * pow(CONSTANT_E, x)     -> exp(x)
2281          * inv(sin(x))            -> csc(x)
2282          * inv(cos(x))            -> sec(x)
2283          * inv(tan(x))            -> cot(x)
2284          */
2285
2286
2287         if(GetOp() == cPow)
2288         {
2289             const SubTree &p0 = getp0();
2290             const SubTree &p1 = getp1();
2291             if(p0->GetOp()    == cImmed
2292             && p0->GetImmed() == CONSTANT_E)
2293             {
2294                 ReplaceWith(cExp, p1);
2295             }
2296             else if(p1->GetOp()    == cImmed
2297                  && p1->GetImmed() == 0.5)
2298             {
2299                 ReplaceWith(cSqrt, p0);
2300             }
2301         }
2302         if(GetOp() == cMul)
2303         {
2304             if(GetArgCount() == 1 && getp0().getsign())
2305             {
2306                 /***/if(getp0()->GetOp() == cSin)ReplaceWith(cCsc, getp0()->getp0());
2307                 else if(getp0()->GetOp() == cCos)ReplaceWith(cSec, getp0()->getp0());
2308                 else if(getp0()->GetOp() == cTan)ReplaceWith(cCot, getp0()->getp0());
2309             }
2310         }
2311         // Separate "if", because op may have just changed
2312         if(GetOp() == cMul)
2313         {
2314             CodeTree *found_log = 0;
2315
2316             ConstList cl = BuildConstList();
2317
2318             for(pit a=GetBegin(); a!=GetEnd(); ++a)
2319             {
2320                 SubTree &pa = *a;
2321                 if(pa->GetOp() == cLog && !pa.getsign())
2322                     found_log = &*pa;
2323             }
2324             if(cl.value == CONSTANT_L10I && found_log)
2325             {
2326                 // Change the log() to log10()
2327                 found_log->SetOp(cLog10);
2328                 // And forget the constant
2329                 KillConst(cl);
2330             }
2331             else if(cl.value == CONSTANT_DR)
2332             {
2333                 OptimizeRedundant();
2334                 ReplaceWith(cDeg, *this);
2335             }
2336             else if(cl.value == CONSTANT_RD)
2337             {
2338                 OptimizeRedundant();
2339                 ReplaceWith(cRad, *this);
2340             }
2341             else FinishConst(cl);
2342         }
2343
2344         SortIfPossible();
2345     }
2346 };
2347
2348 void CodeTreeDataPtr::Shock()
2349 {
2350  /*
2351     PrepareForWrite();
2352     paramlist &p2 = (*this)->args;
2353     for(paramlist::iterator i=p2.begin(); i!=p2.end(); ++i)
2354     {
2355         (*i)->data.Shock();
2356     }
2357  */
2358 }
2359
2360 CodeTree::ConstList CodeTree::BuildConstList()
2361 {
2362     ConstList result;
2363     result.value     =
2364     result.voidvalue = GetOp()==cMul ? 1.0 : 0.0;
2365
2366     list<pit> &cp = result.cp;
2367     for(pit b, a=GetBegin(); a!=GetEnd(); a=b)
2368     {
2369         SubTree &pa = *a;  b=a; ++b;
2370         if(!pa->IsImmed()) continue;
2371
2372         double thisvalue = pa->GetImmed();
2373         if(thisvalue == result.voidvalue)
2374         {
2375             // This value is no good, forget it
2376             Erase(a);
2377             continue;
2378         }
2379         if(GetOp() == cMul)
2380             result.value *= thisvalue;
2381         else
2382             result.value += thisvalue;
2383         cp.push_back(a);
2384     }
2385     if(GetOp() == cMul)
2386     {
2387         /*
2388           Jos joku niistä arvoista on -1 eikä se ole ainoa arvo,
2389           niin joku muu niistä arvoista negatoidaan.
2390         */
2391         for(bool done=false; cp.size() > 1 && !done; )
2392         {
2393             done = true;
2394             for(list<pit>::iterator b,a=cp.begin(); a!=cp.end(); a=b)
2395             {
2396                 b=a; ++b;
2397                 if((**a)->GetImmed() == -1.0)
2398                 {
2399                     Erase(*a);
2400                     cp.erase(a);
2401
2402                     // take randomly something
2403                     (**cp.begin())->data->NegateImmed();
2404                     if(cp.size() < 2)break;
2405                     done = false;
2406                 }
2407             }
2408         }
2409     }
2410     return result;
2411 }
2412
2413 void CodeTree::Assemble
2414    (vector<unsigned> &byteCode,
2415     vector<double>   &immed) const
2416 {
2417     #define AddCmd(op) byteCode.push_back((op))
2418     #define AddConst(v) do { \
2419         byteCode.push_back(cImmed); \
2420         immed.push_back((v)); \
2421     } while(0)
2422
2423     if(IsVar())
2424     {
2425         AddCmd(GetVar());
2426         return;
2427     }
2428     if(IsImmed())
2429     {
2430         AddConst(GetImmed());
2431         return;
2432     }
2433
2434     switch(GetOp())
2435     {
2436         case cAdd:
2437         case cMul:
2438         {
2439             unsigned opcount = 0;
2440             for(pcit a=GetBegin(); a!=GetEnd(); ++a)
2441             {
2442                 const SubTree &pa = *a;
2443
2444                 if(opcount < 2) ++opcount;
2445
2446                 bool pnega = pa.getsign();
2447
2448                 bool done = false;
2449                 if(pa->IsImmed())
2450                 {
2451                     if(GetOp() == cMul
2452                     && pa->data->IsInverted()
2453                     && (pnega || opcount==2)
2454                       )
2455                     {
2456                         CodeTree tmp = *pa;
2457                         tmp.data->InvertImmed();
2458                         tmp.Assemble(byteCode, immed);
2459                         pnega = !pnega;
2460                         done = true;
2461                     }
2462                     else if(GetOp() == cAdd
2463                     && (pa->data->IsNegatedOriginal()
2464                 //     || pa->GetImmed() < 0
2465                        )
2466                     && (pnega || opcount==2)
2467                            )
2468                     {
2469                         CodeTree tmp = *pa;
2470                         tmp.data->NegateImmed();
2471                         tmp.Assemble(byteCode, immed);
2472                         pnega = !pnega;
2473                         done = true;
2474                     }
2475                 }
2476                 if(!done)
2477                     pa->Assemble(byteCode, immed);
2478
2479                 if(opcount == 2)
2480                 {
2481                     unsigned tmpop = GetOp();
2482                     if(pnega) // negate
2483                     {
2484                         tmpop = (tmpop == cMul) ? cDiv : cSub;
2485                     }
2486                     AddCmd(tmpop);
2487                 }
2488                 else if(pnega)
2489                 {
2490                     if(GetOp() == cMul) AddCmd(cInv);
2491                     else AddCmd(cNeg);
2492                 }
2493             }
2494             break;
2495         }
2496         case cIf:
2497         {
2498             // If the parameter amount is != 3, we're screwed.
2499             getp0()->Assemble(byteCode, immed);
2500
2501             unsigned ofs = byteCode.size();
2502             AddCmd(cIf);
2503             AddCmd(0); // code index
2504             AddCmd(0); // immed index
2505
2506             getp1()->Assemble(byteCode, immed);
2507
2508             byteCode[ofs+1] = byteCode.size()+2;
2509             byteCode[ofs+2] = immed.size();
2510
2511             ofs = byteCode.size();
2512             AddCmd(cJump);
2513             AddCmd(0); // code index
2514             AddCmd(0); // immed index
2515
2516             getp2()->Assemble(byteCode, immed);
2517
2518             byteCode[ofs+1] = byteCode.size()-1;
2519             byteCode[ofs+2] = immed.size();
2520
2521             break;
2522         }
2523         case cFCall:
2524         {
2525             // If the parameter count is invalid, we're screwed.
2526             for(pcit a=GetBegin(); a!=GetEnd(); ++a)
2527             {
2528                 const SubTree &pa = *a;
2529                 pa->Assemble(byteCode, immed);
2530             }
2531             AddCmd(GetOp());
2532             AddCmd(data->GetFuncNo());
2533             break;
2534         }
2535         case cPCall:
2536         {
2537             // If the parameter count is invalid, we're screwed.
2538             for(pcit a=GetBegin(); a!=GetEnd(); ++a)
2539             {
2540                 const SubTree &pa = *a;
2541                 pa->Assemble(byteCode, immed);
2542             }
2543             AddCmd(GetOp());
2544             AddCmd(data->GetFuncNo());
2545             break;
2546         }
2547         default:
2548         {
2549             // If the parameter count is invalid, we're screwed.
2550             for(pcit a=GetBegin(); a!=GetEnd(); ++a)
2551             {
2552                 const SubTree &pa = *a;
2553                 pa->Assemble(byteCode, immed);
2554             }
2555             AddCmd(GetOp());
2556             break;
2557         }
2558     }
2559 }
2560
2561 void CodeTree::Optimize()
2562 {
2563     // Phase:
2564     //   Phase 0: Do local optimizations.
2565     //   Phase 1: Optimize each.
2566     //   Phase 2: Do local optimizations again.
2567
2568     for(unsigned phase=0; phase<=2; ++phase)
2569     {
2570         if(phase == 1)
2571         {
2572             // Optimize each parameter.
2573             for(pit a=GetBegin(); a!=GetEnd(); ++a)
2574             {
2575                 (*a)->Optimize();
2576                 CHECKCONSTNEG(*a, GetOp());
2577             }
2578             continue;
2579         }
2580         if(phase == 0 || phase == 2)
2581         {
2582             // Do local optimizations.
2583
2584             OptimizeConstantMath1();
2585             OptimizeLogarithm();
2586             OptimizeFunctionCalls();
2587             OptimizeExponents();
2588             OptimizeLinearExplode();
2589             OptimizePascal();
2590
2591             /* Optimization paths:
2592
2593                doublenegations=
2594                redundant= * doublenegations
2595                conflict= * redundant
2596                addmulflat=
2597                constantmath1= addmulflat * conflict
2598                linearcombine= conflict * addmulflat¹ redundant¹
2599                powmuladd=
2600                exponents= linearcombine * powmuladd conflict¹
2601                logarithm= exponents *
2602                functioncalls= IDLE
2603                linearexplode= IDLE
2604                pascal= IDLE
2605
2606                * = actions here
2607                Â¹ = only if made changes
2608             */
2609         }
2610     }
2611 }
2612
2613
2614 bool CodeTree::operator== (const CodeTree& b) const
2615 {
2616     if(GetOp() != b.GetOp()) return false;
2617     if(IsImmed()) if(GetImmed()  != b.GetImmed())  return false;
2618     if(IsVar())   if(GetVar()    != b.GetVar())    return false;
2619     if(data->IsFunc())
2620         if(data->GetFuncNo() != b.data->GetFuncNo()) return false;
2621     return data->args == b.data->args;
2622 }
2623
2624 bool CodeTree::operator< (const CodeTree& b) const
2625 {
2626     if(GetArgCount() != b.GetArgCount())
2627         return GetArgCount() > b.GetArgCount();
2628
2629     if(GetOp() != b.GetOp())
2630     {
2631         // sort immeds last
2632         if(IsImmed() != b.IsImmed()) return IsImmed() < b.IsImmed();
2633
2634         return GetOp() < b.GetOp();
2635     }
2636
2637     if(IsImmed())
2638     {
2639         if(GetImmed() != b.GetImmed()) return GetImmed() < b.GetImmed();
2640     }
2641     if(IsVar() && GetVar() != b.GetVar())
2642     {
2643         return GetVar() < b.GetVar();
2644     }
2645     if(data->IsFunc() && data->GetFuncNo() != b.data->GetFuncNo())
2646     {
2647         return data->GetFuncNo() < b.data->GetFuncNo();
2648     }
2649
2650     pcit i = GetBegin(), j = b.GetBegin();
2651     for(; i != GetEnd(); ++i, ++j)
2652     {
2653         const SubTree &pa = *i, &pb = *j;
2654
2655         if(!(pa == pb))
2656             return pa < pb;
2657     }
2658     return false;
2659 }
2660
2661 namespace {
2662 bool IsNegate(const SubTree &p1, const SubTree &p2) /*const */
2663 {
2664     if(p1->IsImmed() && p2->IsImmed())
2665     {
2666         return p1->GetImmed() == -p2->GetImmed();
2667     }
2668     if(p1.getsign() == p2.getsign()) return false;
2669     return *p1 == *p2;
2670 }
2671 bool IsInverse(const SubTree &p1, const SubTree &p2) /*const*/
2672 {
2673     if(p1->IsImmed() && p2->IsImmed())
2674     {
2675         // FIXME: potential divide by zero.
2676         return p1->GetImmed() == 1.0 / p2->GetImmed();
2677     }
2678     if(p1.getsign() == p2.getsign()) return false;
2679     return *p1 == *p2;
2680 }
2681 }
2682
2683 SubTree::SubTree() : tree(new CodeTree), sign(false)
2684 {
2685 }
2686
2687 SubTree::SubTree(const SubTree &b) : tree(new CodeTree(*b.tree)), sign(b.sign)
2688 {
2689 }
2690
2691 #define SubTreeDecl(p1, p2) \
2692     SubTree::SubTree p1 : tree(new CodeTree p2), sign(false) { }
2693
2694 SubTreeDecl( (const struct CodeTree &b), (b) )
2695 SubTreeDecl( (double value),             (value) )
2696
2697 #undef SubTreeDecl
2698
2699 SubTree::~SubTree()
2700 {
2701     delete tree; tree=0;
2702 }
2703
2704 const SubTree &SubTree::operator= (const SubTree &b)
2705 {
2706     sign = b.sign;
2707     CodeTree *oldtree = tree;
2708     tree = new CodeTree(*b.tree);
2709     delete oldtree;
2710     return *this;
2711 }
2712 const SubTree &SubTree::operator= (const CodeTree &b)
2713 {
2714     sign = false;
2715     CodeTree *oldtree = tree;
2716     tree = new CodeTree(b);
2717     delete oldtree;
2718     return *this;
2719 }
2720
2721 bool SubTree::operator< (const SubTree& b) const
2722 {
2723     if(getsign() != b.getsign()) return getsign() < b.getsign();
2724     return *tree < *b.tree;
2725 }
2726 bool SubTree::operator== (const SubTree& b) const
2727 {
2728     return sign == b.sign && *tree == *b.tree;
2729 }
2730 void SubTree::Negate() // Note: Parent must be cAdd
2731 {
2732     flipsign();
2733     CheckConstNeg();
2734 }
2735 void SubTree::CheckConstNeg()
2736 {
2737     if(tree->IsImmed() && getsign())
2738     {
2739         tree->NegateImmed();
2740         sign = false;
2741     }
2742 }
2743 void SubTree::Invert() // Note: Parent must be cMul
2744 {
2745     flipsign();
2746     CheckConstInv();
2747 }
2748 void SubTree::CheckConstInv()
2749 {
2750     if(tree->IsImmed() && getsign())
2751     {
2752         tree->InvertImmed();
2753         sign = false;
2754     }
2755 }
2756
2757 void FunctionParser::MakeTree(struct CodeTree *result) const
2758 {
2759     vector<CodeTree> stack(1);
2760
2761     #define GROW(n) do { \
2762         stacktop += n; \
2763         if(stack.size() <= stacktop) stack.resize(stacktop+1); \
2764     } while(0)
2765
2766     #define EAT(n, opcode) do { \
2767         unsigned newstacktop = stacktop-n; \
2768         stack[stacktop].SetOp((opcode)); \
2769         for(unsigned a=0, b=(n); a<b; ++a) \
2770             stack[stacktop].AddParam(stack[newstacktop+a]); \
2771         stack.erase(stack.begin() + newstacktop, \
2772                     stack.begin() + stacktop); \
2773         stacktop = newstacktop; GROW(1); \
2774     } while(0)
2775
2776     #define ADDCONST(n) do { \
2777         stack[stacktop].SetImmed((n)); \
2778         GROW(1); \
2779     } while(0)
2780
2781     unsigned stacktop=0;
2782
2783     list<unsigned> labels;
2784
2785     for(unsigned IP=0, DP=0; ; ++IP)
2786     {
2787         while(labels.size() > 0
2788         && *labels.begin() == IP)
2789         {
2790             // The "else" of an "if" ends here
2791             EAT(3, cIf);
2792             labels.erase(labels.begin());
2793         }
2794
2795         if(IP >= Comp.ByteCodeSize)
2796         {
2797             break;
2798         }
2799
2800         unsigned opcode = Comp.ByteCode[IP];
2801
2802         if(opcode == cIf)
2803         {
2804             IP += 2;
2805         }
2806         else if(opcode == cJump)
2807         {
2808             labels.push_front(Comp.ByteCode[IP+1]+1);
2809             IP += 2;
2810         }
2811         else if(opcode == cImmed)
2812         {
2813             ADDCONST(Comp.Immed[DP++]);
2814         }
2815         else if(opcode < VarBegin)
2816         {
2817             switch(opcode)
2818             {
2819                 // Unary operators
2820                 case cNeg:
2821                 {
2822                     EAT(1, cAdd); // Unary minus is negative adding.
2823                     stack[stacktop-1].getp0().Negate();
2824                     break;
2825                 }
2826                 // Binary operators
2827                 case cSub:
2828                 {
2829                     EAT(2, cAdd); // Minus is negative adding
2830                     stack[stacktop-1].getp1().Negate();
2831                     break;
2832                 }
2833                 case cDiv:
2834                 {
2835                     EAT(2, cMul); // Divide is inverse multiply
2836                     stack[stacktop-1].getp1().Invert();
2837                     break;
2838                 }
2839
2840                 // ADD ALL TWO PARAMETER NON-FUNCTIONS HERE
2841                 case cAdd: case cMul:
2842                 case cMod: case cPow:
2843                 case cEqual: case cLess: case cGreater:
2844                 case cAnd: case cOr:
2845                     EAT(2, opcode);
2846                     break;
2847
2848                 case cFCall:
2849                 {
2850                     unsigned index = Comp.ByteCode[++IP];
2851                     unsigned params = FuncPtrs[index].params;
2852                     EAT(params, opcode);
2853                     stack[stacktop-1].data->SetFuncNo(index);
2854                     break;
2855                 }
2856                 case cPCall:
2857                 {
2858                     unsigned index = Comp.ByteCode[++IP];
2859                     unsigned params = FuncParsers[index]->varAmount;
2860                     EAT(params, opcode);
2861                     stack[stacktop-1].data->SetFuncNo(index);
2862                     break;
2863                 }
2864
2865                 // Converted to cMul on fly
2866                 case cDeg:
2867                     ADDCONST(CONSTANT_DR);
2868                     EAT(2, cMul);
2869                     break;
2870
2871                 // Converted to cMul on fly
2872                 case cRad:
2873                     ADDCONST(CONSTANT_RD);
2874                     EAT(2, cMul);
2875                     break;
2876
2877                 // Functions
2878                 default:
2879                 {
2880                     const FuncDefinition& func = Functions[opcode-cAbs];
2881
2882                     unsigned paramcount = func.params;
2883 #ifndef DISABLE_EVAL
2884                     if(opcode == cEval) paramcount = varAmount;
2885 #endif
2886                     if(opcode == cSqrt)
2887                     {
2888                         // Converted on fly: sqrt(x) = x^0.5
2889                         opcode = cPow;
2890                         paramcount = 2;
2891                         ADDCONST(0.5);
2892                     }
2893                     if(opcode == cExp)
2894                     {
2895                         // Converted on fly: exp(x) = CONSTANT_E^x
2896
2897                         opcode = cPow;
2898                         paramcount = 2;
2899                         // reverse the parameters... kludgey
2900                         stack[stacktop] = stack[stacktop-1];
2901                         stack[stacktop-1].SetImmed(CONSTANT_E);
2902                         GROW(1);
2903                     }
2904                     bool do_inv = false;
2905                     if(opcode == cCot) { do_inv = true; opcode = cTan; }
2906                     if(opcode == cCsc) { do_inv = true; opcode = cSin; }
2907                     if(opcode == cSec) { do_inv = true; opcode = cCos; }
2908
2909                     bool do_log10 = false;
2910                     if(opcode == cLog10)
2911                     {
2912                         // Converted on fly: log10(x) = log(x) * CONSTANT_L10I
2913                         opcode = cLog;
2914                         do_log10 = true;
2915                     }
2916                     EAT(paramcount, opcode);
2917                     if(do_log10)
2918                     {
2919                         ADDCONST(CONSTANT_L10I);
2920                         EAT(2, cMul);
2921                     }
2922                     if(do_inv)
2923                     {
2924                         // Unary cMul, inverted. No need for "1.0"
2925                         EAT(1, cMul);
2926                         stack[stacktop-1].getp0().Invert();
2927                     }
2928                     break;
2929                 }
2930             }
2931         }
2932         else
2933         {
2934             stack[stacktop].SetVar(opcode);
2935             GROW(1);
2936         }
2937     }
2938
2939     if(!stacktop)
2940     {
2941         // ERROR: Stack does not have any values!
2942         return;
2943     }
2944
2945     --stacktop; // Ignore the last element, it is always nop (cAdd).
2946
2947     if(stacktop > 0)
2948     {
2949         // ERROR: Stack has too many values!
2950         return;
2951     }
2952
2953     // Okay, the tree is now stack[0]
2954     *result = stack[0];
2955 }
2956
2957 void FunctionParser::Optimize()
2958 {
2959
2960     CodeTree tree;
2961     MakeTree(&tree);
2962
2963     // Do all sorts of optimizations
2964     tree.Optimize();
2965     // Last changes before assembly
2966     tree.FinalOptimize();
2967
2968     // Now rebuild from the tree.
2969
2970     vector<unsigned> byteCode;
2971     vector<double> immed;
2972
2973 #if 0
2974     byteCode.resize(Comp.ByteCodeSize);
2975     for(unsigned a=0; a<Comp.ByteCodeSize; ++a)byteCode[a] = Comp.ByteCode[a];
2976
2977     immed.resize(Comp.ImmedSize);
2978     for(unsigned a=0; a<Comp.ImmedSize; ++a)immed[a] = Comp.Immed[a];
2979 #else
2980     byteCode.clear(); immed.clear();
2981     tree.Assemble(byteCode, immed);
2982 #endif
2983
2984     delete[] Comp.ByteCode; Comp.ByteCode = 0;
2985     if((Comp.ByteCodeSize = byteCode.size()) > 0)
2986     {
2987         Comp.ByteCode = new unsigned[Comp.ByteCodeSize];
2988         for(unsigned a=0; a<byteCode.size(); ++a)
2989             Comp.ByteCode[a] = byteCode[a];
2990     }
2991
2992     delete[] Comp.Immed; Comp.Immed = 0;
2993     if((Comp.ImmedSize = immed.size()) > 0)
2994     {
2995         Comp.Immed = new double[Comp.ImmedSize];
2996         for(unsigned a=0; a<immed.size(); ++a)
2997             Comp.Immed[a] = immed[a];
2998     }
2999 }
3000
3001
3002 #else /* !SUPPORT_OPTIMIZER */
3003
3004 /* keep the linker happy */
3005 void FunctionParser::MakeTree(struct CodeTree *) const {}
3006 void FunctionParser::Optimize()
3007 {
3008     // Do nothing if no optimizations are supported.
3009 }
3010 #endif