Scribus Developer Documentation - API
logo-pen

Scribus Developer Documentation - API

fpoptimizer.cpp

Go to the documentation of this file.
00001 /*
00002 For general Scribus (>=1.3.2) copyright and licensing information please refer
00003 to the COPYING file provided with the program. Following this notice may exist
00004 a copyright and/or license notice that predates the release of Scribus 1.3.2
00005 for which a new license (GPL+exception) is in place.
00006 */
00007 //===========================================
00008 // Function parser v2.8 optimizer by Bisqwit
00009 //===========================================
00010 
00011 /*
00012  NOTE!
00013  ----
00014  Everything that goes into the #ifndef SUPPORT_OPTIMIZER part
00015  (ie. when SUPPORT_OPTIMIZER is not defined) should be put in
00016  the end of fparser.cc file, not in this file.
00017 
00018  Everything in this file should be inside the #ifdef SUPPORT_OPTIMIZER
00019  block (except the #include "fpconfig.hh" line).
00020 */
00021 
00022 
00023 #include "fpconfig.h"
00024 
00025 #ifdef SUPPORT_OPTIMIZER
00026 
00027 #include "fparser.h"
00028 #include "fptypes.h"
00029 using namespace FUNCTIONPARSERTYPES;
00030 
00031 #include <cmath>
00032 #include <list>
00033 #include <utility>
00034 #include <cassert>
00035 using namespace std;
00036 
00037 #ifndef M_PI
00038 #define M_PI 3.1415926535897932384626433832795
00039 #endif
00040 
00041 #define CONSTANT_E     2.71828182845904509080  // exp(1)
00042 #define CONSTANT_PI    M_PI                    // atan2(0,-1)
00043 #define CONSTANT_L10   2.30258509299404590109  // log(10)
00044 #define CONSTANT_L10I  0.43429448190325176116  // 1/log(10)
00045 #define CONSTANT_L10E  CONSTANT_L10I           // log10(e)
00046 #define CONSTANT_L10EI CONSTANT_L10            // 1/log10(e)
00047 #define CONSTANT_DR    (180.0 / M_PI)          // 180/pi
00048 #define CONSTANT_RD    (M_PI / 180.0)          // pi/180
00049 
00050 namespace {
00051 inline double Min(double d1, double d2)
00052 {
00053     return d1<d2 ? d1 : d2;
00054 }
00055 inline double Max(double d1, double d2)
00056 {
00057     return d1>d2 ? d1 : d2;
00058 }
00059 
00060 class compres
00061 {
00062     // states: 0=false, 1=true, 2=unknown
00063 public:
00064     compres(bool b) : state(b) {}
00065     compres(char v) : state(v) {}
00066     // is it?
00067     operator bool() const { return state != 0; }
00068     // is it not?
00069     bool operator! () const { return state != 1; }
00070     bool operator==(bool b) const { return state != !b; }
00071     bool operator!=(bool b) const { return state != b; }
00072 private:
00073     char state;
00074 };
00075 
00076 const compres maybe = (char)2;
00077 
00078 struct CodeTree;
00079 
00080 class SubTree
00081 {
00082     CodeTree *tree;
00083     bool sign;  // Only possible when parent is cAdd or cMul
00084 
00085     inline void flipsign() { sign = !sign; }
00086 public:
00087     SubTree();
00088     SubTree(double value);
00089     SubTree(const SubTree &b);
00090     SubTree(const CodeTree &b);
00091 
00092     ~SubTree();
00093     const SubTree &operator= (const SubTree &b);
00094     const SubTree &operator= (const CodeTree &b);
00095 
00096     bool getsign() const { return sign; }
00097 
00098     const CodeTree* operator-> () const { return tree; }
00099     const CodeTree& operator* () const { return *tree; }
00100     struct CodeTree* operator-> () { return tree; }
00101     struct CodeTree& operator* () { return *tree; }
00102 
00103     bool operator< (const SubTree& b) const;
00104     bool operator== (const SubTree& b) const;
00105     void Negate(); // Note: Parent must be cAdd
00106     void Invert(); // Note: Parent must be cMul
00107 
00108     void CheckConstNeg();
00109     void CheckConstInv();
00110 };
00111 
00112 bool IsNegate(const SubTree &p1, const SubTree &p2);
00113 bool IsInverse(const SubTree &p1, const SubTree &p2);
00114 
00115 typedef list<SubTree> paramlist;
00116 
00117 struct CodeTreeData
00118 {
00119     paramlist args;
00120 
00121 private:
00122     unsigned op;       // Operation
00123     double value;      // In case of cImmed
00124     unsigned var;      // In case of cVar
00125     unsigned funcno;   // In case of cFCall, cPCall
00126 
00127 public:
00128     CodeTreeData() : op(cAdd) {}
00129     ~CodeTreeData() {}
00130 
00131     void SetOp(unsigned newop)     { op=newop; }
00132     void SetFuncNo(unsigned newno) { funcno=newno; }
00133     unsigned GetFuncNo() const { return funcno; }
00134 
00135     bool IsFunc() const  { return op == cFCall || op == cPCall; }
00136     bool IsImmed() const { return op == cImmed; }
00137     bool IsVar() const   { return op == cVar; }
00138     inline unsigned GetOp() const { return op; }
00139     inline double GetImmed() const
00140     {
00141         return value;
00142     }
00143     inline unsigned GetVar() const
00144     {
00145         return var;
00146     }
00147 
00148     void AddParam(const SubTree &p)
00149     {
00150         args.push_back(p);
00151     }
00152     void SetVar(unsigned v)
00153     {
00154         args.clear();
00155         op  = cVar;
00156         var = v;
00157     }
00158     void SetImmed(double v)
00159     {
00160         args.clear();
00161         op       = cImmed;
00162         value    = orig = v;
00163         inverted = negated = false;
00164     }
00165     void NegateImmed()
00166     {
00167         negated = !negated;
00168         UpdateValue();
00169     }
00170     void InvertImmed()
00171     {
00172         inverted = !inverted;
00173         UpdateValue();
00174     }
00175 
00176     bool IsOriginal() const { return !(IsInverted() || IsNegated()); }
00177     bool IsInverted() const { return inverted; }
00178     bool IsNegated() const { return negated; }
00179     bool IsInvertedOriginal() const { return IsInverted() && !IsNegated(); }
00180     bool IsNegatedOriginal() const { return !IsInverted() && IsNegated(); }
00181 
00182 private:
00183     void UpdateValue()
00184     {
00185         value = orig;
00186         if(IsInverted()) { value = 1.0 / value;
00187                            // FIXME: potential divide by zero.
00188                          }
00189         if(IsNegated()) value = -value;
00190     }
00191 
00192     double orig;
00193     bool inverted;
00194     bool negated;
00195 protected:
00196     // Ensure we don't accidentally copy this
00197     void operator=(const CodeTreeData &b);
00198 };
00199 
00200 
00201 class CodeTreeDataPtr
00202 {
00203     typedef pair<CodeTreeData, unsigned> p_t;
00204     typedef p_t* pp;
00205     mutable pp p;
00206 
00207     void Alloc()   const { ++p->second; }
00208     void Dealloc() const { if(!--p->second) delete p; p = 0; }
00209 
00210     void PrepareForWrite()
00211     {
00212         // We're ready if we're the only owner.
00213         if(p->second == 1) return;
00214 
00215         // Then make a clone.
00216         p_t *newtree = new p_t(p->first, 1);
00217         // Forget the old
00218         Dealloc();
00219         // Keep the new
00220         p = newtree;
00221     }
00222 
00223 public:
00224     CodeTreeDataPtr() : p(new p_t) { p->second = 1; }
00225     CodeTreeDataPtr(const CodeTreeDataPtr &b): p(b.p) { Alloc(); }
00226     ~CodeTreeDataPtr() { Dealloc(); }
00227     const CodeTreeDataPtr &operator= (const CodeTreeDataPtr &b)
00228     {
00229         b.Alloc();
00230         Dealloc();
00231         p = b.p;
00232         return *this;
00233     }
00234     const CodeTreeData *operator-> () const { return &p->first; }
00235     const CodeTreeData &operator*  () const { return p->first; }
00236     CodeTreeData *operator-> () { PrepareForWrite(); return &p->first; }
00237     CodeTreeData &operator*  () { PrepareForWrite(); return p->first; }
00238 
00239     void Shock();
00240 };
00241 
00242 
00243 #define CHECKCONSTNEG(item, op) \
00244     ((op)==cMul) \
00245        ? (item).CheckConstInv() \
00246        : (item).CheckConstNeg()
00247 
00248 struct CodeTree
00249 {
00250     CodeTreeDataPtr data;
00251 
00252 private:
00253     typedef paramlist::iterator pit;
00254     typedef paramlist::const_iterator pcit;
00255 
00256     /*
00257     template<unsigned v> inline void chk() const
00258     {
00259     }
00260     */
00261 
00262 public:
00263     const pcit GetBegin() const { return data->args.begin(); }
00264     const pcit GetEnd()   const { return data->args.end(); }
00265     const pit GetBegin() { return data->args.begin(); }
00266     const pit GetEnd()   { return data->args.end(); }
00267     const SubTree& getp0() const { /*chk<1>();*/pcit tmp=GetBegin();               return *tmp; }
00268     const SubTree& getp1() const { /*chk<2>();*/pcit tmp=GetBegin(); ++tmp;        return *tmp; }
00269     const SubTree& getp2() const { /*chk<3>();*/pcit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; }
00270     unsigned GetArgCount() const { return data->args.size(); }
00271     void Erase(const pit p)      { data->args.erase(p); }
00272 
00273     SubTree& getp0() { /*chk<1>();*/pit tmp=GetBegin();               return *tmp; }
00274     SubTree& getp1() { /*chk<2>();*/pit tmp=GetBegin(); ++tmp;        return *tmp; }
00275     SubTree& getp2() { /*chk<3>();*/pit tmp=GetBegin(); ++tmp; ++tmp; return *tmp; }
00276 
00277     // set
00278     void SetImmed(double v) { data->SetImmed(v); }
00279     void SetOp(unsigned op) { data->SetOp(op); }
00280     void SetVar(unsigned v) { data->SetVar(v); }
00281     // get
00282     double GetImmed() const { return data->GetImmed(); }
00283     unsigned GetVar() const { return data->GetVar(); }
00284     unsigned GetOp() const  { return data->GetOp(); }
00285     // test
00286     bool IsImmed() const { return data->IsImmed(); }
00287     bool IsVar()   const { return data->IsVar(); }
00288     // act
00289     void AddParam(const SubTree &p) { data->AddParam(p); }
00290     void NegateImmed() { data->NegateImmed(); } // don't use when op!=cImmed
00291     void InvertImmed() { data->InvertImmed(); } // don't use when op!=cImmed
00292 
00293     compres NonZero() const { if(!IsImmed()) return maybe;
00294                               return GetImmed() != 0.0; }
00295     compres IsPositive() const { if(!IsImmed()) return maybe;
00296                                  return GetImmed() > 0.0; }
00297 
00298 private:
00299     struct ConstList
00300     {
00301         double voidvalue;
00302         list<pit> cp;
00303         double value;
00304         unsigned size() const { return cp.size(); }
00305     };
00306     struct ConstList BuildConstList();
00307     void KillConst(const ConstList &cl)
00308     {
00309         for(list<pit>::const_iterator i=cl.cp.begin(); i!=cl.cp.end(); ++i)
00310             Erase(*i);
00311     }
00312     void FinishConst(const ConstList &cl)
00313     {
00314         if(cl.value != cl.voidvalue && cl.size() > 1) AddParam(cl.value);
00315         if(cl.value == cl.voidvalue || cl.size() > 1) KillConst(cl);
00316     }
00317 
00318 public:
00319     CodeTree() {}
00320     CodeTree(double v) { SetImmed(v); }
00321 
00322     CodeTree(unsigned op, const SubTree &p)
00323     {
00324         SetOp(op);
00325         AddParam(p);
00326     }
00327     CodeTree(unsigned op, const SubTree &p1, const SubTree &p2)
00328     {
00329         SetOp(op);
00330         AddParam(p1);
00331         AddParam(p2);
00332     }
00333 
00334     bool operator== (const CodeTree& b) const;
00335     bool operator< (const CodeTree& b) const;
00336 
00337 private:
00338     bool IsSortable() const
00339     {
00340         switch(GetOp())
00341         {
00342             case cAdd:  case cMul:
00343             case cEqual:
00344             case cAnd: case cOr:
00345             case cMax: case cMin:
00346                 return true;
00347             default:
00348                 return false;
00349         }
00350     }
00351     void SortIfPossible()
00352     {
00353         if(IsSortable())
00354         {
00355             data->args.sort();
00356         }
00357     }
00358 
00359     void ReplaceWithConst(double value)
00360     {
00361         SetImmed(value);
00362 
00363         /* REMEMBER TO CALL CheckConstInv / CheckConstNeg
00364          * FOR PARENT SubTree, OR MAYHEM HAPPENS
00365          */
00366     }
00367 
00368     void ReplaceWith(const CodeTree &b)
00369     {
00370         // If b is child of *this, mayhem
00371         // happens. So we first make a clone
00372         // and then proceed with copy.
00373         CodeTreeDataPtr tmp = b.data;
00374         tmp.Shock();
00375         data = tmp;
00376     }
00377 
00378     void ReplaceWith(unsigned op, const SubTree &p)
00379     {
00380         ReplaceWith(CodeTree(op, p));
00381     }
00382 
00383     void ReplaceWith(unsigned op, const SubTree &p1, const SubTree &p2)
00384     {
00385         ReplaceWith(CodeTree(op, p1, p2));
00386     }
00387 
00388     void OptimizeConflict()
00389     {
00390         // This optimization does this: x-x = 0, x/x = 1, a+b-a = b.
00391 
00392         if(GetOp() == cAdd || GetOp() == cMul)
00393         {
00394         Redo:
00395             pit a, b;
00396             for(a=GetBegin(); a!=GetEnd(); ++a)
00397             {
00398                 for(b=GetBegin(); ++b != GetEnd(); )
00399                 {
00400                     const SubTree &p1 = *a;
00401                     const SubTree &p2 = *b;
00402 
00403                     if(GetOp() == cMul ? IsInverse(p1,p2)
00404                                        : IsNegate(p1,p2))
00405                     {
00406                         // These parameters complement each others out
00407                         Erase(b);
00408                         Erase(a);
00409                         goto Redo;
00410                     }
00411                 }
00412             }
00413         }
00414         OptimizeRedundant();
00415     }
00416 
00417     void OptimizeRedundant()
00418     {
00419         // This optimization does this: min()=0, max()=0, add()=0, mul()=1
00420 
00421         if(!GetArgCount())
00422         {
00423             if(GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax)
00424                 ReplaceWithConst(0);
00425             else if(GetOp() == cMul)
00426                 ReplaceWithConst(1);
00427             return;
00428         }
00429 
00430         // And this: mul(x) = x, min(x) = x, max(x) = x, add(x) = x
00431 
00432         if(GetArgCount() == 1)
00433         {
00434             if(GetOp() == cMul || GetOp() == cAdd || GetOp() == cMin || GetOp() == cMax)
00435                 if(!getp0().getsign())
00436                     ReplaceWith(*getp0());
00437         }
00438 
00439         OptimizeDoubleNegations();
00440     }
00441 
00442     void OptimizeDoubleNegations()
00443     {
00444         if(GetOp() == cAdd)
00445         {
00446             // Eschew double negations
00447 
00448             // If any of the elements is cMul
00449             // and has a numeric constant, negate
00450             // the constant and negate sign.
00451 
00452             for(pit a=GetBegin(); a!=GetEnd(); ++a)
00453             {
00454                 SubTree &pa = *a;
00455                 if(pa.getsign()
00456                 && pa->GetOp() == cMul)
00457                 {
00458                     CodeTree &p = *pa;
00459                     for(pit b=p.GetBegin();
00460                             b!=p.GetEnd(); ++b)
00461                     {
00462                         SubTree &pb = *b;
00463                         if(pb->IsImmed())
00464                         {
00465                             pb.Negate();
00466                             pa.Negate();
00467                             break;
00468                         }
00469                     }
00470                 }
00471             }
00472         }
00473 
00474         if(GetOp() == cMul)
00475         {
00476             // If any of the elements is cPow
00477             // and has a numeric exponent, negate
00478             // the exponent and negate sign.
00479 
00480             for(pit a=GetBegin(); a!=GetEnd(); ++a)
00481             {
00482                 SubTree &pa = *a;
00483                 if(pa.getsign() && pa->GetOp() == cPow)
00484                 {
00485                     CodeTree &p = *pa;
00486                     if(p.getp1()->IsImmed())
00487                     {
00488                         // negate ok for pow when op=cImmed
00489                         p.getp1().Negate();
00490                         pa.Negate();
00491                     }
00492                 }
00493             }
00494         }
00495     }
00496 
00497     void OptimizeConstantMath1()
00498     {
00499         // This optimization does three things:
00500         //      - For adding groups:
00501         //          Constants are added together.
00502         //      - For multiplying groups:
00503         //          Constants are multiplied together.
00504         //      - For function calls:
00505         //          If all parameters are constants,
00506         //          the call is replaced with constant value.
00507 
00508         // First, do this:
00509         OptimizeAddMulFlat();
00510 
00511         switch(GetOp())
00512         {
00513             case cAdd:
00514             {
00515                 ConstList cl = BuildConstList();
00516                 FinishConst(cl);
00517                 break;
00518             }
00519             case cMul:
00520             {
00521                 ConstList cl = BuildConstList();
00522 
00523                 if(cl.value == 0.0) ReplaceWithConst(0.0);
00524                 else FinishConst(cl);
00525 
00526                 break;
00527             }
00528             #define ConstantUnaryFun(token, fun) \
00529                 case token: { const SubTree &p0 = getp0(); \
00530                     if(p0->IsImmed()) ReplaceWithConst(fun(p0->GetImmed())); \
00531                     break; }
00532             #define ConstantBinaryFun(token, fun) \
00533                 case token: { const SubTree &p0 = getp0(); \
00534                               const SubTree &p1 = getp1(); \
00535                     if(p0->IsImmed() && \
00536                        p1->IsImmed()) ReplaceWithConst(fun(p0->GetImmed(), p1->GetImmed())); \
00537                     break; }
00538 
00539             // FIXME: potential invalid parameters for functions
00540             //        can cause exceptions here
00541 
00542             ConstantUnaryFun(cAbs,   fabs);
00543             ConstantUnaryFun(cAcos,  acos);
00544             ConstantUnaryFun(cAsin,  asin);
00545             ConstantUnaryFun(cAtan,  atan);
00546             ConstantUnaryFun(cCeil,  ceil);
00547             ConstantUnaryFun(cCos,   cos);
00548             ConstantUnaryFun(cCosh,  cosh);
00549             ConstantUnaryFun(cFloor, floor);
00550             ConstantUnaryFun(cLog,   log);
00551             ConstantUnaryFun(cSin,   sin);
00552             ConstantUnaryFun(cSinh,  sinh);
00553             ConstantUnaryFun(cTan,   tan);
00554             ConstantUnaryFun(cTanh,  tanh);
00555             ConstantBinaryFun(cAtan2, atan2);
00556             ConstantBinaryFun(cMax,   Max);
00557             ConstantBinaryFun(cMin,   Min);
00558             ConstantBinaryFun(cMod,   fmod); // not a func, but belongs here too
00559             ConstantBinaryFun(cPow,   pow);
00560 
00561             case cNeg:
00562             case cSub:
00563             case cDiv:
00564                 /* Unreached (nonexistent operator)
00565                  * TODO: internal error here?
00566                  */
00567                 break;
00568 
00569             case cCot:
00570             case cCsc:
00571             case cSec:
00572             case cDeg:
00573             case cRad:
00574             case cLog10:
00575             case cSqrt:
00576             case cExp:
00577                 /* Unreached (nonexistent function)
00578                  * TODO: internal error here?
00579                  */
00580                  break;
00581         }
00582 
00583         OptimizeConflict();
00584     }
00585 
00586     void OptimizeAddMulFlat()
00587     {
00588         // This optimization flattens the topography of the tree.
00589         //   Examples:
00590         //       x + (y+z) = x+y+z
00591         //       x * (y/z) = x*y/z
00592         //       x / (y/z) = x/y*z
00593 
00594         if(GetOp() == cAdd || GetOp() == cMul)
00595         {
00596             // If children are same type as parent add them here
00597             for(pit b, a=GetBegin(); a!=GetEnd(); a=b)
00598             {
00599                 const SubTree &pa = *a;  b=a; ++b;
00600                 if(pa->GetOp() != GetOp()) continue;
00601 
00602                 // Child is same type
00603                 for(pcit c=pa->GetBegin();
00604                          c!=pa->GetEnd();
00605                          ++c)
00606                 {
00607                     const SubTree &pb = *c;
00608                     if(pa.getsign())
00609                     {
00610                         // +a -(+b +c)
00611                         // means b and c will be negated
00612 
00613                         SubTree tmp = pb;
00614                         if(GetOp() == cMul)
00615                             tmp.Invert();
00616                         else
00617                             tmp.Negate();
00618                         AddParam(tmp);
00619                     }
00620                     else
00621                         AddParam(pb);
00622                 }
00623                 Erase(a);
00624 
00625                 // Note: OptimizeConstantMath1() would be a good thing to call next.
00626             }
00627         }
00628     }
00629 
00630     void OptimizeLinearCombine()
00631     {
00632         // This optimization does the following:
00633         //
00634         //   x*x*x*x -> x^4
00635         //   x+x+x+x -> x*4
00636         //   x*x     -> x^2
00637         //   x/z/z   ->
00638         //
00639 
00640         // Remove conflicts first, so we don't have to worry about signs.
00641         OptimizeConflict();
00642 
00643         bool didchanges = false;
00644         if(GetOp() == cAdd || GetOp() == cMul)
00645         {
00646         Redo:
00647             for(pit a=GetBegin(); a!=GetEnd(); ++a)
00648             {
00649                 const SubTree &pa = *a;
00650 
00651                 list<pit> poslist;
00652 
00653                 for(pit b=a; ++b!=GetEnd(); )
00654                 {
00655                     const SubTree &pb = *b;
00656                     if(*pa == *pb)
00657                         poslist.push_back(b);
00658                 }
00659 
00660                 unsigned min = 2;
00661                 if(poslist.size() >= min)
00662                 {
00663                     SubTree arvo = pa;
00664                     bool negate = arvo.getsign();
00665 
00666                     double factor = poslist.size() + 1;
00667 
00668                     if(negate)
00669                     {
00670                         arvo.Negate();
00671                         factor = -factor;
00672                     }
00673 
00674                     CodeTree tmp(GetOp()==cAdd ? cMul : cPow,
00675                                  arvo,
00676                                  factor);
00677 
00678                     list<pit>::const_iterator j;
00679                     for(j=poslist.begin(); j!=poslist.end(); ++j)
00680                         Erase(*j);
00681                     poslist.clear();
00682 
00683                     *a = tmp;
00684                     didchanges = true;
00685                     goto Redo;
00686                 }
00687             }
00688         }
00689         if(didchanges)
00690         {
00691             // As a result, there might be need for this:
00692             OptimizeAddMulFlat();
00693             // And this:
00694             OptimizeRedundant();
00695         }
00696     }
00697 
00698     void OptimizeLogarithm()
00699     {
00700         /*
00701             This is basic logarithm math:
00702               pow(X,Y)/log(Y) = X
00703               log(X)/log(Y) = logY(X)
00704               log(X^Y)      = log(X)*Y
00705               log(X*Y)      = log(X)+log(Y)
00706               exp(log(X)*Y) = X^Y
00707 
00708             This function does these optimizations:
00709                pow(const_E, log(x))   = x
00710                pow(const_E, log(x)*y) = x^y
00711                pow(10,      log(x)*const_L10I*y) = x^y
00712                pow(z,       log(x)/log(z)*y)     = x^y
00713 
00714             And this:
00715                log(x^z)             = z * log(x)
00716             Which automatically causes these too:
00717                log(pow(const_E, x))         = x
00718                log(pow(y,       x))         = x * log(y)
00719                log(pow(pow(const_E, y), x)) = x*y
00720 
00721             And it does this too:
00722                log(x) + log(y) + log(z) = log(x * y * z)
00723                log(x * exp(y)) = log(x) + y
00724 
00725         */
00726 
00727         // Must be already in exponential form.
00728 
00729         // Optimize exponents before doing something.
00730         OptimizeExponents();
00731 
00732         if(GetOp() == cLog)
00733         {
00734             // We should have one parameter for log() function.
00735             // If we don't, we're screwed.
00736 
00737             const SubTree &p = getp0();
00738 
00739             if(p->GetOp() == cPow)
00740             {
00741                 // Found log(x^y)
00742                 SubTree p0 = p->getp0(); // x
00743                 SubTree p1 = p->getp1(); // y
00744 
00745                 // Build the new logarithm.
00746                 CodeTree tmp(GetOp(), p0);  // log(x)
00747 
00748                 // Become log(x) * y
00749                 ReplaceWith(cMul, tmp, p1);
00750             }
00751             else if(p->GetOp() == cMul)
00752             {
00753                 // Redefine &p nonconst
00754                 SubTree &p = getp0();
00755 
00756                 p->OptimizeAddMulFlat();
00757                 p->OptimizeExponents();
00758                 CHECKCONSTNEG(p, p->GetOp());
00759 
00760                 list<SubTree> adds;
00761 
00762                 for(pit b, a = p->GetBegin();
00763                            a != p->GetEnd(); a=b)
00764                 {
00765                     SubTree &pa = *a;  b=a; ++b;
00766                     if(pa->GetOp() == cPow
00767                     && pa->getp0()->IsImmed()
00768                     && pa->getp0()->GetImmed() == CONSTANT_E)
00769                     {
00770                         adds.push_back(pa->getp1());
00771                         p->Erase(a);
00772                         continue;
00773                     }
00774                 }
00775                 if(adds.size())
00776                 {
00777                     CodeTree tmp(cAdd, *this);
00778 
00779                     list<SubTree>::const_iterator i;
00780                     for(i=adds.begin(); i!=adds.end(); ++i)
00781                         tmp.AddParam(*i);
00782 
00783                     ReplaceWith(tmp);
00784                 }
00785             }
00786         }
00787         if(GetOp() == cAdd)
00788         {
00789             // Check which ones are logs.
00790             list<pit> poslist;
00791 
00792             for(pit a=GetBegin(); a!=GetEnd(); ++a)
00793             {
00794                 const SubTree &pa = *a;
00795                 if(pa->GetOp() == cLog)
00796                     poslist.push_back(a);
00797             }
00798 
00799             if(poslist.size() >= 2)
00800             {
00801                 CodeTree tmp(cMul, 1.0); // eek
00802 
00803                 list<pit>::const_iterator j;
00804                 for(j=poslist.begin(); j!=poslist.end(); ++j)
00805                 {
00806                     const SubTree &pb = **j;
00807                     // Take all of its children
00808                     for(pcit b=pb->GetBegin();
00809                              b!=pb->GetEnd();
00810                              ++b)
00811                     {
00812                         SubTree tmp2 = *b;
00813                         if(pb.getsign()) tmp2.Negate();
00814                         tmp.AddParam(tmp2);
00815                     }
00816                     Erase(*j);
00817                 }
00818                 poslist.clear();
00819 
00820                 AddParam(CodeTree(cLog, tmp));
00821             }
00822             // Done, hopefully
00823         }
00824         if(GetOp() == cPow)
00825         {
00826             const SubTree &p0 = getp0();
00827             SubTree &p1 = getp1();
00828 
00829             if(p0->IsImmed() && p0->GetImmed() == CONSTANT_E
00830             && p1->GetOp() == cLog)
00831             {
00832                 // pow(const_E, log(x)) = x
00833                 ReplaceWith(*(p1->getp0()));
00834             }
00835             else if(p1->GetOp() == cMul)
00836             {
00837                 //bool didsomething = true;
00838 
00839                 pit poslogpos; bool foundposlog = false;
00840                 pit neglogpos; bool foundneglog = false;
00841 
00842                 ConstList cl = p1->BuildConstList();
00843 
00844                 for(pit a=p1->GetBegin(); a!=p1->GetEnd(); ++a)
00845                 {
00846                     const SubTree &pa = *a;
00847                     if(pa->GetOp() == cLog)
00848                     {
00849                         if(!pa.getsign())
00850                         {
00851                             foundposlog = true;
00852                             poslogpos   = a;
00853                         }
00854                         else if(*p0 == *(pa->getp0()))
00855                         {
00856                             foundneglog = true;
00857                             neglogpos   = a;
00858                         }
00859                     }
00860                 }
00861 
00862                 if(p0->IsImmed()
00863                 && p0->GetImmed() == 10.0
00864                 && cl.value == CONSTANT_L10I
00865                 && foundposlog)
00866                 {
00867                     SubTree base = (*poslogpos)->getp0();
00868                     p1->KillConst(cl);
00869                     p1->Erase(poslogpos);
00870                     p1->OptimizeRedundant();
00871                     SubTree mul = p1;
00872 
00873                     ReplaceWith(cPow, base, mul);
00874 
00875                     // FIXME: what optimizations should be done now?
00876                     return;
00877                 }
00878 
00879                 // Put back the constant
00880                 FinishConst(cl);
00881 
00882                 if(p0->IsImmed()
00883                 && p0->GetImmed() == CONSTANT_E
00884                 && foundposlog)
00885                 {
00886                     SubTree base = (*poslogpos)->getp0();
00887                     p1->Erase(poslogpos);
00888 
00889                     p1->OptimizeRedundant();
00890                     SubTree mul = p1;
00891 
00892                     ReplaceWith(cPow, base, mul);
00893 
00894                     // FIXME: what optimizations should be done now?
00895                     return;
00896                 }
00897 
00898                 if(foundposlog
00899                 && foundneglog
00900                 && *((*neglogpos)->getp0()) == *p0)
00901                 {
00902                     SubTree base = (*poslogpos)->getp0();
00903                     p1->Erase(poslogpos);
00904                     p1->Erase(neglogpos);
00905 
00906                     p1->OptimizeRedundant();
00907                     SubTree mul = p1;
00908 
00909                     ReplaceWith(cPow, base, mul);
00910 
00911                     // FIXME: what optimizations should be done now?
00912                     return;
00913                 }
00914             }
00915         }
00916     }
00917 
00918     void OptimizeFunctionCalls()
00919     {
00920         /* Goals: sin(asin(x)) = x
00921          *        cos(acos(x)) = x
00922          *        tan(atan(x)) = x
00923          * NOTE:
00924          *   Do NOT do these:
00925          *     asin(sin(x))
00926          *     acos(cos(x))
00927          *     atan(tan(x))
00928          *   Because someone might want to wrap the angle.
00929          */
00930         // FIXME: TODO
00931     }
00932 
00933     void OptimizePowMulAdd()
00934     {
00935         // x^3 * x -> x^4
00936         // x*3 + x -> x*4
00937         // FIXME: Do those
00938 
00939         // x^1 -> x
00940         if(GetOp() == cPow)
00941         {
00942             const SubTree &base     = getp0();
00943             const SubTree &exponent = getp1();
00944 
00945             if(exponent->IsImmed())
00946             {
00947                 if(exponent->GetImmed() == 1.0)
00948                     ReplaceWith(*base);
00949                 else if(exponent->GetImmed() == 0.0
00950                      && base->NonZero())
00951                     ReplaceWithConst(1.0);
00952             }
00953         }
00954     }
00955 
00956     void OptimizeExponents()
00957     {
00958         /* Goals:
00959          *     (x^y)^z   -> x^(y*z)
00960          *     x^y * x^z -> x^(y+z)
00961          */
00962         // First move to exponential form.
00963         OptimizeLinearCombine();
00964 
00965         bool didchanges = false;
00966 
00967     Redo:
00968         if(GetOp() == cPow)
00969         {
00970             // (x^y)^z   -> x^(y*z)
00971 
00972             const SubTree &p0 = getp0();
00973             const SubTree &p1 = getp1();
00974             if(p0->GetOp() == cPow)
00975             {
00976                 CodeTree tmp(cMul, p0->getp1(), p1);
00977                 tmp.Optimize();
00978 
00979                 ReplaceWith(cPow, p0->getp0(), tmp);
00980 
00981                 didchanges = true;
00982                 goto Redo;
00983             }
00984         }
00985         if(GetOp() == cMul)
00986         {
00987             // x^y * x^z -> x^(y+z)
00988 
00989             for(pit a=GetBegin(); a!=GetEnd(); ++a)
00990             {
00991                 const SubTree &pa = *a;
00992 
00993                 if(pa->GetOp() != cPow) continue;
00994 
00995                 list<pit> poslist;
00996 
00997                 for(pit b=a; ++b != GetEnd(); )
00998                 {
00999                     const SubTree &pb = *b;
01000                     if(pb->GetOp() == cPow
01001                     && *(pa->getp0())
01002                     == *(pb->getp0()))
01003                     {
01004                         poslist.push_back(b);
01005                     }
01006                 }
01007 
01008                 if(poslist.size() >= 1)
01009                 {
01010                     poslist.push_back(a);
01011 
01012                     CodeTree base = *(pa->getp0());
01013 
01014                     CodeTree exponent(cAdd, 0.0); //eek
01015 
01016                     // Collect all exponents to cAdd
01017                     list<pit>::const_iterator i;
01018                     for(i=poslist.begin(); i!=poslist.end(); ++i)
01019                     {
01020                         const SubTree &pb = **i;
01021 
01022                         SubTree tmp2 = pb->getp1();
01023                         if(pb.getsign()) tmp2.Invert();
01024 
01025                         exponent.AddParam(tmp2);
01026                     }
01027 
01028                     exponent.Optimize();
01029 
01030                     CodeTree result(cPow, base, exponent);
01031 
01032                     for(i=poslist.begin(); i!=poslist.end(); ++i)
01033                         Erase(*i);
01034                     poslist.clear();
01035 
01036                     AddParam(result); // We're cMul, remember
01037 
01038                     didchanges = true;
01039                     goto Redo;
01040                 }
01041             }
01042         }
01043 
01044         OptimizePowMulAdd();
01045 
01046         if(didchanges)
01047         {
01048             // As a result, there might be need for this:
01049             OptimizeConflict();
01050         }
01051     }
01052 
01053     void OptimizeLinearExplode()
01054     {
01055         // x^2 -> x*x
01056         // But only if x is just a simple thing
01057 
01058         // Won't work on anything else.
01059         if(GetOp() != cPow) return;
01060 
01061         // TODO TODO TODO
01062     }
01063 
01064     void OptimizePascal()
01065     {
01066 #if 0    // Too big, too specific, etc
01067 
01068         // Won't work on anything else.
01069         if(GetOp() != cAdd) return;
01070 
01071         // Must be done after OptimizeLinearCombine();
01072 
01073         // Don't need pascal triangle
01074         // Coefficient for x^a * y^b * z^c = 3! / (a! * b! * c!)
01075 
01076         // We are greedy and want other than just binomials
01077         // FIXME
01078 
01079         // note: partial ones are also nice
01080         //     x*x + x*y + y*y
01081         //   = (x+y)^2 - x*y
01082         //
01083         //     x x * x y * + y y * +
01084         // ->  x y + dup * x y * -
01085 #endif
01086     }
01087 
01088 public:
01089 
01090     void Optimize();
01091 
01092     void Assemble(vector<unsigned> &byteCode,
01093                   vector<double>   &immed) const;
01094 
01095     void FinalOptimize()
01096     {
01097         // First optimize each parameter.
01098         for(pit a=GetBegin(); a!=GetEnd(); ++a)
01099             (*a)->FinalOptimize();
01100 
01101         /* These things are to be done:
01102          *
01103          * x * CONSTANT_DR        -> cDeg(x)
01104          * x * CONSTANT_RD        -> cRad(x)
01105          * pow(x, 0.5)            -> sqrt(x)
01106          * log(x) * CONSTANT_L10I -> log10(x)
01107          * pow(CONSTANT_E, x)     -> exp(x)
01108          * inv(sin(x))            -> csc(x)
01109          * inv(cos(x))            -> sec(x)
01110          * inv(tan(x))            -> cot(x)
01111          */
01112 
01113 
01114         if(GetOp() == cPow)
01115         {
01116             const SubTree &p0 = getp0();
01117             const SubTree &p1 = getp1();
01118             if(p0->GetOp()    == cImmed
01119             && p0->GetImmed() == CONSTANT_E)
01120             {
01121                 ReplaceWith(cExp, p1);
01122             }
01123             else if(p1->GetOp()    == cImmed
01124                  && p1->GetImmed() == 0.5)
01125             {
01126                 ReplaceWith(cSqrt, p0);
01127             }
01128         }
01129         if(GetOp() == cMul)
01130         {
01131             if(GetArgCount() == 1 && getp0().getsign())
01132             {
01133                 /***/if(getp0()->GetOp() == cSin)ReplaceWith(cCsc, getp0()->getp0());
01134                 else if(getp0()->GetOp() == cCos)ReplaceWith(cSec, getp0()->getp0());
01135                 else if(getp0()->GetOp() == cTan)ReplaceWith(cCot, getp0()->getp0());
01136             }
01137         }
01138         // Separate "if", because op may have just changed
01139         if(GetOp() == cMul)
01140         {
01141             CodeTree *found_log = 0;
01142 
01143             ConstList cl = BuildConstList();
01144 
01145             for(pit a=GetBegin(); a!=GetEnd(); ++a)
01146             {
01147                 SubTree &pa = *a;
01148                 if(pa->GetOp() == cLog && !pa.getsign())
01149                     found_log = &*pa;
01150             }
01151             if(cl.value == CONSTANT_L10I && found_log)
01152             {
01153                 // Change the log() to log10()
01154                 found_log->SetOp(cLog10);
01155                 // And forget the constant
01156                 KillConst(cl);
01157             }
01158             else if(cl.value == CONSTANT_DR)
01159             {
01160                 OptimizeRedundant();
01161                 ReplaceWith(cDeg, *this);
01162             }
01163             else if(cl.value == CONSTANT_RD)
01164             {
01165                 OptimizeRedundant();
01166                 ReplaceWith(cRad, *this);
01167             }
01168             else FinishConst(cl);
01169         }
01170 
01171         SortIfPossible();
01172     }
01173 };
01174 
01175 void CodeTreeDataPtr::Shock()
01176 {
01177  /*
01178     PrepareForWrite();
01179     paramlist &p2 = (*this)->args;
01180     for(paramlist::iterator i=p2.begin(); i!=p2.end(); ++i)
01181     {
01182         (*i)->data.Shock();
01183     }
01184  */
01185 }
01186 
01187 CodeTree::ConstList CodeTree::BuildConstList()
01188 {
01189     ConstList result;
01190     result.value     =
01191     result.voidvalue = GetOp()==cMul ? 1.0 : 0.0;
01192 
01193     list<pit> &cp = result.cp;
01194     for(pit b, a=GetBegin(); a!=GetEnd(); a=b)
01195     {
01196         SubTree &pa = *a;  b=a; ++b;
01197         if(!pa->IsImmed()) continue;
01198 
01199         double thisvalue = pa->GetImmed();
01200         if(thisvalue == result.voidvalue)
01201         {
01202             // This value is no good, forget it
01203             Erase(a);
01204             continue;
01205         }
01206         if(GetOp() == cMul)
01207             result.value *= thisvalue;
01208         else
01209             result.value += thisvalue;
01210         cp.push_back(a);
01211     }
01212     if(GetOp() == cMul)
01213     {
01214         /*
01215           Jos joku niistä arvoista on -1 eikä se ole ainoa arvo,
01216           niin joku muu niistä arvoista negatoidaan.
01217         */
01218         for(bool done=false; cp.size() > 1 && !done; )
01219         {
01220             done = true;
01221             for(list<pit>::iterator b,a=cp.begin(); a!=cp.end(); a=b)
01222             {
01223                 b=a; ++b;
01224                 if((**a)->GetImmed() == -1.0)
01225                 {
01226                     Erase(*a);
01227                     cp.erase(a);
01228 
01229                     // take randomly something
01230                     (**cp.begin())->data->NegateImmed();
01231                     if(cp.size() < 2)break;
01232                     done = false;
01233                 }
01234             }
01235         }
01236     }
01237     return result;
01238 }
01239 
01240 void CodeTree::Assemble
01241    (vector<unsigned> &byteCode,
01242     vector<double>   &immed) const
01243 {
01244     #define AddCmd(op) byteCode.push_back((op))
01245     #define AddConst(v) do { \
01246         byteCode.push_back(cImmed); \
01247         immed.push_back((v)); \
01248     } while(0)
01249 
01250     if(IsVar())
01251     {
01252         AddCmd(GetVar());
01253         return;
01254     }
01255     if(IsImmed())
01256     {
01257         AddConst(GetImmed());
01258         return;
01259     }
01260 
01261     switch(GetOp())
01262     {
01263         case cAdd:
01264         case cMul:
01265         {
01266             unsigned opcount = 0;
01267             for(pcit a=GetBegin(); a!=GetEnd(); ++a)
01268             {
01269                 const SubTree &pa = *a;
01270 
01271                 if(opcount < 2) ++opcount;
01272 
01273                 bool pnega = pa.getsign();
01274 
01275                 bool done = false;
01276                 if(pa->IsImmed())
01277                 {
01278                     if(GetOp() == cMul
01279                     && pa->data->IsInverted()
01280                     && (pnega || opcount==2)
01281                       )
01282                     {
01283                         CodeTree tmp = *pa;
01284                         tmp.data->InvertImmed();
01285                         tmp.Assemble(byteCode, immed);
01286                         pnega = !pnega;
01287                         done = true;
01288                     }
01289                     else if(GetOp() == cAdd
01290                     && (pa->data->IsNegatedOriginal()
01291                 //     || pa->GetImmed() < 0
01292                        )
01293                     && (pnega || opcount==2)
01294                            )
01295                     {
01296                         CodeTree tmp = *pa;
01297                         tmp.data->NegateImmed();
01298                         tmp.Assemble(byteCode, immed);
01299                         pnega = !pnega;
01300                         done = true;
01301                     }
01302                 }
01303                 if(!done)
01304                     pa->Assemble(byteCode, immed);
01305 
01306                 if(opcount == 2)
01307                 {
01308                     unsigned tmpop = GetOp();
01309                     if(pnega) // negate
01310                     {
01311                         tmpop = (tmpop == cMul) ? cDiv : cSub;
01312                     }
01313                     AddCmd(tmpop);
01314                 }
01315                 else if(pnega)
01316                 {
01317                     if(GetOp() == cMul) AddCmd(cInv);
01318                     else AddCmd(cNeg);
01319