15241916852

Committed 25 May 2025 08:59PM UTC coverage: 47.908% (-2.8%) from 50.743%

Build # 15241916852

Build Type

push

github

Committed by

tueda

Commit Message

ci: build arm64-windows binaries

Run Details

39009 of 81425 relevant lines covered (47.91%)

1079780.1 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

66.28

/sources/optimize.cc

/** @file optimize.cc
 *
 *        experimental routines for the optimization of FORTRAN or C output.
 */

/* #[ License : */
/*
 *   Copyright (C) 1984-2023 J.A.M. Vermaseren
 *   When using this file you are requested to refer to the publication
 *   J.A.M.Vermaseren "New features of FORM" math-ph/0010025
 *   This is considered a matter of courtesy as the development was paid
 *   for by FOM the Dutch physics granting agency and we would like to
 *   be able to track its scientific use to convince FOM of its value
 *   for the community.
 *
 *   This file is part of FORM.
 *
 *   FORM is free software: you can redistribute it and/or modify it under the
 *   terms of the GNU General Public License as published by the Free Software
 *   Foundation, either version 3 of the License, or (at your option) any later
 *   version.
 *
 *   FORM is distributed in the hope that it will be useful, but WITHOUT ANY
 *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 *   details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with FORM.  If not, see <http://www.gnu.org/licenses/>.
 */
/*
          #] License : 
          #[ includes :
*/

//#define DEBUG
//#define DEBUG_MORE
//#define DEBUG_MCTS
//#define DEBUG_GREEDY

#ifdef HAVE_CONFIG_H
#ifndef CONFIG_H_INCLUDED
#define CONFIG_H_INCLUDED
#include <config.h>
#endif
#endif

#include <vector>
#include <stack>
#include <algorithm>
#include <set>
#include <map>
#include <climits>
#include <cmath>
#include <string>
#include <algorithm>
#include <iostream>

#ifdef APPLE64
#ifndef HAVE_UNORDERED_MAP
#define HAVE_UNORDERED_MAP
#endif
#ifndef HAVE_UNORDERED_SET
#define HAVE_UNORDERED_SET
#endif
#endif

#ifdef HAVE_UNORDERED_MAP
        #include <unordered_map>
        using std::unordered_map;
#elif !defined(HAVE_TR1_UNORDERED_MAP) && defined(HAVE_BOOST_UNORDERED_MAP_HPP)
        #include <boost/unordered_map.hpp>
        using boost::unordered_map;
#else
        #include <tr1/unordered_map>
        using std::tr1::unordered_map;
#endif
#ifdef HAVE_UNORDERED_SET
        #include <unordered_set>
        using std::unordered_set;
#elif !defined(HAVE_TR1_UNORDERED_SET) && defined(HAVE_BOOST_UNORDERED_SET_HPP)
        #include <boost/unordered_set.hpp>
        using boost::unordered_set;
#else
        #include <tr1/unordered_set>
        using std::tr1::unordered_set;
#endif

#if defined(HAVE_BUILTIN_POPCOUNT)
        static inline int popcount(unsigned int x) {
                return __builtin_popcount(x);
        }
#elif defined(HAVE_POPCNT)
        #include <intrin.h>
        static inline int popcount(unsigned int x) {
                return __popcnt(x);
        }
#else
        static inline int popcount(unsigned int x) {
                int count = 0;
                while (x > 0) {
                        if ((x & 1) == 1) count++;
                        x >>= 1;
                }
                return count;
        }
#endif

extern "C" {
#include "form3.h"
}

//#ifdef DEBUG
#include "mytime.h"
//#endif

using namespace std;

// operators
const WORD OPER_ADD = -1;
const WORD OPER_MUL = -2;
const WORD OPER_COMMA = -3;

// class for a node of the MCTS tree
class tree_node {
public:
        vector<tree_node> childs;
        double sum_results;
        int num_visits;
        WORD var;
        bool finished;
        PADPOINTER(1,1,1,1);

        tree_node (int _var=0):
                sum_results(0), num_visits(0), var(_var), finished(false) {}

        tree_node(const tree_node& other):
                childs(other.childs),
                sum_results(other.sum_results),
                num_visits(other.num_visits),
                var(other.var),
                finished(other.finished) {}

        tree_node& operator=(const tree_node& other) {
                if (this != &other) {
                        childs = other.childs;
                        sum_results = other.sum_results;
                        num_visits = other.num_visits;
                        var = other.var;
                        finished = other.finished;
                }
                return *this;
        }
};

// global variables for multithreading
WORD *optimize_expr;
vector<vector<WORD> > optimize_best_Horner_schemes;
int optimize_num_vars;
int optimize_best_num_oper;
vector<WORD> optimize_best_instr;
vector<WORD> optimize_best_vars;

// global variables for MCTS
bool mcts_factorized, mcts_separated;
vector<WORD> mcts_vars;
tree_node mcts_root;
int mcts_expr_score;
set<pair<int,vector<WORD> > > mcts_best_schemes;

#ifdef WITHPTHREADS
pthread_mutex_t optimize_lock;
#endif

/*
          #] includes : 
          #[ print_instr :
*/

void print_instr (const vector<WORD> &instr, WORD num)
{
        const WORD *tbegin = &*instr.begin();
        const WORD *tend = tbegin+instr.size();
        for (const WORD *t=tbegin; t!=tend; t+=*(t+2)) {
                MesPrint("<%d> %a",num,t[2],t);
        }
}

/*
          #] print_instr : 
          #[ my_random_shuffle :
*/

/**  Random shuffle
 *
 *         Description
 *         ===========
 *         Randomly permutes elements in the range [fr,to). Functionality is
 *         the same as C++'s "random_shuffle", but it uses Form's "wranf".
 */
template <class RandomAccessIterator>
void my_random_shuffle (PHEAD RandomAccessIterator fr, RandomAccessIterator to) {
  for (int i=to-fr-1; i>0; --i)
                swap (fr[i],fr[wranf(BHEAD0) % (i+1)]);
}

/*
          #] my_random_shuffle : 
          #[ get_expression :
*/

static WORD comlist[3] = {TYPETOPOLYNOMIAL,3,DOALL};

/**  Get expression
 *
 *         Description
 *         ===========
 *         Reads an expression from the input file into a buffer (called
 *         optimize_expr). This buffer is used during the optimization
 *         process. Non-symbols are removed by ConvertToPoly and are put in
 *         temporary symbols.
 *
 *         The return value is the length of the expression in WORDs, or a
 *         negative number if it failed.
 */
LONG get_expression (int exprnr) {

        GETIDENTITY;

  AR.NoCompress = 1;

  NewSort(BHEAD0);
  EXPRESSIONS e = Expressions+exprnr;
  SetScratch(AR.infile,&(e->onfile));

  // get header term
  WORD *term = AT.WorkPointer;
        GetTerm(BHEAD term);

  NewSort(BHEAD0);

  // get terms
  while (GetTerm(BHEAD term) > 0) {
        AT.WorkPointer = term + *term;
        WORD *t1 = term;
        WORD *t2 = term + *term;
        if (ConvertToPoly(BHEAD t1,t2,comlist,1) < 0) return -1;
        int n = *t2;
        NCOPY(t1,t2,n);
        AT.WorkPointer = term + *term;
        if (StoreTerm(BHEAD term)) return -1;
  }

  // sort and store in buffer
  LONG len = EndSort(BHEAD (WORD *)((VOID *)(&optimize_expr)),2);
  LowerSortLevel();
  AT.WorkPointer = term;

        return len;
}

/*
          #] get_expression : 
          #[ PF_get_expression :
*/
#ifdef WITHMPI

// get_expression for ParFORM
LONG PF_get_expression (int exprnr) {
        LONG len;
        if (PF.me == MASTER) {
                len = get_expression(exprnr);
        }
        if (PF.numtasks > 1) {
                PF_BroadcastBuffer(&optimize_expr, &len);
        }
        return len;
}

// replace get_expression called in Optimize
#define get_expression PF_get_expression

#endif
/*
          #] PF_get_expression : 
          #[ get_brackets :
*/

/**  Get brackets
 *
 *         Description
 *         ===========
 *         Checks whether the input expression (stored in optimize_expr)
 *         contains brackets. If so, this method replaces terms outside
 *         brackets by powers of SEPERATESYMBOL (equal brackets have equal
 *         powers) and the brackets are returned. If not, the result is
 *         empty.
 *
 *         Brackets are used for simultaneous optimization. The symbol
 *         SEPARATESYMBOL is always the first one used in a Horner scheme.
 */

vector<vector<WORD> > get_brackets () {

  // check for brackets in expression
  bool has_brackets = false;
  for (WORD *t=optimize_expr; *t!=0; t+=*t) {
        WORD *tend=t+*t; tend-=ABS(*(tend-1));
        for (WORD *t1=t+1; t1<tend; t1+=*(t1+1))
          if (*t1 == HAAKJE)
                has_brackets=true;
  }

  // replace brackets by SEPARATESYMBOL
  vector<vector<WORD> > brackets;

  if (has_brackets) {
                int exprlen=10;        // we need potential space for an empty bracket
                for (WORD *t=optimize_expr; *t!=0; t+=*t)
                        exprlen += *t;
        WORD *newexpr = (WORD *)Malloc1(exprlen*sizeof(WORD), "optimize newexpr");

        int i=0;
        int sep_power = 0;

        for (WORD *t=optimize_expr; *t!=0; t+=*t) {
          WORD *t1 = t+1;

                        // scan for bracket
          vector<WORD> bracket;
          for (; *t1!=HAAKJE; t1+=*(t1+1))
                bracket.insert(bracket.end(), t1, t1+*(t1+1));

          if (brackets.size()==0 || bracket!=brackets.back()) {
                sep_power++;
                brackets.push_back(bracket);
          }
          t1+=*(t1+1);

          WORD left = t + *t - t1;
          bool more_symbols = (left != ABS(*(t+*t-1)));

                        // add power of SEPARATESYMBOL
          newexpr[i++] = 1 + left + (more_symbols ? 2 : 4);
          newexpr[i++] = SYMBOL;
          newexpr[i++] = (more_symbols ? *(t1+1) + 2 : 4);
          newexpr[i++] = SEPARATESYMBOL;
          newexpr[i++] = sep_power;

                        // add remaining terms
          if (more_symbols) {
                t1+=2;
                left-=2;
          }
          while (left-->0)
                newexpr[i++] = *(t1++);
        }
/*
        We insert here an empty bracket that is zero.
        It is used for the case that there is only a single bracket which is
        outside the notation for trees at a later stage.

        There may be a problem now in that in the case of sep_power==1
        newexpr is bigger than optimize_expr. We have to check that.
*/
        if ( sep_power == 1 )
        {
          WORD *t;
          vector<WORD> bracket;
          bracket.push_back(0);
          bracket.push_back(0);
          bracket.push_back(0);
          bracket.push_back(0);
          sep_power++;
          brackets.push_back(bracket);
          newexpr[i++] = 8;
          newexpr[i++] = SYMBOL;
          newexpr[i++] = 4;
          newexpr[i++] = SEPARATESYMBOL;
          newexpr[i++] = sep_power;
          newexpr[i++] = 1;
          newexpr[i++] = 1;
          newexpr[i++] = 3;
          newexpr[i++] = 0;
          for (t=optimize_expr; *t!=0; t+=*t) {}
          if ( t-optimize_expr+1 < i ) {  // We have to redo this
                M_free(optimize_expr,"$-sort space");
                optimize_expr = (WORD *)Malloc1(i*sizeof(WORD),"$-sort space");
          }
        }
        else {
                newexpr[i++] = 0;
        }
        memcpy(optimize_expr, newexpr, i*sizeof(WORD));
        M_free(newexpr,"optimize newexpr");

        // if factorized, replace SEP by FAC and remove brackets
        if (brackets[0].size()>0 && brackets[0][2]==FACTORSYMBOL) {
                for (WORD *t=optimize_expr; *t!=0; t+=*t) {
                        if (*t == ABS(*(t+*t-1))+1) continue;
                        if (t[1]==SYMBOL)
                                for (int i=3; i<t[2]; i+=2)
                                        if (t[i]==SEPARATESYMBOL) t[i]=FACTORSYMBOL;
                }
                return vector<vector<WORD> >();
        }
  }

  return brackets;
}

/*
          #] get_brackets : 
          #[ count_operators :
*/

/**  Count operators
 *
 *         Description
 *         ===========
 *         Counts the number of operators in a Form-style expression.
 */
int count_operators (const WORD *expr, bool print=false) {

        int n=0;
        while (*(expr+n)!=0) n+=*(expr+n);

        int cntpow=0, cntmul=0, cntadd=0, sumpow=0;
        WORD maxpowfac=1, maxpowsep=1;

        for (const WORD *t=expr; *t!=0; t+=*t) {
                if (t!=expr) cntadd++;                                // new term
                if (*t==ABS(*(t+*t-1))+1) continue; // only coefficient

                int cntsym=0;

                if (t[1]==SYMBOL)
                        for (int i=3; i<t[2]; i+=2) {
                                if (t[i]==FACTORSYMBOL) {
                                        maxpowfac = max(maxpowfac, t[i+1]);
                                        continue;
                                }
                                if (t[i]==SEPARATESYMBOL) {
                                        maxpowsep = max(maxpowsep, t[i+1]);
                                        continue;
                                }
                                if (t[i+1]>2) {                  // (extra)symbol power>2
                                        cntpow++;
                                        sumpow += (int)floor(log(t[i+1])/log(2.0)) + popcount(t[i+1]) - 1;
                                }
                                if (t[i+1]==2) cntmul++; // (extra)symbol squared
                                cntsym++;
                        }

                if (ABS(*(t+*t-1))!=3 || *(t+*t-2)!=1 || *(t+*t-3)!=1) cntsym++; // non +/-1 coefficient

                if (cntsym > 0)        cntmul+=cntsym-1;
        }

        cntadd -= maxpowfac-1;
        cntmul += maxpowfac-1;

        cntadd -= maxpowsep-1;

        if (print)
                MesPrint ("*** STATS: original        %lP %lM %lA : %l", cntpow,cntmul,cntadd,sumpow+cntmul+cntadd);

        return sumpow+cntmul+cntadd;
}

/**  Count operators
 *
 *         Description
 *         ===========
 *         Counts the number of operators in a vector of instructions
 */
int count_operators (const vector<WORD> &instr, bool print=false) {

        int cntpow=0, cntmul=0, cntadd=0, sumpow=0;

        const WORD *ebegin = &*instr.begin();
        const WORD *eend = ebegin+instr.size();

        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                for (const WORD *t=e+3; *t!=0; t+=*t) {
                        if (t!=e+3) {
                                if (*(e+1)==OPER_ADD) cntadd++;         // new term
                                if (*(e+1)==OPER_MUL) cntmul++;         // new term
                        }
                        if (*t == ABS(*(t+*t-1))+1) continue;        // only coefficient
                        if (*(t+1)==SYMBOL || *(t+1)==EXTRASYMBOL) {
                                if (*(t+4)==2) cntmul++;                        // (extra)symbol squared
                                if (*(t+4)>2) {                                         // (extra)symbol power>2
                                        cntpow++;
                                        sumpow += (int)floor(log(*(t+4))/log(2.0)) + popcount(*(t+4)) - 1;
                                }
                        }
                        if (ABS(*(t+*t-1))!=3 || *(t+*t-2)!=1 || *(t+*t-3)!=1) cntmul++; // non +/-1 coefficient
                }
        }

        if (print)
                MesPrint ("*** STATS: optimized %lP %lM %lA : %l", cntpow,cntmul,cntadd,sumpow+cntmul+cntadd);

        return sumpow+cntmul+cntadd;
}

/*
          #] count_operators : 
          #[ occurrence_order :
*/

/**  Occurrence order
 *
 *         Description
 *         ===========
 *         Extracts all variables from an expression and sorts them with
 *         most occurring first (or last, with rev=true)
 */
vector<WORD> occurrence_order (const WORD *expr, bool rev) {

        // count the number of occurrences of variables
        map<WORD,int> cnt;
        for (const WORD *t=expr; *t!=0; t+=*t) {
                if (*t == ABS(*(t+*t-1))+1) continue; // skip constant term
                if (t[1] == SYMBOL)
                        for (int i=3; i<t[2]; i+=2)
                                cnt[t[i]]++;
        }

        bool is_fac=false, is_sep=false;
        if (cnt.count(FACTORSYMBOL)) {
                is_fac=true;
                cnt.erase(FACTORSYMBOL);
        }
        if (cnt.count(SEPARATESYMBOL)) {
                is_sep=true;
                cnt.erase(SEPARATESYMBOL);
        }

        // determine the order of the variables
        vector<pair<int,WORD> > cnt_order;
        for (map<WORD,int>::iterator i=cnt.begin(); i!=cnt.end(); i++)
                cnt_order.push_back(make_pair(i->second, i->first));
        sort(cnt_order.rbegin(), cnt_order.rend());

        // create resulting order
        vector<WORD> order;
        for (int i=0; i<(int)cnt_order.size(); i++)
                order.push_back(cnt_order[i].second);

        if (rev) reverse(order.begin(),order.end());

        // add FACTORSYMBOL/SEPARATESYMBOL
        if (is_fac) order.insert(order.begin(), FACTORSYMBOL);
        if (is_sep) order.insert(order.begin(), SEPARATESYMBOL);

        return order;
}

/*
          #] occurrence_order : 
          #[ Horner_tree :
*/

/**  Horner tree building
 *
 *         Description
 *         ===========
 *         Given a Form-style expression (in a buffer in memory), this
 *         builds an expression tree. The tree is determined by a
 *         multivariate Horner scheme, i.e., something of the form:
 *
 *                1+y+x*(2+y*(1+y)+x*(3-y*(...)))
 *
 *         The order of the variables is given to the method "Horner_tree",
 *         which renumbers ad reorders the terms of the expression. Next,
 *         the recursive method "build_Horner_tree" does the actual tree
 *         construction.
 *
 *         The tree is represented in postfix notation. Tokens are of the
 *         following forms:
 *
 *         - SNUMBER tokenlength num den coefflength
 *         - SYMBOL tokenlength variable power
 *         - OPER_ADD or OPER_MUL
 *
 *         Note
 *         ====
 *         Sets AN.poly_num_vars and allocates AN.poly_vars. The latter
 *         should be freed later.
 */

/**  Get power of variable (helper function for build_Horner_tree)
 *
 *         Description
 *         ===========
 *         Returns the power of the variable "var", which is at position
 *         "pos" in this term, if it is present.
 */
WORD getpower (const WORD *term, int var, int pos) {
        if (*term == ABS(*(term+*term-1))+1) return 0; // constant term
        if (2*pos+2 >= term[2]) return 0;                           // too few symbols
        if (term[2*pos+3] != var) return 0;                    // incorrect symbol
        return term[2*pos+4];                                                   // return power
}

/**  Call GcdLong/DivLong with leading zeroes
 *
 *         Description
 *         ===========
 *         These method remove leading zeroes, so that GcdLong and DivLong
 *         can safely be called.
 */
void fixarg (UWORD *t, WORD &n) {
        int an=ABS(n), sn=SGN(n);
        while (*(t+an-1)==0) an--;
        n=an*sn;
}

void GcdLong_fix_args (PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc) {
        fixarg(a,na);
        fixarg(b,nb);
        GcdLong(BHEAD a,na,b,nb,c,nc);
}

void DivLong_fix_args(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c,        WORD *nc, UWORD *d, WORD *nd) {
        fixarg(a,na);
        fixarg(b,nb);
        DivLong(a,na,b,nb,c,nc,d,nd);
}

/**  Build the Horner tree
 *
 *         Description
 *         ===========
 *         Constructs the Horner tree. The method processes one variable and
 *         continues recursively until the Horner scheme is completed.
 *
 *         "terms" is a pointer to the starts of the terms. "numterms" is
 *         the number of terms to be processed. "var" is the next variable
 *         to be processed (index between 0 and #maxvar) and "maxvar" is the
 *         last variable to be processed, so that partial Horner trees can
 *         also be constructed. "pos" is the position that the power of
 *         "var" should be in (one level further in the recursion, "pos" has
 *         increased by 0 or 1 depending on whether the previous power was 0
 *         or not). The result is written at the pointer "res".
 *
 *         This method also factors out gcds of the coefficients. The result
 *         should end with "gcd OPER_MUL" at all times, so that one level
 *         higher gcds can be extracted again.
 */
void build_Horner_tree (const WORD **terms, int numterms, int var, int maxvar, int pos, vector<WORD> *res) {

        GETIDENTITY;

        if (var == maxvar) {
                // Horner tree is finished, so add remaining terms unfactorized
                // (note: since only complete Horner schemes seem to be useful, numterms=1 here

                for (int fr=0; fr<numterms; fr++) {

                        bool empty = true;
                        const WORD *t = terms[fr];

                        // add symbols
                        if (*t != ABS(*(t+*t-1))+1)
                                for (int i=3+2*pos; i<t[2]; i+=2) {
                                        res->push_back(SYMBOL);
                                        res->push_back(4);
                                        res->push_back(t[i]);
                                        res->push_back(t[i+1]);
                                        if (!empty) res->push_back(OPER_MUL);
                                        empty = false;
                                }

                        // if empty, add a 1, since the result should look like "... coeff *"
                        if (empty) {
                                res->push_back(SNUMBER);
                                res->push_back(5);
                                res->push_back(1);
                                res->push_back(1);
                                res->push_back(3);
                        }

                        // add coefficient
                        res->push_back(SNUMBER);
                        WORD n = ABS(*(t+*t-1));
                        res->push_back(n+2);
                        for (int i=0; i<n; i++)
                                res->push_back(*(t+*t-n+i));
                        res->push_back(OPER_MUL);

                        if (fr>0) res->push_back(OPER_ADD);
                }

                // result should end with gcd of the terms; right now this never
                // triggers, but if one would allow for incomplete Horner schemes,
                // one should extract the gcd here
                if (numterms > 1) {
                        res->push_back(SNUMBER);
                        res->push_back(5);
                        res->push_back(1);
                        res->push_back(1);
                        res->push_back(3);
                        res->push_back(OPER_MUL);
                }
        }
        else {
                // extract variable "var" and the gcd and proceed recursively

                WORD nnum = 0, nden = 0, ntmp = 0, ndum = 0;
                UWORD *num = NumberMalloc("build_Horner_tree");
                UWORD *den = NumberMalloc("build_Horner_tree");
                UWORD *tmp = NumberMalloc("build_Horner_tree");
                UWORD *dum = NumberMalloc("build_Horner_tree");

                // previous coefficient for gcd extraction or coefficient multiplication
                int prev_coeff_idx = -1;

                for (int fr=0; fr<numterms;) {

                        // find power of current term
                        WORD pow = getpower(terms[fr], var, pos);

                        // find all terms with that power
                        int to=fr+1;
                        while (to<numterms && getpower(terms[to], var, pos) == pow) to++;

                        // recursively build Horner tree of all terms proportional to var^pow
                        build_Horner_tree (terms+fr, to-fr, var+1, maxvar, pow==0?pos:pos+1, res);

                        if (AN.poly_vars[var] != FACTORSYMBOL && AN.poly_vars[var] != SEPARATESYMBOL) {
                                // if normal symbol, find gcd(numerators) and gcd(denominators)
                                WORD  n1 = res->at(res->size()-2) / 2;
                                WORD *t1 = &res->at(res->size()-2-2*ABS(n1));

                                WORD *t2 = fr==0 ? t1 : &res->at(prev_coeff_idx);
                                WORD  n2 = fr==0 ? n1 : *(t2+*(t2+1)-1) / 2;
                                if (fr>0) t2+=2;

                                GcdLong_fix_args(BHEAD (UWORD *)t1,n1,(UWORD *)t2,n2,num,&nnum);
                                GcdLong_fix_args(BHEAD (UWORD *)t1+ABS(n1),ABS(n1),(UWORD *)t2+ABS(n2),ABS(n2),den,&nden);

                                // divide out gcds; note: leading zeroes can be added here
                                for (int i=0; i<2; i++) {
                                        if (i==1 && fr==0) break;

                                        WORD *t = i==0 ? t1 : t2;
                                        WORD n = i==0 ? n1 : n2;

                                        DivLong_fix_args((UWORD *)t, n, num, nnum, tmp, &ntmp, dum, &ndum);
                                        for (int j=0; j<ABS(ntmp); j++) *t++ = tmp[j];
                                        for (int j=0; j<ABS(n)-ABS(ntmp); j++) *t++ = 0;

                                        if (SGN(ntmp) != SGN(n)) n=-n;

                                        DivLong_fix_args((UWORD *)t, n, den, nden, tmp, &ntmp, dum, &ndum);
                                        for (int j=0; j<ABS(ntmp); j++) *t++ = tmp[j];
                                        for (int j=0; j<ABS(n)-ABS(ntmp); j++) *t++ = 0;

                                        *t++ = SGN(n) * (2*ABS(n)+1);
                                }

                                // add the addition operator
                                if (fr>0) res->push_back(OPER_ADD);

                                // add the power of "var"
                                WORD nextpow = (to==numterms ? 0 : getpower(terms[to], var, pos));

                                if (pow-nextpow > 0) {
                                        res->push_back(SYMBOL);
                                        res->push_back(4);
                                        res->push_back(var);
                                        res->push_back(pow-nextpow);
                                        res->push_back(OPER_MUL);
                                }

                                // add the extracted gcd
                                res->push_back(SNUMBER);
                                WORD n = MaX(ABS(nnum),nden);
                                res->push_back(n*2+3);
                                for (int i=0; i<ABS(nnum); i++) res->push_back(num[i]);
                                for (int i=0; i<n-ABS(nnum); i++) res->push_back(0);
                                for (int i=0; i<nden; i++) res->push_back(den[i]);
                                for (int i=0; i<n-ABS(nden); i++) res->push_back(0);
                                res->push_back(SGN(nnum)*(2*n+1));
                                res->push_back(OPER_MUL);

                                prev_coeff_idx = res->size() - ABS(res->at(res->size()-2)) - 3;
                        }
                        else if (AN.poly_vars[var]==FACTORSYMBOL) {

                                // if factorsymbol, multiply overall integer contents

                                if (fr>0) {
                                        WORD  n1 = res->at(res->size()-2) / 2;
                                        WORD *t1 = &res->at(res->size()-2-2*ABS(n1));
                                        WORD *t2 = &res->at(prev_coeff_idx);
                                        WORD  n2 = *(t2+*(t2+1)-1) / 2;
                                        t2+=2;

                                        MulRat(BHEAD (UWORD *)t1,n1,(UWORD *)t2,n2,tmp,&ntmp);

                                        // replace previous coefficient with 1
                                        n2=ABS(n2);
                                        for (int i=0; i<ABS(n2); i++)
                                                t2[i] = t2[n2+i] = i==0 ? 1 : 0;
                                        t2[2*n2] = 2*n2+1;

                                        // remove this coefficient
                                        for (int i=0; i<2*ABS(n1)+2; i++)
                                                res->pop_back();

                                        // add product
                                        res->back() = 2*ABS(ntmp)+3;                                   // adjust size of term
                                        res->insert(res->end(), tmp, tmp+2*ABS(ntmp)); // num/den coefficient
                                        res->push_back(SGN(ntmp)*(2*ABS(ntmp)+1));           // size of coefficient
                                        res->push_back(OPER_MUL);                                           // operator
                                }

                                prev_coeff_idx = res->size() - ABS(res->at(res->size()-2)) - 3;

                                // multiply previous factors with this factor
                                if (fr>0)
                                        res->push_back(OPER_MUL);
                        }
                        else { // AN.poly_vars[var]==SEPARATESYMBOL
                                if (fr>0)
                                        res->push_back(OPER_COMMA);
                                prev_coeff_idx = -1;
                        }

                        fr=to;
                }

                // cleanup
                NumberFree(dum,"build_Horner_tree");
                NumberFree(tmp,"build_Horner_tree");
                NumberFree(den,"build_Horner_tree");
                NumberFree(num,"build_Horner_tree");
        }
}

/**  Term compare (helper function for Horner_tree)
 *
 *         Description
 *         ===========
 *         Compares two terms of the form "L SYM 4 x n coeff" or "L
 *         coeff". Lower powers of lower-indexed symbols come first. This is
 *         used to sort the terms in correct order.
 */
bool term_compare (const WORD *a, const WORD *b) {
        if (*a == ABS(*(a+*a-1))+1) return true; // coefficient comes first
        if (*b == ABS(*(b+*b-1))+1) return false;
        if (a[1]!=SYMBOL) return true;
        if (b[1]!=SYMBOL) return false;
        for (int i=3; i<a[2] && i<b[2]; i+=2) {
                if (a[i]  !=b[i]  ) return a[i]  >b[i];
                if (a[i+1]!=b[i+1]) return a[i+1]<b[i+1];
        }
        return a[2]<b[2];
}

/**  Prepare Horner tree building
 *
 *         Description
 *         ===========
 *         This method renumbers the variables to 0...#vars-1 according to
 *         the specified order. Next, it stored pointer to individual terms
 *         and sorts the terms with higher power first. Then the sorted
 *         lists of power is used for the construction of the Horner tree.
 */
vector<WORD> Horner_tree (const WORD *expr, const vector<WORD> &order) {
#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: Horner_tree(%a)", thetime_str().c_str(), order.size(), &order[0]);
#endif

        GETIDENTITY;

        // find the renumbering scheme (new numbers are 0,1,...,#vars-1)
        map<WORD,WORD> renum;
        AN.poly_num_vars = order.size();
        AN.poly_vars = (WORD *)Malloc1(AN.poly_num_vars*sizeof(WORD), "AN.poly_vars");
        for (int i=0; i<AN.poly_num_vars; i++) {
                AN.poly_vars[i] = order[i];
                renum[order[i]] = i;
        }

        // sort variables in individual terms using bubble sort
        WORD* sorted;
        WORD* dynamicAlloc = 0;
        LONG sumsize = 0;

        for (const WORD *t=expr; *t!=0; t+=*t) {
                sumsize += *t;
        }

        // We actually need sumsize + 1 WORDS available, due to the "*sorted = 0;"
        // at the end of the sort loop.
        if ( AT.WorkPointer + sumsize + 1 > AT.WorkTop ) {
                dynamicAlloc = (WORD*)Malloc1(sizeof(*dynamicAlloc)*(sumsize+1), "Horner_tree alloc");
                sorted = dynamicAlloc;
        }
        else {
                sorted = AT.WorkPointer;
        }

        for (const WORD *t=expr; *t!=0; t+=*t) {
                memcpy(sorted, t, *t*sizeof(WORD));

                if (*t != ABS(*(t+*t-1))+1) {
                        // non-constant term

                        // renumber variables, adding new variables if necessary
                        for (int i=3; i<sorted[2]; i+=2) {
                                if (!renum.count(sorted[i])) {
                                        renum[sorted[i]] = AN.poly_num_vars;

                                        WORD *new_poly_vars = (WORD *)Malloc1((AN.poly_num_vars+1)*sizeof(WORD), "AN.poly_vars");
                                        memcpy(new_poly_vars, AN.poly_vars, AN.poly_num_vars*sizeof(WORD));
                                        M_free(AN.poly_vars,"poly_vars");
                                        AN.poly_vars = new_poly_vars;
                                        AN.poly_vars[AN.poly_num_vars] = sorted[i];
                                        AN.poly_num_vars++;
                                }
                                sorted[i] = renum[sorted[i]];
                        }
                        // order variables
                        for (int i=0; i<sorted[2]/2; i++)
                                for (int j=3; j+2<sorted[2]; j+=2)
                                        if (sorted[j] > sorted[j+2]) {
                                                swap(sorted[j]        , sorted[j+2]);
                                                swap(sorted[j+1], sorted[j+3]);
                                        }
                }

                sorted += *sorted;
        }

        *sorted = 0;
        if ( dynamicAlloc ) {
                sorted = dynamicAlloc;
        }
        else {
                sorted = AT.WorkPointer;
        }

        // find pointers to all terms and sort them efficiently
        vector<const WORD *> terms;
        for (const WORD *t=sorted; *t!=0; t+=*t)
                terms.push_back(t);
        sort(terms.rbegin(),terms.rend(),term_compare);

        // construct the Horner tree
        vector<WORD> res;
        build_Horner_tree(&terms[0], terms.size(), 0, AN.poly_num_vars, 0, &res);

        // remove leading zeroes in coefficients
        int j=0;
        for (int i=0; i<(int)res.size();) {
                if (res[i]==OPER_ADD || res[i]==OPER_MUL || res[i]==OPER_COMMA)
                        res[j++] = res[i++];
                else if (res[i]==SYMBOL) {
                        memmove(&res[j], &res[i], res[i+1]*sizeof(WORD));
                        i+=res[j+1];
                        j+=res[j+1];
                }
                else if (res[i]==SNUMBER) {
                        int n = (res[i+1]-2)/2;
                        int dn = 0;
                        while (res[i+1+n-dn]==0 && res[i+1+2*n-dn]==0) dn++;
                        res[j++] = SNUMBER;
                        res[j++] = 2*(n-dn) + 3;
                        memmove(&res[j], &res[i+2], (n-dn)*sizeof(WORD));
                        j+=n-dn;
                        memmove(&res[j], &res[i+n+2], (n-dn)*sizeof(WORD));
                        j+=n-dn;
                        res[j++] = SGN(res[i+2*n+2])*(2*(n-dn)+1);
                        i+=2*n+3;
                }
        }
        res.resize(j);

        if ( dynamicAlloc ) {
                M_free(dynamicAlloc, "Horner_tree alloc");
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: Horner_tree(%a)", thetime_str().c_str(), order.size(), &order[0]);
#endif

        return res;
}

/*
          #] Horner_tree : 
          #[ print_tree :
*/

// print Horner tree (for debugging)
void print_tree (const vector<WORD> &tree) {

        GETIDENTITY;

        for (int i=0; i<(int)tree.size();) {
                if (tree[i]==OPER_ADD) {
                        MesPrint("+%");
                        i++;
                }
                else if (tree[i]==OPER_MUL) {
                        MesPrint("*%");
                        i++;
                }
                else if (tree[i]==OPER_COMMA) {
                        MesPrint(",%");
                        i++;
                }
                else if (tree[i]==SNUMBER) {
                        UBYTE buf[100];
                        int n = tree[i+tree[i+1]-1]/2;
                        PrtLong((UWORD *)&tree[i+2], n, buf);
                        int l = strlen((char *)buf);
                        buf[l]='/';
                        n=ABS(n);
                        PrtLong((UWORD *)&tree[i+2+n], n, buf+l+1);
                        MesPrint("%s%",buf);
                        i+=tree[i+1];
                }
                else if (tree[i]==SYMBOL) {
                        if (AN.poly_vars[tree[i+2]] < 10000)
                                MesPrint("%s^%d%", VARNAME(symbols,AN.poly_vars[tree[i+2]]), tree[i+3]);
                        else
                                MesPrint("Z%d^%d%", MAXVARIABLES-AN.poly_vars[tree[i+2]], tree[i+3]);
                        i+=tree[i+1];
                }
                else {
                        MesPrint("error");
                        exit(1);
                }

                MesPrint(" %");
        }

        MesPrint("");
}

/*
          #] print_tree : 
          #[ generate_instructions :
*/


struct CSEHash {
        size_t operator()(const vector<WORD>& n) const {
                return n[0];
        }
};

struct CSEEq {
        bool operator()(const vector<WORD>& lhs, const vector<WORD>& rhs) const {
                        if (lhs.size() != rhs.size()) return false;
                        for (unsigned int i = 1; i < lhs.size(); i++) {
                                if (lhs[i] != rhs[i]) return false;
                        }
                        return true;
        }
};


template<typename T> size_t hash_range(T* array, int size) {
        size_t hash = 0;

        for (int i = 0; i < size; i++) {
                hash ^= array[i] + 0x9e3779b9 + (hash << 6) + (hash >> 2);
        }

        return hash;
}

/**  Generate instructions
 *
 *         Description
 *         ===========
 *         Converts the expression tree to a list of instructions that
 *         directly translate to code. Instructions are of the form:
 *
 *                expr.nr operator length [operands]+ trailing.zero
 *
 *         The operands are of the form:
 *
 *                length [(EXTRA)SYMBOL length variable power] coeff
 *
 *         This method only generates binary operators. Merging is done
 *         later. The method also checks for common subexpressions and
 *         eliminates them and the flag "do_CSE" is set.
 *
 *         Implementation details
 *         ======================
 *         The map "ID" keeps track of which subexpressions already
 *         exist. The key is formatted as one of the following:
 *
 *                SYMBOL x n
 *                SNUMBER coeff
 *                OPERATOR LHS RHS
 *
 *         with LHS/RHS formatted as one of the following:
 *
 *                SNUMBER idx 0
 *                (EXTRA)SYMBOL x n
 *
 *         ID[symbol] or ID[operator] equals a subexpression
 *         number. ID[coeff] equals the position of the number in the input.
 *
 *         The stack s is used the process the postfix expression
 *         tree. Three-word tokens of the form:
 *
 *                SNUMBER idx.of.coeff 0
 *                SYMBOL x n
 *                EXTRASYMBOL x 1
 *
 *         are pushed onto it. Operators pop two operands and push the
 *         resulting expression.
 *
 *         (Extra)symbols are 1-indexed, because -X is also needed to
 *         represent -1 times this term.
 *
 *         There is currently a bug. The notation cannot tell if there is a single
 *         bracket and then ignores the bracket.
 *
 *         TODO: check if this method performs properly if do_CSE=false
 */
vector<WORD> generate_instructions (const vector<WORD> &tree, bool do_CSE) {
#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: generate_instructions(cse=%d)",
                                                thetime_str().c_str(), do_CSE?1:0);
#endif

        typedef unordered_map<vector<WORD>, int, CSEHash, CSEEq> csemap;
        csemap ID;

        // reserve lots of space, to prevent later rehashes
        // TODO: what if this is too large? make a parameter?
        if (do_CSE) {
                        ID.rehash(mcts_expr_score * 2);
        }

        // s is a stack of operands to process when you encounter operators
        // in the postfix expression tree. Operands consist of three WORDs,
        // formatted as the LHS/RHS of the keys in ID.
        stack<int> s;
        vector<WORD> instr;
        WORD numinstr = 0;
        vector<WORD> x;

        // process the expression tree
        for (int i=0; i<(int)tree.size();) {
                x.clear();

                if (tree[i]==SNUMBER) {
                        WORD hash = hash_range(&tree[i], tree[i + 1]);
                        x.push_back(hash);
                        x.push_back(SNUMBER);
                        x.insert(x.end(),&tree[i],&tree[i]+tree[i+1]);
                        int sign = SGN(x.back());
                        x.back() *= sign;

                        std::pair<csemap::iterator, bool> suc = ID.insert(csemap::value_type(x, i + 1));
                        s.push(0);
                        s.push(sign * suc.first->second);
                        s.push(SNUMBER);
                        s.push(hash);
                        i+=tree[i+1];
                }
                else if (tree[i]==SYMBOL) {
                        WORD hash = hash_range(&tree[i], tree[i + 1]);
                        if (tree[i+3]>1) {
                                x.push_back(hash);
                                x.push_back(SYMBOL);
                                x.push_back(tree[i+2]+1); // variable (1-indexed)
                                x.push_back(tree[i+3]);   // power

                                csemap::iterator it = ID.find(x);
                                if (do_CSE && it != ID.end()) {
                                        // already-seen power of a symbol
                                        s.push(1);
                                        s.push(it->second);
                                        s.push(EXTRASYMBOL);
                                } else {
                                        //MesPrint("strange");
                                        if (numinstr == MAXPOSITIVE) {
                                                MesPrint((char *)"ERROR: too many temporary variables needed in optimization");
                                                Terminate(-1);
                                        }

                                        // new power greater than 1 of a symbol
                                        instr.push_back(numinstr);         // expr.nr
                                        instr.push_back(OPER_MUL);         // operator
                                        instr.push_back(9+OPTHEAD);  // length total
                                        instr.push_back(8);                  // length operand
                                        instr.push_back(SYMBOL);         // SYMBOL
                                        instr.push_back(4);                  // length symbol
                                        instr.push_back(tree[i+2]);  // variable
                                        instr.push_back(tree[i+3]);  // power
                                        instr.push_back(1);                  // numerator
                                        instr.push_back(1);                  // denominator
                                        instr.push_back(3);                  // length coeff
                                        instr.push_back(0);                  // trailing 0

                                        ID[x] = ++numinstr;
                                        s.push(1);
                                        s.push(numinstr);
                                        s.push(EXTRASYMBOL);
                                }
                        }
                        else {
                                // power of 1
                                s.push(tree[i+3]);         // power
                                s.push(tree[i+2]+1); // variable (1-indexed)
                                s.push(SYMBOL);
                        }

                        s.push(hash); // push hash
                        i+=tree[i+1];
                }
                else { // tree[i]==OPERATOR
                        int oper = tree[i];
                        i++;

                        x.push_back(0); // placeholder for hash
                        x.push_back(oper);
                        vector<WORD> hash;
                        hash.push_back(oper);

                        // get two operands from the stack
                        for (int operand=0; operand<2; operand++) {
                                hash.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                        }

                        x[0] = hash_range(&hash[0], 3);

                        // get rid of multiplications by +/-1
                        if (oper==OPER_MUL) {
                                bool do_continue = false;

                                for (int operand=0; operand<2; operand++) {
                                        int idx_oper1 = operand==0 ? 2 : 5;
                                        int idx_oper2 = operand==0 ? 5 : 2;

                                        // check whether operand 1 equals +/-1
                                        if (x[idx_oper1]==SNUMBER) {

                                                int idx = ABS(x[idx_oper1+1])-1;
                                                if (tree[idx+2]==1 && tree[idx+3]==1 && ABS(tree[idx+4])==3) {
                                                        // push +/- other operand and continue
                                                        s.push(x[idx_oper2+2]);
                                                        s.push(x[idx_oper2+1]*SGN(x[idx_oper1+1]));
                                                        s.push(x[idx_oper2]);
                                                        s.push(hash[1 + (operand + 1) % 2]);
                                                        do_continue = true;
                                                        break;
                                                }
                                        }
                                }

                                if (do_continue) continue;
                        }

                        // check whether this subexpression has been seen before
                        // if not, generate instruction to define it
                        csemap::iterator it = ID.find(x);
                        if (!do_CSE || it == ID.end()) {
                                if (numinstr == MAXPOSITIVE) {
                                        MesPrint((char *)"ERROR: too many temporary variables needed in optimization");
                                        Terminate(-1);
                                }

                                instr.push_back(numinstr); // expr.nr.
                                instr.push_back(x[1]);           // operator
                                instr.push_back(3);            // length

                                ID[x] = ++numinstr;

                                int lenidx = instr.size()-1;

                                for (int j=0; j<2; j++)
                                        if (x[3*j+2]==SYMBOL || x[3*j+2]==EXTRASYMBOL) {
                                                instr.push_back(8);                                            // length total
                                                instr.push_back(x[3*j+2]);                                   // (extra)symbol
                                                instr.push_back(4);                                            // length (extra)symbol
                                                instr.push_back(ABS(x[3*j+3])-1);                   // variable (0-indexed)
                                                instr.push_back(x[3*j+4]);                                   // power
                                                instr.push_back(1);                                            // numerator
                                                instr.push_back(1);                                            // denominator
                                                instr.push_back(3*SGN(x[3*j+3]));                   // length coeff
                                                instr[lenidx] += 8;
                                        }
                                        else { // x[3*j+1]==SNUMBER
                                                int t = ABS(x[3*j+3])-1;
                                                instr.push_back(tree[t+1]-1);                                                           // length number
                                                instr.insert(instr.end(), &tree[t+2], &tree[t]+tree[t+1]); // digits
                                                instr.back() *= SGN(instr.back()) * SGN(x[3*j+3]);
                                                instr[lenidx] += tree[t+1]-1;
                                        }

                                instr.push_back(0); // trailing 0
                                instr[lenidx]++;
                        }

                        // push new expression on the stack
                        s.push(1);
                        s.push(ID[x]);
                        s.push(EXTRASYMBOL);
                        s.push(x[0]); // push hash
                }
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: generate_instructions(cse=%d) : numoper=%d",
                                                thetime_str().c_str(), do_CSE?1:0, count_operators(instr));
#endif

        return instr;
}

/*
          #] generate_instructions : 
          #[ count_operators_cse :
*/


/**
* Count number of operators in a binary tree, while removing CSEs on the fly.
* The instruction set is not created, which makes this method slightly faster.
*
* A hash is created on the fly and is passed through the stack.
* TODO: find better hash functions
*/
int count_operators_cse (const vector<WORD> &tree) {
        //MesPrint ("*** [%s] Starting CSEE", thetime_str().c_str());

        typedef unordered_map<vector<WORD>, int, CSEHash, CSEEq> csemap;
        csemap ID;

        // reserve lots of space, to prevent later rehashes
        // TODO: what if this is too large? make a parameter?
        ID.rehash(mcts_expr_score * 2);

        // s is a stack of operands to process when you encounter operators
        // in the postfix expression tree. Operands consist of three WORDs,
        // formatted as the LHS/RHS of the keys in ID.
        stack<int> s;
        int numinstr = 0, numcommas = 0;
        vector<WORD> x;

        // process the expression tree
        for (int i=0; i<(int)tree.size();) {
                x.clear();

                if (tree[i]==SNUMBER) {
                        WORD hash = hash_range(&tree[i], tree[i + 1]);
                        x.push_back(hash);
                        x.push_back(SNUMBER);
                        x.insert(x.end(),&tree[i],&tree[i]+tree[i+1]);
                        int sign = SGN(x.back());
                        x.back() *= sign;

                        std::pair<csemap::iterator, bool> suc = ID.insert(csemap::value_type(x, i + 1));
                        s.push(0);
                        s.push(sign * suc.first->second);
                        s.push(SNUMBER);
                        s.push(hash);
                        i+=tree[i+1];
                }
                else if (tree[i]==SYMBOL) {
                        WORD hash = hash_range(&tree[i], tree[i + 1]);
                        if (tree[i+3]>1) {
                                x.push_back(hash);
                                x.push_back(SYMBOL);
                                x.push_back(tree[i+2]+1); // variable (1-indexed)
                                x.push_back(tree[i+3]);   // power

                                csemap::iterator it = ID.find(x);
                                if (it != ID.end()) {
                                        // already-seen power of a symbol
                                        s.push(1);
                                        s.push(it->second);
                                        s.push(EXTRASYMBOL);
                                } else {
                                        if (tree[i + 3] == 2)
                                                numinstr++;
                                        else
                                                numinstr += (int)floor(log(tree[i+3])/log(2.0)) + popcount(tree[i+3]) - 1;

                                        ID[x] = numinstr;
                                        s.push(1);
                                        s.push(numinstr);
                                        s.push(EXTRASYMBOL);
                                }
                        }
                        else {
                                // power of 1
                                s.push(tree[i+3]);         // power
                                s.push(tree[i+2]+1); // variable (1-indexed)
                                s.push(SYMBOL);
                        }

                        s.push(hash); // push hash

                        i+=tree[i+1];
                }
                else { // tree[i]==OPERATOR
                        int oper = tree[i];
                        i++;

                        x.push_back(0); // placeholder for hash
                        x.push_back(oper);
                        vector<WORD> hash;
                        hash.push_back(oper);

                        // get two operands from the stack
                        for (int operand=0; operand<2; operand++) {
                                hash.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                                x.push_back(s.top()); s.pop();
                        }


                        x[0] = hash_range(&hash[0], 3);

                        // get rid of multiplications by +/-1
                        if (oper==OPER_MUL) {
                                bool do_continue = false;

                                for (int operand=0; operand<2; operand++) {
                                        int idx_oper1 = operand==0 ? 2 : 5;
                                        int idx_oper2 = operand==0 ? 5 : 2;

                                        // check whether operand 1 equals +/-1
                                        if (x[idx_oper1]==SNUMBER) {

                                                int idx = ABS(x[idx_oper1+1])-1;
                                                if (tree[idx+2]==1 && tree[idx+3]==1 && ABS(tree[idx+4])==3) {
                                                        // push +/- other operand and continue
                                                        s.push(x[idx_oper2+2]);
                                                        s.push(x[idx_oper2+1]*SGN(x[idx_oper1+1]));
                                                        s.push(x[idx_oper2]);
                                                        s.push(hash[1 + (operand + 1) % 2]);
                                                        do_continue = true;
                                                        break;
                                                }
                                        }
                                }

                                if (do_continue) continue;
                        }

                        // check whether this subexpression has been seen before
                        // if not, generate instruction to define it
                        csemap::iterator it = ID.find(x);
                        if (it == ID.end()) {
                                if (numinstr == MAXPOSITIVE) {
                                        MesPrint((char *)"ERROR: too many temporary variables needed in optimization");
                                        Terminate(-1);
                                }

                                if (oper == OPER_COMMA) numcommas++;
                                ID[x] = ++numinstr;
                                s.push(1);
                                s.push(numinstr);
                                s.push(EXTRASYMBOL);
                        } else {
                                // push new expression on the stack
                                s.push(1);
                                s.push(it->second);
                                s.push(EXTRASYMBOL);
                        }

                        s.push(x[0]); // push hash
                }
        }

        //MesPrint ("*** [%s] Stopping CSEE", thetime_str().c_str());
        return numinstr - numcommas;
}

/*
          #] count_operators_cse : 
          #[ count_operators_cse_topdown :
*/

typedef struct node {
        const WORD* data;
        struct node* l;
        struct node* r; // TODO: add l,r to data?
        WORD sign; // TODO: use data for this?
        UWORD hash;

        node() : l(NULL), r(NULL), sign(1), hash(0) {};
        node(const WORD* data) : data(data), l(NULL), r(NULL), sign(1), hash(0) {};

        // a minus sign in the tree should only count as a different entry if
        // it is a compound expression: a = -a, but T+-V != T+V
        int cmp(const struct node* rhs) const {
                if (this == rhs) return 0;

                if (data[0] != rhs->data[0]) return data[0] < rhs->data[0] ? -1 : 1;
                int mod = data[0] == SNUMBER ? -1 : 0; // don't check sign, for numbers
                if (data[0] == SYMBOL || data[0] == SNUMBER) {
                        for (int i = 0; i < data[1] + mod; i++) {
                                if (data[i] != rhs->data[i]) return data[i] < rhs->data[i] ? -1 : 1;
                                }
                } else {
                        int lv = l->cmp(rhs->l);
                        if (lv != 0) return lv;
                        int rv = r->cmp(rhs->r);
                        if (rv != 0) return rv;

                        // TODO: only for ADD operation
                        if (l->sign != rhs->l->sign) return l->sign < rhs->l->sign ? -1 : 1;
                        if (r->sign != rhs->r->sign) return r->sign < rhs->r->sign ? -1 : 1;
                }

                return 0;
        }

                // less than operator
        bool operator() (const struct node* lhs, const struct node* rhs) const
        {
                return lhs->cmp(rhs) < 0;
        }

        void calcHash() {
                int mod = data[0] == SNUMBER ? -1 : 0; // don't check sign, for numbers
                if (data[0] == SYMBOL || data[0] == SNUMBER) {
                        hash = hash_range(data, data[1] + mod);
                } else {
                        if (l->hash == 0) l->calcHash();
                        if (r->hash == 0) r->calcHash();

                        // signs only matter for compound expressions
                        size_t newr[] = {(size_t)data[0], l->hash, (size_t)l->sign, r->hash, (size_t)r->sign};
                        hash = hash_range(newr, 5);
                }
        }
} NODE;

struct NodeHash {
        size_t operator()(const NODE* n) const {
                return n->hash; // already computed
        }
};

struct NodeEq {
        bool operator()(const NODE* lhs, const NODE* rhs) const {
                return lhs->cmp(rhs) == 0;
        }
};

NODE* buildTree(vector<WORD> &tree) {
        //MesPrint ("*** [%s] Starting CSEE topdown", thetime_str().c_str());

        // allocate spaces for the tree, cannot be more nodes than tree size
        NODE* ar = (NODE*)Malloc1(tree.size() * sizeof(NODE), "CSE tree");
        NODE* c = 0;
        unsigned int curIndex = 0;

        stack<NODE*> st;
        for (int i=0; i<(int)tree.size();) {
                c = ar + curIndex;
                new (c) NODE(&tree[i]); // placement new
                curIndex++;

                if (tree[i]==SYMBOL || tree[i] == SNUMBER) {
                        // extract the sign to a new class member
                        if (tree[i] == SNUMBER) {
                                c->sign = SGN(tree[i + tree[i + 1] -1]);
                        }

                        c->calcHash();
                        st.push(c);
                        i+=tree[i+1];
                } else {
                        c->r = st.top(); st.pop();
                        c->l = st.top(); st.pop();

                        // filter *1 and *-1
                        // TODO: also multiply if there are two numbers?
                        if (c->data[0] == OPER_MUL) {
                                NODE* ch[] = {c->r, c->l};
                                for (int j = 0; j < 2; j++)
                                        if (ch[j]->data[0] == SNUMBER && ch[j]->data[1] == 5 && ch[j]->data[2]==1 && ch[j]->data[3]==1) {
                                                        ch[(j+1)%2]->sign *= ch[j]->sign; // transfer sign
                                                        c = ch[(j+1)%2];
                                                        break;
                                        }
                        }

                        c->calcHash();
                        st.push(c);
                        i++;
                }
        }

        // TODO: reallocate to smaller size? Could save memory
        //MesPrint("Memory difference: %d vs %d", curIndex, tree.size());

        // we want to make the root of the tree the first element
        // so that we can easily free the array.
        // we swap the first element with the root
        // we need to change the pointer in the operator node that has this element as a child
        // TODO: check performance
        for (unsigned int i = 0; i < curIndex; i++) {
                if (ar[i].l == ar) ar[i].l = st.top();
                if (ar[i].r == ar) ar[i].r = st.top();
        }

        swap(ar[0], *st.top());        
        return ar;
}

int count_operators_cse_topdown (vector<WORD> &tree) {
        typedef unordered_set<NODE*, NodeHash, NodeEq> nodeset;
        nodeset ID;

        // reserve lots of space, to prevent later rehashes
        // TODO: what if this is too large? make a parameter?
        ID.rehash(mcts_expr_score * 2);

        int numinstr = 0;

        NODE* root = buildTree(tree);

        stack<NODE*> stack;
        stack.push(root);
        while (!stack.empty())
        {
                NODE* c = stack.top();
                stack.pop();

                if (c->data[0] == SYMBOL) {
                        if (c->data[3] > 1) {
                                std::pair<nodeset::iterator, bool> suc = ID.insert(c);
                                if (suc.second) { // new
                                        if (c->data[3] == 2)
                                                numinstr++;
                                        else
                                                numinstr += (int)floor(log(c->data[3])/log(2.0)) + popcount(c->data[3]) - 1;
                                }
                        }
                } else {
                        if (c->data[0] != SNUMBER) {
                                // operator
                                std::pair<nodeset::iterator, bool> suc = ID.insert(c);
                                if (suc.second) {
                                        stack.push(c->r);
                                        stack.push(c->l);

                                        // ignore OPER_COMMA
                                        if (c->data[0] == OPER_MUL || c->data[0] == OPER_ADD)
                                                numinstr++;
                                }
                        }
                }
        }

        //MesPrint ("*** [%s] Stopping CSEE", thetime_str().c_str());
        M_free(root, "CSE tree");

        return numinstr;
}

/*
          #] count_operators_cse_topdown : 
          #[ simulated_annealing :
*/
vector<WORD> simulated_annealing() {
        float minT = AO.Optimize.saMinT.fval;
        float maxT = AO.Optimize.saMaxT.fval;
        float T = maxT;
        float coolrate = pow(minT / maxT, 1 / (float)AO.Optimize.saIter);

        GETIDENTITY;

        // create a valid state where FACTORSYMBOL/SEPARATESYMBOL remains first
        vector<WORD> state = occurrence_order(optimize_expr, false);
        int startindex = 0;
        if ((state.size() > 0 && state[0] == SEPARATESYMBOL) || (state.size() > 1 && state[1] == FACTORSYMBOL)) startindex++;
        if (state.size() > 1 && state[1] == FACTORSYMBOL) startindex++;

        my_random_shuffle(BHEAD state.begin() + startindex, state.end()); // start from random scheme

        vector<WORD> tree = Horner_tree(optimize_expr, state);
        int curscore = count_operators_cse_topdown(tree);

        // clean poly_vars, that are allocated by Horner_tree
        AN.poly_num_vars = 0;
        M_free(AN.poly_vars,"poly_vars");

        std::vector<WORD> best = state; // best state
        int bestscore = curscore;

        if (startindex == (int)state.size() || (int)state.size() - startindex < 2) {
                return best;
        }
        
        for (int o = 0; o < AO.Optimize.saIter; o++) {
                int inda = iranf(BHEAD state.size() - startindex) + startindex;
                int indb = iranf(BHEAD state.size() - startindex) + startindex;

                if (inda == indb) {
                        continue;
                }

                swap(state[inda], state[indb]); // swap works best for Horner

                vector<WORD> tree = Horner_tree(optimize_expr, state);
                int newscore = count_operators_cse_topdown(tree);

                // clean poly_vars, that are allocated by Horner_tree
                AN.poly_num_vars = 0;
                M_free(AN.poly_vars,"poly_vars");

                if (newscore <= curscore || 2.0 * wranf(BHEAD0) / (float)(UWORD)(-1) < exp((curscore - newscore) / T)) {
                        curscore = newscore;

                        if (curscore < bestscore) {
                                bestscore = curscore;
                                best = state;
                        }
                } else {
                        swap(state[inda], state[indb]);
                }

#ifdef DEBUG_SA
        MesPrint("Score at step %d: %d", o, curscore);
#endif
                T *= coolrate;
        }

#ifdef DEBUG_SA
        MesPrint("Simulated annealing score: %d", bestscore);
#endif

        return best;
}

/*
          #] simulated_annealing : 
          #[ printpstree :
*/

/*
// print MCTS tree with LaTeX/pstricks (for analysis)
void printpstree_rec (tree_node x, string pre="") {

        if (x.num_visits==1) {
                MesPrint("%s\\TR{%d}",pre.c_str(),x.var);
        }
        else {
                MesPrint("%s\\pstree%s{\\TR{%d}}{",pre.c_str(),
                                                 pre=="  "?"[nodesep=0, levelsep=40]":"",
                                                 x.var);
                for (int i=0; i<(int)x.childs.size(); i++)
                        if (x.childs[i].num_visits>0)
                                printpstree_rec(x.childs[i], pre+"        ");
                MesPrint("%s}",pre.c_str());
        }
}

void printpstree () {
        // draw tree with pstricks
        MesPrint ("\\documentclass{article}");
        MesPrint ("\\usepackage{pstricks,pst-node,pst-tree,graphicx}");
        MesPrint ("\\begin{document}");
        MesPrint ("\\scalebox{0.02}{");
        printpstree_rec(mcts_root,"  ");
        MesPrint ("}");
        MesPrint ("\\end{document}");
}
*/

/*
          #] printpstree : 
          #[ find_Horner_MCTS_expand_tree :
*/

/**  Expand MCTS tree
 *
 *         Description
 *         ===========
 *         This method does one MCTS step: it selects the most-promising
 *         node, expands it, randomly completes the Horner scheme and
 *         backpropagates the results.
 *
 *         Selection is done according to the UCT formula:
 *
 *                UCT(i) = <x(i)> + C * sqrt(2*log(N)/n(i)),
 *
 *         where <x(i)> is the average result of child i, n(i) is the number
 *         of time child i is visited, N=SUM(n(i)) and C is a constant to be
 *         determined experimentally (can be set via mctsconstant).
 *
 *         A "virtual loss" is added once a node is selected. This is
 *         relevant to avoid duplicate work in the parallel version.
 *
 *         Notes
 *         =====
 *         - The method is called from "find_Horner_MCTS" in Form and from
 *           "RunThread" via "find_Horner_MCTS_expand_tree_threaded" in
 *           TForm.
 *         - The code is divided into three functions: "next_MCTS_scheme",
 *           "try_MCTS_scheme" and "update_MCTS_scheme". In this way, the
 *           source code is shared with ParForm; "try_MCTS_scheme" is
 *           assumed to run on workers, while the others are assumed to run
 *           on the master.
 */

/*
                 #[ next_MCTS_scheme :
*/

// find a Horner scheme to be used for the next simulation
inline static void next_MCTS_scheme (PHEAD vector<WORD> *porder, vector<WORD> *pscheme, vector<tree_node *> *ppath) {

        vector<WORD> &order = *porder;
        vector<WORD> &schemev = *pscheme;
        vector<tree_node *> &path = *ppath;
        int depth = 0, nchild0;
        float slide_down_factor = 1.0;

        order.clear();
        path.clear();

        // MCTS step I: select
        tree_node *select = &mcts_root;
        path.push_back(select);
        nchild0 = select->childs.size();
        while (select->childs.size() > 0) {
                // add virtual loss
                select->num_visits++;
                select->sum_results+=mcts_expr_score;

//-------------------------------------------------------------------
                                switch ( AO.Optimize.mctsdecaymode ) {
                                        case 1:  // Based on http://arxiv.org/abs/arXiv:1312.0841
                                                slide_down_factor = 1.0-(1.0*AT.optimtimes)/(1.0*AO.Optimize.mctsnumexpand);
                                                break;
                                        case 2:  // This gives a bit more cleanup time at the end.
                                                if ( 2*AT.optimtimes < AO.Optimize.mctsnumexpand ) {
                                                        slide_down_factor = 1.0*(AO.Optimize.mctsnumexpand-2*AT.optimtimes);
                                                        slide_down_factor /= 1.0*AO.Optimize.mctsnumexpand;
                                                }
                                                else {
                                                        slide_down_factor = 0.0001;
                                                }
                                                break;
                                        case 3:  // depth dependent factor combined with case 1
                                                float dd = 1.0-(1.0*depth)/(1.0*nchild0);
                                                slide_down_factor = 1.0-(1.0*AT.optimtimes)/(1.0*AO.Optimize.mctsnumexpand);
                                                if ( dd <= 0.000001 ) slide_down_factor = 1.0;
                                                else slide_down_factor /= dd;
                                                if ( slide_down_factor > 1.0 ) slide_down_factor = 1.0;
                                                break;
                                }
//-------------------------------------------------------------------

#ifdef DEBUG_MCTS
                MesPrint("select %d",select->var);
#endif

                // find most-promising node
                double best=0;
                tree_node *next=NULL;
                for (vector<tree_node>::iterator p=select->childs.begin(); p<select->childs.end(); p++) {
                        double score;
                        if (p->num_visits >= 1) {

                                // there are results calculated, so select with the UCT formula
                                score = mcts_expr_score / (p->sum_results/p->num_visits) +
//-------------------------------------------------------------------------
                                        slide_down_factor *
//-------------------------------------------------------------------------
                                        2 * AO.Optimize.mctsconstant.fval * sqrt(2*log(select->num_visits) / p->num_visits);

#ifdef DEBUG_MCTS
                                printf("%d: %.2lf [x=%.2lf n=%d fin=%i]\n",p->var,score,mcts_expr_score / (p->sum_results/p->num_visits),
                                                         p->num_visits,p->finished?1:0);
                                fflush(stdout);
#endif
                        }
                        else {
                // no results yet, so select this node by setting score=infinite
                                score = 1e100;

#ifdef DEBUG_MCTS
                                printf("%d: inf\n",p->var); fflush(stdout);
#endif
                        }

                        // update best candidate
                        if (!p->finished && score>best) {
                                best=score;
                                next=&*p;
                        }
                }

                // if no node is found, this node is finished
                if (next==NULL) {
                        select->finished=true;
                        break;
                }

                // traverse down the tree
                select = next;
                path.push_back(select);
                order.push_back(select->var);
                depth++;
        }

        // MCTS step II: expand

#ifdef DEBUG_MCTS
        MesPrint("expand %d",select->var);
#endif

        // variables used so far
        set<WORD> var_used;

        for (int i=0; i<(int)order.size(); i++)
                var_used.insert(ABS(order[i])-1);

        // if this a new node, create node and add children
        if (!select->finished && select->childs.size()==0) {
                tree_node new_node(select->var);
                int sign = (order.size() == 0) ? 1 : SGN(order.back());
                for (int i=0; i<(int)mcts_vars.size(); i++)
                        if (!var_used.count(mcts_vars[i])) {
                                new_node.childs.push_back(tree_node(sign*(mcts_vars[i]+1)));
                                if (AO.Optimize.hornerdirection==O_FORWARDANDBACKWARD)
                                        new_node.childs.push_back(tree_node(-sign*(mcts_vars[i]+1)));
                        }
                my_random_shuffle(BHEAD new_node.childs.begin(), new_node.childs.end());

                // here locking is necessary, since operator=(tree_node) is a
                // non-atomic operation (using pointers makes this lock obsolete)
                LOCK(optimize_lock);
                *select = new_node;
                UNLOCK(optimize_lock);
        }
        // set finished if necessary
        if (select->childs.size()==0)
                select->finished = true;

        // add virtual loss of number of operators in original expression
        select->num_visits++;
        select->sum_results+=mcts_expr_score;

        // MCTS step III: simulation

        // create complete Horner scheme
        deque<WORD> scheme;

        for (int i=0; i<(int)mcts_vars.size(); i++)
                if (!var_used.count(mcts_vars[i]))
                        scheme.push_back(mcts_vars[i]);
        my_random_shuffle(BHEAD scheme.begin(), scheme.end());

        for (int i=(int)order.size()-1; i>=0; i--) {
                if (order[i] > 0)
                        scheme.push_front(order[i]-1);
                else
                        scheme.push_back(-order[i]-1);
        }

        // add FACTORSYMBOL/SEPARATESYMBOL is necessary
        if (mcts_factorized)
                scheme.push_front(FACTORSYMBOL);
        if (mcts_separated)
                scheme.push_front(SEPARATESYMBOL);

        // Horner scheme as a vector
        schemev = vector<WORD>(scheme.begin(), scheme.end());

}

/*
                 #] next_MCTS_scheme : 
                 #[ try_MCTS_scheme :
*/

// count the number of operators in the given Horner scheme
inline static void try_MCTS_scheme (PHEAD const vector<WORD> &scheme, int *pnum_oper) {

        // do Horner, CSE and count the number of operators
        vector<WORD> tree = Horner_tree(optimize_expr, scheme);
        //vector<WORD> instr = generate_instructions(tree, true);
        //int num_oper = count_operators(instr);
        //int num_oper2 = count_operators_cse(tree);
        //int num_oper2 = count_operators_cse_topdown(tree);
        //MesPrint("%d %d", num_oper, num_oper2);

        int num_oper = count_operators_cse_topdown(tree);

        // clean poly_vars, that is allocated by Horner_tree
        AN.poly_num_vars = 0;
        M_free(AN.poly_vars,"poly_vars");

        *pnum_oper = num_oper;

}

/*
                 #] try_MCTS_scheme : 
                 #[ update_MCTS_scheme :
*/

// update the best score and the search tree
inline static void update_MCTS_scheme (int num_oper, const vector<WORD> &scheme, vector<tree_node *> *ppath) {

        vector<tree_node *> &path = *ppath;

        // update the (global) list of best Horner scheme
        if ((int)mcts_best_schemes.size() < AO.Optimize.mctsnumkeep ||
                        (--mcts_best_schemes.end())->first > num_oper) {
                // here locking is necessary, for otherwise best_schemes may
                // become corrupted; lock can be prevented if each thread keeps
                // track of it's own list and those lists are merged in the end,
                // but this seems not useful to implement
                LOCK(optimize_lock);
                mcts_best_schemes.insert(make_pair(num_oper,scheme));
                if ((int)mcts_best_schemes.size() > AO.Optimize.mctsnumkeep)
                        mcts_best_schemes.erase(--mcts_best_schemes.end());
                UNLOCK(optimize_lock);
        }

        // MCTS step IV: backpropagate

        // add number of operator and subtract mcts_expr_score, which
        // behaves as a "virtual loss"
        for (vector<tree_node *>::iterator p=path.begin(); p<path.end(); p++)
                (*p)->sum_results += num_oper - mcts_expr_score;

}

/*
                 #] update_MCTS_scheme : 
*/

void find_Horner_MCTS_expand_tree () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: find_Horner_MCTS_expand_tree", thetime_str().c_str());
#endif

        GETIDENTITY;

        // the order for the Horner scheme up to the selected node, with signs
        // indicating forward or backward.
        vector<WORD> order;

        // complete Horner scheme
        vector<WORD> scheme;

        // path to the selected node
        vector<tree_node *> path;

        // the number of operations obtained by the simulation
        int num_oper;

        next_MCTS_scheme(BHEAD &order, &scheme, &path);
        try_MCTS_scheme(BHEAD scheme, &num_oper);
#ifdef DEBUG_MCTS
        // Actually "order" is needed only for this debug output.
        MesPrint ("{%a} -> {%a} -> %d", order.size(), &order[0], scheme.size(), &scheme[0], num_oper);
#endif
        update_MCTS_scheme(num_oper, scheme, &path);

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: find_Horner_MCTS_expand_tree(%a-> %d)",
                                                thetime_str().c_str(), scheme.size(), &scheme[0], num_oper);
#endif

}

/*
          #] find_Horner_MCTS_expand_tree : 
          #[ PF_find_Horner_MCTS_expand_tree :
*/
#ifdef WITHMPI

// To remember which task is assigned to each slave. A task is represented by
// a pair of a Horner scheme and the selected path for the scheme.
// The index range is from 1 to PF.numtasks-1.
vector<pair<vector<WORD>, vector<tree_node *> > > PF_opt_MCTS_tasks;

// Initialization.
void PF_find_Horner_MCTS_expand_tree_master_init () {

        PF_opt_MCTS_tasks.resize(PF.numtasks);
        for (int i = 1; i < PF.numtasks; i++) {
                pair<vector<WORD>, vector<tree_node *> > &p = PF_opt_MCTS_tasks[i];
                p.first.clear();
                p.second.clear();
        }

}

// Wait for an idle slave and return the process number.
int PF_find_Horner_MCTS_expand_tree_master_next () {

        // Find an idle slave.
        int next;
        PF_Receive(PF_ANY_SOURCE, PF_OPT_MCTS_MSGTAG, &next, NULL);

        // Check if the slave had a task.
        pair<vector<WORD>, vector<tree_node *> > &p = PF_opt_MCTS_tasks[next];
        if (!p.first.empty()) {
                // If so, update the result.
                int num_oper;
                PF_Unpack(&num_oper, 1, PF_INT);
                update_MCTS_scheme(num_oper, p.first, &p.second);

                // Clear the task.
                p.first.clear();
                p.second.clear();
        }

        return next;

}

// The main function on the master.
void PF_find_Horner_MCTS_expand_tree_master () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_find_Horner_MCTS_expand_tree_master", thetime_str().c_str());
#endif

        vector<WORD> order;
        vector<WORD> scheme;
        vector<tree_node *> path;

        next_MCTS_scheme(BHEAD &order, &scheme, &path);

        // Find an idle slave.
        int next = PF_find_Horner_MCTS_expand_tree_master_next();

        // Send a new task to the slave.
        PF_PrepareLongSinglePack();
        int len = scheme.size();
        PF_LongSinglePack(&len, 1, PF_INT);
        PF_LongSinglePack(&scheme[0], len, PF_WORD);
        PF_LongSingleSend(next, PF_OPT_MCTS_MSGTAG);

        // Remember the task.
        pair<vector<WORD>, vector<tree_node *> > &p = PF_opt_MCTS_tasks[next];
        p.first = scheme;
        p.second = path;

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_find_Horner_MCTS_expand_tree_master", thetime_str().c_str());
#endif

}

// Wait for all the slaves to finish their tasks.
void PF_find_Horner_MCTS_expand_tree_master_wait () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_find_Horner_MCTS_expand_tree_master_wait", thetime_str().c_str());
#endif

        // Wait for all the slaves.
        for (int i = 1; i < PF.numtasks; i++) {
                int next = PF_find_Horner_MCTS_expand_tree_master_next();
                // Send a null task.
                PF_PrepareLongSinglePack();
                int len = 0;
                PF_LongSinglePack(&len, 1, PF_INT);
                PF_LongSingleSend(next, PF_OPT_MCTS_MSGTAG);
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_find_Horner_MCTS_expand_tree_master_wait", thetime_str().c_str());
#endif

}

// The main function on the slaves.
void PF_find_Horner_MCTS_expand_tree_slave () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_find_Horner_MCTS_expand_tree_slave", thetime_str().c_str());
#endif

        vector<WORD> scheme;

        {
                // Send the first message to the master, which indicates I am idle.
                PF_PreparePack();
                int dummy = 0;
                PF_Pack(&dummy, 1, PF_INT);
                PF_Send(MASTER, PF_OPT_MCTS_MSGTAG);
        }

        for (;;) {
                // Get a task from the master.
                PF_LongSingleReceive(MASTER, PF_OPT_MCTS_MSGTAG, NULL, NULL);

                // Length of the task.
                int len;
                PF_LongSingleUnpack(&len, 1, PF_INT);

                // No task remains.
                if (len == 0) break;

                // Perform the given task.
                scheme.resize(len);
                PF_LongSingleUnpack(&scheme[0], len, PF_WORD);
                int num_oper;
                try_MCTS_scheme(scheme, &num_oper);

                // Send the result to the master.
                PF_PreparePack();
                PF_Pack(&num_oper, 1, PF_INT);
                PF_Send(MASTER, PF_OPT_MCTS_MSGTAG);
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_find_Horner_MCTS_expand_tree_slave", thetime_str().c_str());
#endif

}

#endif
/*
          #] PF_find_Horner_MCTS_expand_tree : 
          #[ find_Horner_MCTS :
*/

/**  Find best Horner schemes using MCTS
 *
 *         Description
 *         ===========
 *         The method governs the MCTS for the best Horner schemes. It does
 *         some pre-processing, calls "find_Horner_MCTS_expand_tree" a
 *         number of times and does some post-processing.
 */
//vector<vector<WORD> > find_Horner_MCTS () {
void find_Horner_MCTS () {

#ifdef WITHMPI
        if (PF.me != MASTER) {
                if (PF.numtasks <= 1) return;
                PF_find_Horner_MCTS_expand_tree_slave();
                return;
        }
#endif

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: find_Horner_MCTS", thetime_str().c_str());
#endif

        GETIDENTITY;

        LONG start_time = TimeWallClock(1);

        // initialize the used global variables
        mcts_expr_score = count_operators(optimize_expr);
        mcts_root = tree_node();

        // extract all symbols from the expression
        set<WORD> var_set;
        for (WORD *t=optimize_expr; *t!=0; t+=*t)
                if (t[1] == SYMBOL)
                        for (int i=3; i<t[2]; i+=2)
                                var_set.insert(t[i]);

        // check for factorized/separated expression and make sure that
        // FACTORSYMBOL/SEPARATESYMBOL isn't included in the MCTS
        mcts_factorized = var_set.count(FACTORSYMBOL);
        if (mcts_factorized) var_set.erase(FACTORSYMBOL);
        mcts_separated = var_set.count(SEPARATESYMBOL);
        if (mcts_separated) var_set.erase(SEPARATESYMBOL);

        mcts_vars = vector<WORD>(var_set.begin(), var_set.end());
        optimize_num_vars = (int)mcts_vars.size();
        // initialize MCTS tree root
        for (int i=0; i<(int)mcts_vars.size(); i++) {
                if (AO.Optimize.hornerdirection != O_BACKWARD)
                        mcts_root.childs.push_back(tree_node(+(mcts_vars[i]+1)));
                if (AO.Optimize.hornerdirection != O_FORWARD)
                        mcts_root.childs.push_back(tree_node(-(mcts_vars[i]+1)));
        }
        my_random_shuffle(BHEAD mcts_root.childs.begin(), mcts_root.childs.end());

#if defined(WITHMPI)
        PF_find_Horner_MCTS_expand_tree_master_init();
#endif

        // initialize a potential variable mctsconstant scheme.
        AT.optimtimes = 0;

        // call expand_tree until it is called "mctsnumexpand" times, the
        // time limit is reached or the tree is fully finished
        for (int times=0; times<AO.Optimize.mctsnumexpand && !mcts_root.finished &&
                                 (AO.Optimize.mctstimelimit==0 ||        (TimeWallClock(1)-start_time)/100 < AO.Optimize.mctstimelimit);
                         times++) {
                AT.optimtimes = times;
        // call expand_tree routine depending on threading mode
#if defined(WITHPTHREADS)
                if (AM.totalnumberofthreads > 1)
                        find_Horner_MCTS_expand_tree_threaded();
                else
#elif defined(WITHMPI)
                if (PF.numtasks > 1)
                        PF_find_Horner_MCTS_expand_tree_master();
                else
#endif
                        find_Horner_MCTS_expand_tree();
        }

        // if TForm or ParForm, wait for everyone to finish
#ifdef WITHPTHREADS
        MasterWaitAll();
#endif
#ifdef WITHMPI
        PF_find_Horner_MCTS_expand_tree_master_wait();
#endif
#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: find_Horner_MCTS", thetime_str().c_str());
#endif

}

/*
          #] find_Horner_MCTS : 
          #[ merge_operators :
*/

/**  Merge operators
 *
 *         Description
 *         ===========
 *         The input instructions form a binary DAG. This method merges
 *         expressions like
 *
 *                Z1 = a+b;
 *                Z2 = Z1+c;
 *
 *         into
 *
 *                Z2 = a+b+c;
 *
 *         An instruction is merged iff it only has one parent and the
 *         operator equals its parent's operator.
 *
 *         This still leaves some freedom: where should the coefficients end
 *         up in cases as:
 *
 *                Z1 = Z2 + x    <=>         Z1 = 2*Z2 + x
 *                Z2 = 2*x*y                         Z2 = x*y
 *
 *         Both are relevant, e.g. for CSE of the form "2*x" and "2*Z2". The
 *         flag "move_coeff" moves coefficients from LHS-like expressions to
 *         RHS-like expressions.
 *
 *         Furthermore, this method removes empty equation (Z1=0), that are
 *         introduced by some "optimize_greedy" substitutions.
 *
 *         Implementation details
 *         ======================
 *         Expressions are mostly traversed via a stack, so that parents are
 *         evaluated before their children.
 *
 *         With "move_coeff" set coefficients are moved, but this leads to
 *         some tricky cases, e.g.
 *
 *           Z1 = Z2 + x
 *           Z2 = 2*y
 *
 *         Here Z2 reduces to the trivial equation Z2=y, which should be
 *         eliminated. Here the array skip[i] comes in.
 *
 *         Furthermore in the case
 *
 *           Z1 = Z2 + x
 *           Z2 = 2*Z3
 *           Z3 = x*Z4
 *           Z4 = y*z
 *
 *         after substituting Z1 = 2*Z3 + x, the parent expression for Z4
 *         becomes Z3 instead of Z2. This is where renum_par[i] comes in.
 *
 *         Finally, once a coefficient has been moved, skip_coeff[i] is set
 *         and this coefficient is copied into the new expression anymore.
 */
vector<WORD> merge_operators (const vector<WORD> &all_instr, bool move_coeff) {

#ifdef DEBUG_MORE
        MesPrint ("*** [%s, w=%w] CALL: merge_operators", thetime_str().c_str());
#endif

        // get starting positions of instructions
        vector<const WORD *> instr;
        const WORD *tbegin = &*all_instr.begin();
        const WORD *tend = tbegin+all_instr.size();
        // copy all instructions to temp space. There will be n of them in instr.
        for (const WORD *t=tbegin; t!=tend; t+=*(t+2)) {
                instr.push_back(t);
        }
        int n = instr.size();

        // find parents and number of parents of instructions
        vector<int> par(n), numpar(n,0);
        for (int i=0; i<n; i++) par[i]=i;

        for (int i=0; i<n; i++) {
                for (const WORD *t=instr[i]+OPTHEAD; *t!=0; t+=*t) {
                        if ( *(t+1)==EXTRASYMBOL && *t!=1+ABS(*(t+*t-1)) ) {
                                // extra symbol t[3] is referred to in instr i.
                                par[*(t+3)]=i;
                                numpar[*(t+3)]++;

                                // if coefficient isn't +/-1 or power > 1, increase numpar,
                                // so that this is not merged
                                if (*(t+*t-3)!=1 || *(t+*t-2)!=1 || ABS(*(t+*t-1))!=3) numpar[*(t+3)]++;
                                if (*(t+4)>1) numpar[*(t+3)]++;
                        }
                }
        }
        // determine which instructions to merge
        stack<int> s;
        s.push(n-1);
        vector<bool> vis(n,false);

        while (!s.empty()) {

                int i=s.top(); s.pop();
                if (vis[i]) continue;
                vis[i]=true;

                for (const WORD *t=instr[i]+OPTHEAD; *t!=0; t+=*t)
                        if ( *(t+1)==EXTRASYMBOL && *t!=1+ABS(*(t+*t-1)) )
                                s.push(*(t+3));

                // condition: one parent and equal operator as parent
                if (numpar[i]==1 && *(instr[i]+1)==*(instr[par[i]]+1))
                        par[i] = par[par[i]]; // The expr into which we subst par[i] to get i
                else
                        par[i] = i;
        }

        // merge instructions into new instructions
        vector<WORD> newinstr;

        // stack of new expressions, all 0-indexed
        stack<WORD> new_expr;
        new_expr.push(n-1);
        vis = vector<bool>(n,false);

        // skip empty equations (might be introduced by greedy optimizations)
        vector<bool> skip(n,false), skipcoeff(n,false);
        for (int i=0; i<n; i++)
                if (*(instr[i]+OPTHEAD) == 0) skip[i]=skipcoeff[i]=true;

        // for renumbering merged parents
        vector<int> renum_par(n);
        for (int i=0; i<n; i++)
                renum_par[i]=i;

        while (!new_expr.empty()) {
                int x = new_expr.top(); new_expr.pop();
                if (vis[x]) continue;
                vis[x] = true;

                // find all instructions with parent=x and copy their arguments
                // into a new expression
                bool first_copy=true;
                int lenidx = newinstr.size()+2;

                // 1-indexed, since signs may occur
                stack<WORD> this_expr;
                this_expr.push(x+1);

                while (!this_expr.empty()) {
                        // pop from stack, determine expr.nr and sign
                        int i = this_expr.top(); this_expr.pop();
                        int sign = SGN(i);
                        i = ABS(i)-1;
                        for (const WORD *t=instr[i]+OPTHEAD; *t!=0; t+=*t) { // terms in i
                                // don't copy a term if:
                                // (1) skip=true, since then it's merged into the parent
                                // (2) extrasymbol with parent=x, because its children should be copied
                                // (3) coefficient with skipcoeff=true, since it's already copied
                                bool copy = !skip[i];
                                if (*t!=1+ABS(*(t+*t-1)) && *(t+1)==EXTRASYMBOL) {
                                        if (par[*(t+3)] == x) {
                                                // parent of term refers to x. we push it with its sign if no skip is true
                                                // and the sign of the expr.
                                                this_expr.push(sign * (skip[i]||skipcoeff[i] ? 1 : SGN(*(t+*t-1))) * (1+*(t+3)));
                                                if (*(instr[i]+1) == OPER_MUL) sign=1;
                                                copy=false;
                                        }
                                        else {
                                                new_expr.push(*(t+3));
                                        }
                                }

                                if (*t == 1+ABS(*(t+*t-1)) && skipcoeff[i]) {
                                        copy=false;
                                }

                                if (copy) {
                                        // first term, so add header
                                        if (first_copy) {
                                                newinstr.push_back(renum_par[x]);  // expr.nr.
                                                newinstr.push_back(*(instr[x]+1)); // operator
                                                newinstr.push_back(3);                           // length   OPTHEAD?
                                                first_copy=false;
                                        }

                                        // copy term and adjust sign
                                        int thislenidx = newinstr.size();
                                        newinstr.insert(newinstr.end(), t, t+*t); // Put the whole term in newinstr
                                        newinstr.back() *= sign;
                                        if (*(instr[i]+1) == OPER_MUL) sign=1;
                                        newinstr[lenidx] += *t;

                                        // check for moving coefficients up
                                        // necessary condition: MUL-expression with 1 parent
                                        if (move_coeff && *t!=1+ABS(*(t+*t-1)) && *(instr[i]+1)!=OPER_COMMA &&
                                                        *(t+1)==EXTRASYMBOL && numpar[*(t+3)]==1 && *(instr[*(t+3)]+1)==OPER_MUL) {

                                                // coefficient is always the first term (that's how Horner+generate works)
                                                const WORD *t1 = instr[*(t+3)]+OPTHEAD;
                                                const WORD *t2 = t1+*t1;

                                                if (*t1 == 1+ABS(*(t1+*t1-1))) {
                                                        // t1 pointer to a coefficient, so move it

                                                        // remove old coefficient of 1
                                                        WORD *t3 = &*newinstr.end();                                         //
                                                        int sign2 = SGN(t3[-1]);                                                  //
                                                        newinstr.erase(newinstr.end()-3, newinstr.end());
                                                        // count number of arguments; iff it is 2 move the (extra)symbol too
                                                        int numargs=0;
                                                        for (const WORD *tt=t1; *tt!=0; tt+=*tt) {
                                                                numargs++;
                                                        }
                                                        if (numargs==2 && *(t2+4)==1) {
                                                                // replace (extra)symbol
                                                                newinstr[newinstr.size()-4] = *(t2+1);
                                                                newinstr[newinstr.size()-3] = *(t2+2);
                                                                newinstr[newinstr.size()-2] = *(t2+3);
                                                                newinstr[newinstr.size()-1] = *(t2+4);
                                                                sign2 *= SGN(*(t2+*t2-1));                        // was t2[7]

                                                                // ignore this expression from now on
                                                                skip[*(t+3)]=true;
                                                                if (*(t2+1)==EXTRASYMBOL)
                                                                        renum_par[*(t+3)] = *(t2+3);
                                                        }
                                                        else {
                                                                // otherwise, ignore coefficient from now on
                                                                // we need to collect the signs of the terms
                                                                // first and set them to one. This was forgotten
                                                                // before. Gave occasional errors.
                                                                if ( numargs > 2 || ( numargs == 2 && t2[4] > 1 ) ) {
                                                                        for (WORD *tt=(WORD *)t2; *tt!=0; tt+=*tt) {
                                                                                if ( tt[*tt-1] < 0 ) {
                                                                                        tt[*tt-1] = -tt[*tt-1];
                                                                                        sign2 = -sign2;
                                                                                }
                                                                        }
                                                                }
                                                                skipcoeff[*(t+3)]=true;
                                                        }

                                                        // add new coefficient
                                                        newinstr.insert(newinstr.end(), t1+1, t1+*t1);
                                                        newinstr.back() *= sign2;
                                                        newinstr[thislenidx] += ABS(newinstr.back()) - 3;
                                                        newinstr[lenidx] += ABS(newinstr.back()) - 3;
                                                }
                                        }
                                }
                        }
                }

                // if something has been copied, add trailing zero
                if (!first_copy) {
                        newinstr.push_back(0);
                        newinstr[lenidx]++;
                }
        }

        // renumber the expressions to 0,1,2,..,; only keep expressions with
        // skip=false which are their own parent after a renumbering in case
        // of moved coefficients

        // find renumber scheme
        vector<int> renum(n,-1);
        int next=0;
        for (int i=0; i<n; i++)
                if (!skip[i] && renum_par[par[i]]==i) renum[renum_par[i]]=next++;

        // find new instruction index
        tbegin = &*newinstr.begin();
        tend = tbegin+newinstr.size();
        for (const WORD *t=tbegin; t!=tend; t+=*(t+2))
                instr[renum[*t]] = t;

        // renumbering expressions and sort them lexicographically (in
        // new_instr they are in the preorder of the expression tree)
        vector<WORD> sortinstr;
        for (int i=0; i<next; i++) {
                int idx = sortinstr.size();
                sortinstr.insert(sortinstr.end(), instr[i], instr[i]+*(instr[i]+2));

                sortinstr[idx] = renum[sortinstr[idx]];

                // renumber content of an expression
                for (WORD *t2=&sortinstr[idx]+3; *t2!=0; t2+=*t2)
                        if (*t2!=1+ABS(*(t2+*t2-1)) && *(t2+1)==EXTRASYMBOL)
                                *(t2+3) = renum[*(t2+3)];
        }

#ifdef DEBUG_MORE
        MesPrint ("*** [%s, w=%w] DONE: merge_operators", thetime_str().c_str());
#endif

        return sortinstr;
}

/*
          #] merge_operators : 
          #[ class Optimization :
*/

/**  class Optimization
 *
 *         Description
 *         ===========
 *         This object represents an optimization. Its type is a number in
 *         the range 0 to 5. Depending on this type, the variables arg1, arg2
 *         and coeff indicate:
 *
 *         type==0 : optimization of the form x[arg1] ^ arg2 (coeff=empty)
 *         type==1 : optimization of the form x[arg1] * x[arg2] (coeff=empty)
 *         type==2 : optimization of the form x[arg1] * coeff (arg2=0)
 *         type==3 : optimization of the form x[arg1] + coeff (arg2=0)
 *         type==4 : optimization of the form x[arg1] + x[arg2] (coeff=empty)
 *         type==5 : optimization of the form x[arg1] - x[arg2] (coeff=empty)
 *
 *         Here, "x[arg]" represents a symbol (if positive) or an
 *         extrasymbol (if negative). The represented symbol's id is
 *         ABS(x[arg])-1.
 *
 *         "eqns" is a list of equation, where this optimization can be
 *         performed.
 *
 *         "improve" is the total improvement of this optimization.
 */
class optimization {
public:
        int type, arg1, arg2, improve;
        vector<WORD> coeff;
        vector<int> eqnidxs;

        bool operator< (const optimization &a) const {
                if (arg1 != a.arg1) return arg1 < a.arg1;
                if (arg2 != a.arg2) return arg2 < a.arg2;
                if (type != a.type) return type < a.type;
                return coeff < a.coeff;
        }
};

/*
          #] class Optimization : 
          #[ find_optimizations :
*/

/**  Find optimizations
 *
 *         Description
 *         ===========
 *         This method find all optimization of the form described in "class
 *         Optimization". It process every equation, looking for possible
 *         optimizations and stores them in a fast-access data structure to
 *         count the total improvement of an optimization.
 */
vector<optimization> find_optimizations (const vector<WORD> &instr) {

#ifdef DEBUG_GREEDY
        MesPrint ("*** [%s, w=%w] CALL: find_optimizations", thetime_str().c_str());
#endif

//        #[ Startup :
        // the resulting vector of optimizations
        vector<optimization> res;

        // a map to count the improvement of an optimization; the
        // improvement is stored as a vector<int> with equation numbers
        map<optimization, vector<int> > cnt;

        // a map to identify coefficients
        map<vector<int>,int> idx_coeff;

        const WORD *ebegin = &*instr.begin();
        const WORD *eend = ebegin+instr.size();

        for (int optim_type=0; optim_type<=4; optim_type++) {

                cnt.clear();
                idx_coeff.clear();

          optimization optim;
                optim.type = optim_type;
//        #] Startup : 

//        #[ type 0 :          find optimizations of the form z=x^n (optim.type==0)
                if (optim_type == 0) {
                        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                                if (*(e+1) != OPER_MUL) continue;
                                for (const WORD *t=e+3; *t!=0; t+=*t) {
                                        if (*t == ABS(*(t+*t-1))+1) continue;
                                        if (*(t+4) > 1) {
                                                optim.arg1 = (*(t+1)==SYMBOL ? 1 : -1) * (*(t+3) + 1);
                                                optim.arg2 = *(t+4);
                                                cnt[optim].push_back(e-ebegin);
                                        }
                                }
                        }
                }
//        #] type 0 : 
//        #[ type 1 :         find optimizations of the form z=x*y (optim.type==1)
                if (optim_type == 1) {
                        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                                if (*(e+1) != OPER_MUL) continue;

                                for (const WORD *t1=e+3; *t1!=0; t1+=*t1) {
                                        if (*t1 == ABS(*(t1+*t1-1))+1) continue;
                                        int x1 = (*(t1+1)==SYMBOL ? 1 : -1) * (*(t1+3) + 1);

                                        for (const WORD *t2=t1+*t1; *t2!=0; t2+=*t2) {
                                                if (*t2 == ABS(*(t2+*t2-1))+1) continue;
                                                int x2 = (*(t2+1)==SYMBOL ? 1 : -1) * (*(t2+3) + 1);

                                                if (*(t1+4) == *(t2+4)) {
                                                        optim.arg1 = x1;
                                                        optim.arg2 = x2;
                                                        if (optim.arg1 > optim.arg2)
                                                                swap(optim.arg1, optim.arg2);

                                                        if (*(t1+4) == 1)
                                                                cnt[optim].push_back(e-ebegin);
                                                        else {
                                                                // E=x^n*y^n -> z=x*y; E=z^n is double improvement
                                                                cnt[optim].push_back(e-ebegin);
                                                                cnt[optim].push_back(e-ebegin);
                                                        }
                                                }
                                        }
                                }
                        }
                }
//        #] type 1 : 
//        #[ type 2 :         find optimizations of the form z=c*x (optim.type==2)
                if (optim_type == 2) {
                        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {

                                if (*(e+1) == OPER_ADD) {
                                        // in ADD-equation

                                        for (const WORD *t=e+3; *t!=0; t+=*t) {
                                                if (*t == ABS(*(t+*t-1))+1) continue;

                                                if (*(t+4)==1) {
                                                        if (ABS(*(t+*t-1))==3 && *(t+*t-2)==1 && *(t+*t-3)==1) continue;

                                                        optim.coeff = vector<WORD>(t+*t-ABS(*(t+*t-1)), t+*t);
                                                        optim.coeff.back() = ABS(optim.coeff.back());

                                                        optim.arg1 = (*(t+1)==SYMBOL ? 1 : -1) * (*(t+3) + 1);
                                                        optim.arg2 = 0;

                                                        cnt[optim].push_back(e-ebegin);
                                                }
                                        }
                                }
                                else if (*(e+1) == OPER_MUL) {
                                        // in MUL-equation
                                        optim.coeff.clear();

                                        for (const WORD *t=e+3; *t!=0; t+=*t)
                                                if (*t == ABS(*(t+*t-1))+1) {
                                                        if (ABS(*(t+*t-1))==3 && *(t+*t-2)==1 && *(t+*t-3)==1) continue;
                                                        optim.coeff = vector<WORD>(t+*t-ABS(*(t+*t-1)), t+*t);
                                                        optim.coeff.back() = ABS(optim.coeff.back());
                                                }

                                        if (!optim.coeff.empty())
                                                for (const WORD *t=e+3; *t!=0; t+=*t) {
                                                        if (*t == ABS(*(t+*t-1))+1) continue;
                                                        if (*(t+4) != 1) continue;
                                                        optim.arg1 = (*(t+1)==SYMBOL ? 1 : -1) * (*(t+3) + 1);
                                                        optim.arg2 = 0;
                                                        cnt[optim].push_back(e-ebegin);
                                                }
                                }
                        }
                }
//        #] type 2 : 
//        #[ type 3 :         find optimizations of the form z=x+c (optim.type==3)
                if (optim_type == 3) {
                        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {

                                if (*(e+1) != OPER_ADD) continue;

                                optim.coeff.clear();

                                for (const WORD *t=e+3; *t!=0; t+=*t)
                                        if (*t == ABS(*(t+*t-1))+1) {
                                                optim.coeff = vector<WORD>(t+*t-ABS(*(t+*t-1)), t+*t);
                                        }

                                if (!optim.coeff.empty()) {
                                        for (const WORD *t=e+3; *t!=0; t+=*t) {
                                                if (*t == ABS(*(t+*t-1))+1) continue;
                                                if (ABS(*(t+*t-1))!=3 || *(t+*t-2)!=1 || *(t+*t-3)!=1) continue;
                                                if (*(t+*t-1)==-3) optim.coeff.back() *= -1;

                                                optim.arg1 = (*(t+1)==SYMBOL ? 1 : -1) * (*(t+3) + 1);
                                                optim.arg2 = 0;
                                                cnt[optim].push_back(e-ebegin);
                                                if (*(t+*t-1)==-3) optim.coeff.back() *= -1;
                                        }
                                }
                        }
                }
//        #] type 3 : 
//        #[ type 4,5 : find optimizations of the form z=x+y or z=x-y (optim.type==4 or 5)
                if (optim_type == 4) {
                        for (const WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                                if (*(e+1) != OPER_ADD) continue;

                                for (const WORD *t1=e+3; *t1!=0; t1+=*t1) {
                                        if (*t1 == ABS(*(t1+*t1-1))+1) continue;
                                        int x1 = (*(t1+1)==SYMBOL ? 1 : -1) * (*(t1+3) + 1);

                                        for (const WORD *t2=t1+*t1; *t2!=0; t2+=*t2) {
                                                if (*t2 == ABS(*(t2+*t2-1))+1) continue;
                                                int x2 = (*(t2+1)==SYMBOL ? 1 : -1) * (*(t2+3) + 1);

                                                int sign1 = SGN(*(t1+*t1-1));
                                                int sign2 = SGN(*(t2+*t2-1));

                                                if (BigLong((UWORD *)t1+5, ABS(*(t1+*t1-1))-1, (UWORD *)t2+5, ABS(*(t2+*t2-1))-1) == 0) {

                                                        optim.type = (sign1 * sign2 == 1 ? 4 : 5); // optimization type
                                                        optim.arg1 = x1;
                                                        optim.arg2 = x2;
                                                        if (optim.arg1 > optim.arg2) {
                                                                swap(optim.arg1, optim.arg2);
                                                        }

                                                        if (ABS(*(t1+*t1-1))==3 && *(t1+*t1-2)==1 && *(t1+*t1-3)==1)
                                                                cnt[optim].push_back(e-ebegin);
                                                        else {
                                                                // E=2x+2y -> z=x+y; E=2z is improvement bby itself
                                                                cnt[optim].push_back(e-ebegin);
                                                                cnt[optim].push_back(e-ebegin);
                                                        }
                                                }
                                        }
                                }
                        }
                }
//        #] type 4,5 : 
//        #[ add :

                // add optimizations with positive improvement to the result
                for (map<optimization, vector<int> >::iterator i=cnt.begin(); i!=cnt.end(); i++) {
                        int improve = i->second.size() - 1;
                        if (improve > 0) {
                                res.push_back(i->first);
                                res.back().improve = improve;
                                res.back().eqnidxs = i->second;

                                // remove duplicates, that were add to get the correct improvement
                                res.back().eqnidxs.erase(unique(res.back().eqnidxs.begin(), res.back().eqnidxs.end()), res.back().eqnidxs.end());
                        }
                }
        }
//        #] add : 

#ifdef DEBUG_GREEDY
        MesPrint ("*** [%s, w=%w] DONE: find_optimizations",thetime_str().c_str());
#endif

        return res;
}

/*
          #] find_optimizations : 
          #[ do_optimization :
*/

/**  Do optimization
 *
 *         Description
 *         ===========
 *         This method performs an optimization. It scans through the
 *         equations of "optim.eqnidxs" and looks in which this optimization
 *         can still be performed (due to other performed optimizations this
 *         isn't always the case). If possible, it substitutes the common
 *         subexpression by a new extra symbol numbered "newid". Finally,
 *         the new extrasymbol is defined accordingly.
 *
 *         Substitutions may lead to trivial equations of the form "Zi=Zj",
 *         but these are removed in the end of the method. The method returns
 *         whether the substitution has been done once or more (or not).
 */

bool do_optimization (const optimization optim, vector<WORD> &instr, int newid) {

//        #[ Debug code :
#ifdef DEBUG_GREEDY
        if (optim.type==0)
                MesPrint ("*** [%s, w=%w] CALL: do_optimization(improve=%d, %c%d^%d)",
                  thetime_str().c_str(), optim.improve,
                optim.arg1>0?'x':'Z', ABS(optim.arg1)-1, optim.arg2);
        else if (optim.type==1 || optim.type>=4)
                MesPrint ("*** [%s, w=%w] CALL: do_optimization(improve=%d, %c%d%c%c%d)",
                thetime_str().c_str(), optim.improve,
                optim.arg1>0?'x':'Z', ABS(optim.arg1)-1,
                optim.type==1 ? '*' : optim.type==4 ? '+' : '-',
                optim.arg2>0?'x':'Z', ABS(optim.arg2)-1);
        else {
        WORD n = optim.coeff.back()/2;
        UBYTE num[BITSINWORD*ABS(n)], den[BITSINWORD*ABS(n)];
        PrtLong((UWORD *)&optim.coeff[0], n, num);
        PrtLong((UWORD *)&optim.coeff[ABS(n)], ABS(n), den);
                MesPrint ("*** [%s, w=%w] CALL: do_optimization(improve=%d, %c%d%c%s/%s)",
                thetime_str().c_str(), optim.improve,
                optim.arg1>0?'x':'Z', ABS(optim.arg1)-1,
                optim.type==2 ? '*' : '+', num,den);
  }
#endif
//        #] Debug code : 

        bool substituted = false;
        WORD *ebegin = &*instr.begin();

//        #[ type 0 : substitution of the form z=x^n (optim.type==0)
        if (optim.type == 0) {

                int vartypex = optim.arg1>0 ? SYMBOL : EXTRASYMBOL;
                int varnumx  = ABS(optim.arg1) - 1;
                int n = optim.arg2;

                for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                        WORD *e = ebegin + optim.eqnidxs[i];
                        if (*(e+1) != OPER_MUL) continue;

                        // scan through equation
                        for (WORD *t=e+3; *t!=0; t+=*t) {
                                if (*t == ABS(*(t+*t-1))+1) continue;
                                if (*(t+1)==vartypex &&
                                                *(t+3)==varnumx &&
                                                *(t+4) % n == 0) {

                                        // substitute
                                        *(t+1) = EXTRASYMBOL;
                                        *(t+3) = newid;
                                        *(t+4) /= n;

                                        substituted = true;
                                }
                        }
                }

                if (!substituted) {
#ifdef DEBUG_GREEDY
                        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=false", thetime_str().c_str(), optim.improve);
#endif
                        return false;
                }

                // add extra equation (Tnew = x^n)
                instr.push_back(newid);    // eqn.nr
                instr.push_back(OPER_MUL); // operator
                instr.push_back(12);           // total length
                instr.push_back(8);            // term length
                instr.push_back(vartypex); // (extra)symbol
                instr.push_back(4);            // symbol length
                instr.push_back(varnumx);  // symbol id
                instr.push_back(n);            // power
                instr.push_back(1);
                instr.push_back(1);            // coefficient 1
                instr.push_back(3);
                instr.push_back(0);            // trailing 0
        }
//        #] type 0 : 
//        #[ type 1 : substitution of the form z=x*y (optim.type==1)
        if (optim.type == 1) {

                int vartypex = optim.arg1>0 ? SYMBOL : EXTRASYMBOL;
                int varnumx  = ABS(optim.arg1) - 1;
                int vartypey = optim.arg2>0 ? SYMBOL : EXTRASYMBOL;
                int varnumy  = ABS(optim.arg2) - 1;

                for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                        WORD *e = ebegin + optim.eqnidxs[i];
                        if (*(e+1) != OPER_MUL) continue;

                        // scan through equation
                        int powx=0, powy=0;
                        for (WORD *t=e+3; *t!=0; t+=*t) {
                                if (*t == ABS(*(t+*t-1))+1) continue;
                                if (*(t+1)==vartypex && *(t+3)==varnumx) powx = *(t+4);
                                if (*(t+1)==vartypey && *(t+3)==varnumy) powy = *(t+4);
                        }

                        // substitute if found
                        if (powx>0 && powy>0 && powx==powy) {

                                WORD sign = 1;
                                WORD *newt = e+3;

                                for (WORD *t=e+3; *t!=0;) {
                                        int dt=*t;

                                        if (*t == ABS(*(t+*t-1))+1 ||
                                                        (!(*(t+1)==vartypex && *(t+3)==varnumx) &&
                                                         !(*(t+1)==vartypey && *(t+3)==varnumy))) {
                                                memmove(newt, t, *t*sizeof(WORD));
                                                newt += dt;
                                        }
                                        else {
                                                sign *= SGN(*(t+*t-1));
                                        }

                                        t+=dt;
                                }

                                *newt++ = 8;                   // term length
                                *newt++ = EXTRASYMBOL; // extrasymbol
                                *newt++ = 4;                   // symbol length
                                *newt++ = newid;           // symbol id
                                *newt++ = powx;            // power
                                *newt++ = 1;
                                *newt++ = 1;                   // coefficient +/-1
                                *newt++ = 3*sign;
                                *newt++ = 0;                   // trailing 0

                                substituted = true;
                        }
                }

                if (!substituted) {
#ifdef DEBUG_GREEDY
                        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=false", thetime_str().c_str(), optim.improve);
#endif
                        return false;
                }

                // add extra equation (Tnew = x*y)
                instr.push_back(newid);    // eqn.nr
                instr.push_back(OPER_MUL); // operator
                instr.push_back(20);           // total length
                instr.push_back(8);            // LHS length
                instr.push_back(vartypex); // (extra)symbol
                instr.push_back(4);            // symbol length
                instr.push_back(varnumx);  // symbol id
                instr.push_back(1);            // power 1
                instr.push_back(1);
                instr.push_back(1);            // coefficient 1
                instr.push_back(3);
                instr.push_back(8);            // RHS length
                instr.push_back(vartypey); // (extra)symbol
                instr.push_back(4);            // symbol length
                instr.push_back(varnumy);  // symbol id
                instr.push_back(1);            // power 1
                instr.push_back(1);
                instr.push_back(1);            // coefficient 1
                instr.push_back(3);
                instr.push_back(0);            // trailing 0
        }
//        #] type 1 : 
//        #[ type 2 : substitution of the form z=c*x (optim.type==2)

        if (optim.type == 2) {
                int vartype = optim.arg1>0 ? SYMBOL : EXTRASYMBOL;
                int varnum        = ABS(optim.arg1) - 1;

                WORD ncoeff = optim.coeff.back();

                // scan through equations
                for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                        WORD *e = ebegin + optim.eqnidxs[i];

                        if (*(e+1) == OPER_ADD) {
                                // scan through ADD-equation
                                for (WORD *t=e+3; *t!=0; t+=*t) {

                                        if (*t == ABS(*(t+*t-1))+1) continue;

                                        if (*(t+1)==vartype && ABS(*(t+3))==varnum && *(t+4)==1 &&
                                                        BigLong((UWORD *)&optim.coeff[0],ncoeff-1,
                                                                                        (UWORD *)t+*t-ABS(*(t+*t-1)),ABS(*(t+*t-1))-1) == 0) {
                                                // substitute

                                                int sign = SGN(*(t+*t-1));

                                                WORD *tend = t;
                                                while (*tend!=0) tend+=*tend;
                                                WORD nmove = tend - t - *t;
                                                memmove(t, t+*t, nmove*sizeof(WORD));
                                                t += nmove;

                                                *t++ = 8;                        // term length
                                                *t++ = EXTRASYMBOL; // (extra)symbol
                                                *t++ = 4;                        // symbol length
                                                *t++ = newid;                // symbol id
                                                *t++ = 1;                        // power of 1
                                                *t++ = 1;
                                                *t++ = 1;                        // coefficient of +/-1
                                                *t++ = 3 * sign;
                                                *t++ = 0;                        // trailing 0

                                                substituted = true;
                                                break;
                                        }
                                }
                        }
                        else if (*(e+1) == OPER_MUL) {

                                bool coeff_match=false, var_match=false;
                                int sign = 1;

                                // scan through MUL-equation
                                for (WORD *t=e+3; *t!=0; t+=*t) {
                                        if (*t == ABS(*(t+*t-1))+1 && BigLong((UWORD *)&optim.coeff[0],ncoeff-1,
                                                                                                                                                                                                (UWORD *)t+*t-ABS(*(t+*t-1)),ABS(*(t+*t-1))-1) == 0) {
                                                coeff_match = true;
                                                sign *= SGN(*(t+*t-1));
                                        }
                                        else if (*(t+1)==vartype && ABS(*(t+3))==varnum && *(t+4)==1) {
                                                var_match = true;
                                                sign *= SGN(*(t+*t-1));
                                        }
                                }

                                // substitute if found
                                if (coeff_match && var_match) {

                                        WORD *newt = e+3;

                                        for (WORD *t=e+3; *t!=0;) {

                                                int dt=*t;

                                                if (*t!=ABS(*(t+*t-1))+1 && !(*(t+1)==vartype && ABS(*(t+3))==varnum && *(t+4)==1)) {
                                                        memmove(newt, t, dt*sizeof(WORD));
                                                        newt += dt;
                                                }

                                                t+=dt;
                                        }

                                        *newt++ = 8;                   // term length
                                        *newt++ = EXTRASYMBOL; // extrasymbol
                                        *newt++ = 4;                   // symbol length
                                        *newt++ = newid;           // symbol id
                                        *newt++ = 1;                   // power of 1
                                        *newt++ = 1;
                                        *newt++ = 1;                   // coefficient of +/-1
                                        *newt++ = 3 * sign;
                                        *newt++ = 0;                   // trailing 0

                                        substituted = true;
                                }
                        }
                }

                if (!substituted) {
#ifdef DEBUG_GREEDY
                        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=false", thetime_str().c_str(), optim.improve);
#endif
                        return false;
                }

                // add extra equation (Tnew = c*y)
                instr.push_back(newid);                    // eqn.nr
                instr.push_back(OPER_ADD);                   // operator
                instr.push_back(9+ABS(ncoeff));    // total length
                instr.push_back(5+ABS(ncoeff));    // term length
                instr.push_back(vartype);                   // (extra)symbol
                instr.push_back(4);                            // symbol length
                instr.push_back(varnum);                   // symbol id
                instr.push_back(1);                            // power of 1
                for (int i=0; i<ABS(ncoeff); i++)  // coefficient
                        instr.push_back(optim.coeff[i]);
                instr.push_back(0);                            // trailing 0
        }
//        #] type 2 : 
//        #[ type 3 : substitution of the form z=x+c (optim.type==3)
        if (optim.type == 3) {
                int vartype = optim.arg1>0 ? SYMBOL : EXTRASYMBOL;
                int varnum        = ABS(optim.arg1) - 1;

                WORD ncoeff = optim.coeff.back();

                // scan through equation
                for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                        WORD *e = ebegin + optim.eqnidxs[i];

                        if (*(e+1) != OPER_ADD) continue;

                        int coeff_match=0, var_match=0;

                        for (WORD *t=e+3; *t!=0; t+=*t) {
                                if (*t == ABS(*(t+*t-1))+1 && BigLong((UWORD *)&optim.coeff[0],ABS(ncoeff)-1,
                                                                                                                                                                                        (UWORD *)t+*t-ABS(*(t+*t-1)),ABS(*(t+*t-1))-1) == 0)
                                        coeff_match = SGN(ncoeff) * SGN(*(t+*t-1));
                                else if (*(t+1)==vartype && ABS(*(t+3))==varnum && *(t+4)==1)
                                        var_match = SGN(*(t+7));
                        }

                        // substitute if found (x+c and -x-c and matches)
                        if (coeff_match * var_match == 1) {

                                WORD *newt = e+3;

                                for (WORD *t=e+3; *t!=0;) {
                                        int dt=*t;
                                        if (*t!=ABS(*(t+*t-1))+1 && !(*(t+1)==vartype && ABS(*(t+3))==varnum && *(t+4)==1)) {
                                                memmove(newt, t, dt*sizeof(WORD));
                                                newt += dt;
                                        }
                                        t+=dt;
                                }

                                *newt++ = 8;                        // term length
                                *newt++ = EXTRASYMBOL;        // extrasymbol
                                *newt++ = 4;                        // symbol length
                                *newt++ = newid;                // symbol id
                                *newt++ = 1;                        // power of 1
                                *newt++ = 1;
                                *newt++ = 1;                        // coefficient of +/-1
                                *newt++ = 3*coeff_match;
                                *newt++ = 0;                        // trailing zero

                                substituted = true;
                        }
                }

                if (!substituted) {
#ifdef DEBUG_GREEDY
                        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=false", thetime_str().c_str(), optim.improve);
#endif
                        return false;
                }

                // add extra equation (Tnew = x+c)
                instr.push_back(newid);                   // eqn.nr
                instr.push_back(OPER_ADD);                  // operator
                instr.push_back(13+ABS(ncoeff));  // total length
                instr.push_back(8);                           // x-term length
                instr.push_back(vartype);                  // (extra)symbol
                instr.push_back(4);                           // symbol length
                instr.push_back(varnum);                  // symbol id
                instr.push_back(1);                           // power of 1
                instr.push_back(1);
                instr.push_back(1);                           // coefficient of 1
                instr.push_back(3);
                instr.push_back(ABS(ncoeff)+1);   // c-term length
                for (int i=0; i<ABS(ncoeff); i++) // coefficient
                        instr.push_back(optim.coeff[i]);
                instr.push_back(0);                           // trailing zero
        }
//        #] type 3 : 
//        #[ type 4,5 : substitution of the form z=x+y or z=x-y (optim.type=4 or 5)
        if (optim.type >= 4) {

                int vartypex = optim.arg1>0 ? SYMBOL : EXTRASYMBOL;
                int varnumx  = ABS(optim.arg1) - 1;
                int vartypey = optim.arg2>0 ? SYMBOL : EXTRASYMBOL;
                int varnumy  = ABS(optim.arg2) - 1;

                // scan through equations
                for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                        WORD *e = ebegin + optim.eqnidxs[i];

                        if (*(e+1) != OPER_ADD) continue;

                        const WORD *coeffx=NULL, *coeffy=NULL;
                        WORD ncoeffx=0,ncoeffy=0;

                        // looks for terms
                        for (WORD *t=e+3; *t!=0; t+=*t) {
                                if (*t == ABS(*(t+*t-1))+1) continue; // constant
                                if (*(t+1)==vartypex && *(t+3)==varnumx && *(t+4)==1) {
                                        coeffx = t+5;
                                        ncoeffx = *(t+*t-1);
                                }
                                if (*(t+1)==vartypey && *(t+3)==varnumy && *(t+4)==1) {
                                        coeffy = t+5;
                                        ncoeffy = *(t+*t-1);
                                }
                        }

                        // check signs (type=4: x+y and -x-y, type=5: x-y and -x+y) ??????
                        // check signs (type=4: x+y, type=5: x-y) !!!!!!!!!!
                        if (SGN(ncoeffx) * SGN(ncoeffy) * (optim.type==4 ? 1 : -1) == 1) {
                                // check absolute value of coefficients
                                if (BigLong((UWORD *)coeffx, ABS(ncoeffx)-1, (UWORD *)coeffy, ABS(ncoeffy)-1) == 0) {
                                        // substitute
                                        vector<WORD> coeff(coeffx, coeffx+ABS(ncoeffx));

                                        WORD *newt = e+3;
/*
if ( optim.type == 5 ) {
        while ( *newt ) newt+=*newt;
        int i = (newt - e) - 3;
        MesPrint("        < %a",i,e+3);
        newt = e+3;
}
*/
                                        for (WORD *t=e+3; *t!=0;) {
                                                int dt=*t;
                                                if (*t == ABS(*(t+*t-1))+1 ||
                                                                (!(*(t+1)==vartypex && *(t+3)==varnumx) &&
                                                                 !(*(t+1)==vartypey && *(t+3)==varnumy))) {
                                                        memmove(newt, t, dt*sizeof(WORD));
                                                        newt += dt;
                                                }
                                                t+=dt;
                                        }

                                        *newt++ = 5 + ABS(ncoeffx);           // term length
                                        *newt++ = EXTRASYMBOL;                          // extrasymbol
                                        *newt++ = 4;                                          // symbol length
                                        *newt++ = newid;                                  // symbol id
                                        *newt++ = 1;                                          // power of 1
                                        for (int j=0; j<ABS(ncoeffx); j++)// coefficient
                                                *newt++ = coeff[j];
                                        *newt++ = 0;                                          // trailing 0
                                        substituted = true;
/*
if ( optim.type == 5 ) {
        int i = (newt - e) - 4;
        MesPrint("        > %a",i,e+3);
}
*/
                                }
                        }
                }

                if (!substituted) {
#ifdef DEBUG_GREEDY
                        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=false", thetime_str().c_str(), optim.improve);
#endif
                        return false;
                }
/*
if ( optim.type == 5 )
MesPrint ("improve=%d, %c%d%c%c%d)", optim.improve,
        optim.arg1>0?'x':'Z', ABS(optim.arg1)-1,
        optim.type==1 ? '*' : optim.type==4 ? '+' : '-',
        optim.arg2>0?'x':'Z', ABS(optim.arg2)-1);
*/
                // add extra equation (Tnew = x+/-y)
                instr.push_back(newid);    // eqn.nr
                instr.push_back(OPER_ADD); // operator
                instr.push_back(20);           // total length
                instr.push_back(8);            // term length
                instr.push_back(vartypex); // (extra)symbol
                instr.push_back(4);            // symbol length
                instr.push_back(varnumx);  // symbol id
                instr.push_back(1);            // power of 1
                instr.push_back(1);
                instr.push_back(1);            // coefficient of 1
                instr.push_back(3);
                instr.push_back(8);            // term length
                instr.push_back(vartypey); // (extra)symbol
                instr.push_back(4);            // symbol length
                instr.push_back(varnumy);  // symbol id
                instr.push_back(1);            // power of 1
                instr.push_back(1);
                instr.push_back(1);            // coefficient of +/-1
                instr.push_back(3*(optim.type==4?1:-1));
                instr.push_back(0);            // trailing 0
        }
//        #] type 4,5 : 
//        #[ trivial :  remove trivial equations of the form Zi = +/-Zj
        vector<int> renum(newid+1, 0);
        bool do_renum=false;

        // vector may be moved when it is extended
        ebegin = &*instr.begin();

        for (int i=0; i<(int)optim.eqnidxs.size(); i++) {
                WORD *e = ebegin + optim.eqnidxs[i];
                WORD *t = e+3;
                if (*t==0) continue;          // empty (removed) equation
                if (*(t+*t)!=0) continue; // more than 1 term
                if (*t!=8) continue;          // term not of correct form
                if (*(t+4)!=1) continue;  // power != 1
                if (*(t+5)!=1 || *(t+6)!=1) continue; // coeff != +/-1

                // trivial term, so renumber this one
                renum[*e] = SGN(*(t+7)) * (*(t+3) + 1);
                do_renum = true;

                // remove equation
                *t=0;
        }
//        #] trivial : 
//        #[ renumbering :

        // there are renumberings to be done, so loop through all equations
        if (do_renum) {
                WORD *eend = ebegin+instr.size();

                for (WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                        for (WORD *t=e+3; *t!=0; t+=*t) {
                                if (*t == ABS(*(t+*t-1))+1) continue;
                                if (*(t+1)==EXTRASYMBOL && renum[*(t+3)]!=0) {
                                        *(t+*t-1) *= SGN(renum[*(t+3)]);
                                        *(t+3) = ABS(renum[*(t+3)]) - 1;
                                }
                        }
                }
        }
//        #] renumbering : 

#ifdef DEBUG_GREEDY
        MesPrint ("*** [%s, w=%w] DONE: do_optimization : res=true", thetime_str().c_str(), optim.improve);
#endif

        return true;
}

/*
          #] do_optimization : 
          #[ partial_factorize :
*/

/**  Partial factorization of instructions
 *
 *         Description
 *         ===========
 *         This method performs partial factorization of instructions. In
 *         particular the following instructions
 *
 *           Z1 = x*a*b
 *           Z2 = x*c*d*e
 *           Z3 = 2*x + Z1 + Z2 + more
 *
 *         are replaced by
 *
 *           Z1 = a*b
 *           Z2 = c*d*e
 *           Z3 = Zj + more
 *           Zi = 2 + Z1 + Z2
 *           Zj = x*Zi
 *
 *         Here it is necessary that no other equations refer to Z1 and
 *         Z2. The generation of trivial instructions (Zi=Zj or Zi=x) is
 *         prevented.
 */
int partial_factorize (vector<WORD> &instr, int n, int improve) {

#ifdef DEBUG_GREEDY
        MesPrint ("*** [%s, w=%w] CALL: partial_factorize (n=%d)", thetime_str().c_str(), n);
#endif

        GETIDENTITY;

        // get starting positions of instructions
        vector<int> instr_idx(n);
        WORD *ebegin = &*instr.begin();
        WORD *eend = ebegin+instr.size();
        for (WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                instr_idx[*e] = e - ebegin;
        }

        // get reference counts
/*
 *        The next construction replaces the vector construction which is
 *        rather costly for valgrind (and maybe also in normal running)
 */
        int nmax = 2*n;
        WORD *numpar = (WORD *)Malloc1(nmax*sizeof(WORD),"numpar");
        for ( int i = 0; i < nmax; i++ ) numpar[i] = 0;
//        vector<WORD> numpar(n);
        for (WORD *e=ebegin; e!=eend; e+=*(e+2))
                for (WORD *t=e+3; *t!=0; t+=*t) {
                        if (*t == ABS(*(t+*t-1))+1) continue;
                        if (*(t+1) == EXTRASYMBOL) numpar[*(t+3)]++;
                }

        // find factorizable expressions
        for (int i=0; i<n; i++) {
                WORD *e = &*instr.begin() + instr_idx[i];
                if (*(e+1) != OPER_ADD) continue;

                // count symbol occurrences
                map<WORD,WORD> cnt; // 1-indexed, <0:EXTRASYMBOL, >0:SYMBOL

                for (WORD *t=e+3; *t!=0; t+=*t) {
                        if (*t==ABS(*(t+*t-1))+1) continue;

                        // count symbols in t
                        if (*(t+4)==1)
                                cnt[(*(t+1)==SYMBOL ? 1 : -1) * (*(t+3)+1)]++;

                        // count symbols in extrasymbols of t
                        if (*(t+1)==EXTRASYMBOL && *(t+4)==1 && numpar[*(t+3)]==1) {
                                WORD *t2 = &*instr.begin() + instr_idx[*(t+3)];
                                if (*(t2+1) != OPER_MUL) continue;
                                for (t2+=3; *t2!=0; t2+=*t2) {
                                        if (*t2 == ABS(*(t2+*t2-1))+1) continue;
                                        if (*(t2+4)==1)
                                                cnt[(*(t2+1)==SYMBOL ? 1 : -1) * (*(t2+3)+1)]++;
                                }
                        }
                }

                // find most-occurring symbol
                WORD x=0, best=0;
                for (map<WORD,WORD>::iterator it=cnt.begin(); it!=cnt.end(); it++)
                        if (it->second > best) { x=it->first; best=it->second; }

                // occurrence>=2 and occurrence>improve, so factorize

                if (best>=2 && best>improve) {
                        // initialize new equation (Zi from example above)
                        vector<WORD> new_eqn;
                        new_eqn.push_back(n);
                        new_eqn.push_back(OPER_ADD);
                        new_eqn.push_back(0); // length

                        WORD dt;
                        WORD *newt=e+3;
                        for (WORD *t=e+3; *t!=0; t+=dt) {
                                dt = *t;
                                bool keep=true;

                                if (*t!=ABS(*(t+*t-1))+1) {

                                        // factorized symbol is in t itself
                                        if (*(t+4)==1) {
                                                WORD y = (*(t+1)==SYMBOL ? 1 : -1) * (*(t+3)+1);
                                                if (y==x) {
                                                        new_eqn.push_back(*t-4);
                                                        new_eqn.insert(new_eqn.end(), t+5, t+dt);
                                                        keep=false;
                                                }
                                        }

                                        // look in extrasymbol of t with ref.count=1
                                        if (*(t+1)==EXTRASYMBOL && *(t+4)==1 && numpar[*(t+3)]==1) {
                                                WORD *t2 = &*instr.begin() + instr_idx[*(t+3)];
                                                if (*(t2+1) == OPER_MUL) {
                                                        bool has_x=false;
                                                        for (t2+=3; *t2!=0; t2+=*t2) {
                                                                if (*t2 == ABS(*(t2+*t2-1))+1) continue;
                                                                WORD y = (*(t2+1)==SYMBOL ? 1 : -1) * (*(t2+3)+1);
                                                                // extrasymbol has factorized symbol
                                                                if (y==x && *(t2+4)==1) {
                                                                        has_x=true;
                                                                        // copy remaining part
                                                                        WORD *tend=t2+*t2;
                                                                        WORD sign = SGN(*(tend-1));
                                                                        while (*tend!=0) tend+=*tend;
                                                                        int dt2 = tend - (t2+*t2);
                                                                        memmove(t2, t2+*t2, (dt2+1)*sizeof(WORD));
                                                                        t2 += dt2;
                                                                        *(t2-1) *= sign;
                                                                        break;
                                                                }
                                                        }
                                                        if (has_x) {
                                                                // extrasymbol has x, so add it to new equation
                                                                keep=false;
                                                                int thisidx=new_eqn.size();
                                                                new_eqn.insert(new_eqn.end(), t, t+dt);
                                                                t2 = &*instr.begin() + instr_idx[*(t+3)] + 3;
                                                                // if becomes trivial, substitute the term
                                                                if (*(t2+*t2)==0) {
                                                                        // it's a number
                                                                        if (*t2 == ABS(*(t2+*t2-1))+1) {
                                                                                if (ABS(new_eqn[new_eqn.size()-1])==3 && new_eqn[new_eqn.size()-2]==1 && new_eqn[new_eqn.size()-3]==1) {
                                                                                        // original equation has coefficient of +/-1, so replace it
                                                                                        WORD sign = SGN(new_eqn.back());
                                                                                        new_eqn.erase(new_eqn.begin()+thisidx, new_eqn.end());
                                                                                        new_eqn.insert(new_eqn.end(), t2, t2+*t2);
                                                                                        new_eqn.back() *= sign;
                                                                                        *t2 = 0;
                                                                                }
                                                                                else {
                                                                                        // two non-trivial coefficients, so multiply them
                                                                                        // note: untested code (found no way to trigger it)
                                                                                        UWORD *tmp = NumberMalloc("partial_factorize");
                                                                                        WORD ntmp=0;
                                                                                        MulRat(BHEAD (UWORD *)t2+*t2-ABS(*(t2+*t2-1)), *(t2+*t2-1),
                                                                                                                 (UWORD *)&*(new_eqn.end()-ABS(new_eqn.back())), new_eqn.back(),
                                                                                                                 tmp, &ntmp);
                                                                                        new_eqn.erase(new_eqn.begin()+thisidx, new_eqn.end());
                                                                                        new_eqn.push_back(ABS(ntmp)+1);
                                                                                        new_eqn.insert(new_eqn.end(), tmp, tmp+ABS(ntmp));
                                                                                        NumberFree(tmp,"partial_factorize");
                                                                                        *t2 = 0;
                                                                                }
                                                                        }
                                                                        else if (*(t2+4)==1) {
                                                                                // it's a variable
                                                                                new_eqn.back() *= SGN(*(t2+*t2-1));
                                                                                new_eqn[thisidx+1] = *(t2+1);
                                                                                new_eqn[thisidx+2] = *(t2+2);
                                                                                new_eqn[thisidx+3] = *(t2+3);
                                                                                new_eqn[thisidx+4] = *(t2+4);
                                                                                *t2 = 0;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }

                                // no x, so copy it
                                if (keep) {
                                        memmove(newt, t, dt*sizeof(WORD));
                                        newt += dt;
                                }
                        }

                        // finalize new equation
                        new_eqn.push_back(0);
                        new_eqn[2] = new_eqn.size();

                        bool empty = newt == e+3;
                        if ( n+1 >= nmax ) {
                                int i, newnmax = nmax*2;
                                WORD *newnumpar = (WORD *)Malloc1(newnmax*sizeof(WORD),"newnumpar");
                                for ( i = 0; i < n; i++ ) newnumpar[i] = numpar[i];
                                for ( ; i < newnmax; i++ ) newnumpar[i] = 0;
                                M_free(numpar,"numpar");
                                numpar = newnumpar;
                                nmax = newnmax;
                        }
//                        numpar.push_back(0);
                        n++;

                        // if original is not empty, add new equation (Zj) to it
                        // otherwise replace it later
                        if (!empty) {
                                *newt++ = 8;
                                *newt++ = EXTRASYMBOL;
                                *newt++ = 4;
                                *newt++ = n;
                                *newt++ = 1;
                                *newt++ = 1;
                                *newt++ = 1;
                                *newt++ = 3;
                                *newt++ = 0;
                        }

                        // add new equation to instructions
                        instr_idx.push_back(instr.size());
                        instr.insert(instr.end(), new_eqn.begin(), new_eqn.end());

                        // generate another new equation (Zj=x*Zi)
                        new_eqn.clear();
                        new_eqn.push_back(n);
                        new_eqn.push_back(OPER_MUL);
                        new_eqn.push_back(20);
                        new_eqn.push_back(8);

                        // add factorized symbol
                        if (x>0) {
                                new_eqn.push_back(SYMBOL);
                                new_eqn.push_back(4);
                                new_eqn.push_back(x-1);
                                new_eqn.push_back(1);
                        }
                        else {
                                new_eqn.push_back(EXTRASYMBOL);
                                new_eqn.push_back(4);
                                new_eqn.push_back(-x-1);
                                new_eqn.push_back(1);
                        }
                        new_eqn.push_back(1);
                        new_eqn.push_back(1);
                        new_eqn.push_back(3);
                        new_eqn.push_back(8);
                        // add new equation (Zi)
                        new_eqn.push_back(EXTRASYMBOL);
                        new_eqn.push_back(4);
                        new_eqn.push_back(n-1);
                        new_eqn.push_back(1);
                        new_eqn.push_back(1);
                        new_eqn.push_back(1);
                        new_eqn.push_back(3);
                        new_eqn.push_back(0);

                        if (!empty) {
                                // add new equation (Zj) to instructions
                                instr_idx.push_back(instr.size());
                                instr.insert(instr.end(), new_eqn.begin(), new_eqn.end());
                                if ( n+1 >= nmax ) {
                                        int i, newnmax = nmax*2;
                                        WORD *newnumpar = (WORD *)Malloc1(newnmax*sizeof(WORD),"newnumpar");
                                        for ( i = 0; i < n; i++ ) newnumpar[i] = numpar[i];
                                        for ( ; i < newnmax; i++ ) newnumpar[i] = 0;
                                        M_free(numpar,"numpar");
                                        numpar = newnumpar;
                                        nmax = newnmax;
                                }
//                                numpar.push_back(0);
                                n++;
                        }
                        else {
                                // replace e with Zj
                                e = &*instr.begin() + instr_idx[i];
                                e[1] = OPER_MUL;
                                memcpy(e+3, &new_eqn[3], (new_eqn.size()-3)*sizeof(WORD));
                        }

                        // decrease i, so this expression is factorized again if possible
                        i--;
                }
        }

#ifdef DEBUG_GREEDY
        MesPrint ("*** [%s, w=%w] DONE: partial_factorize (n=%d)", thetime_str().c_str(), n);
#endif
        M_free(numpar,"numpar");
        return n;
}

/*
          #] partial_factorize : 
          #[ optimize_greedy :
*/

/**  Optimize instructions greedily
 *
 *         Description
 *         ===========
 *         This method optimizes an expression greedily. It calls
 *         "find_optimizations" to obtain candidates and performs the best
 *         one(s) by calling "do_optimization".
 *
 *         How many different optimization are done, before
 *         "find_optimization" is called again, is determined by the
 *         settings "greedyminnum" and "greedymaxperc".
 *
 *         During the optimization process, sequences of zeroes are
 *         introduced in the instructions, since moving all instructions
 *         when one gets optimized, is very costly. Therefore, in the end,
 *         the instructions are "compressed" again to remove these extra
 *         zeroes.
 */
vector<WORD> optimize_greedy (vector<WORD> instr, LONG time_limit) {

#ifdef DEBUG
        int old_num_oper = count_operators(instr);
        MesPrint ("*** [%s, w=%w] CALL: optimize_greedy(numoper=%d)",
                                                thetime_str().c_str(), old_num_oper);
#endif

        LONG start_time = TimeWallClock(1);

        WORD *ebegin = &*instr.begin();
        WORD *eend = ebegin+instr.size();

        // store final equation, since it must be the last equation later
        int final_eqn_idx = 0;
        int next_eqn = 0;

        for (WORD *e=ebegin; e!=eend; e+=*(e+2)) {
                next_eqn = *e + 1;
                final_eqn_idx = e-ebegin;
        }
        // optimize instructions
        while (TimeWallClock(1)-start_time < time_limit) {
                int old_next_eqn = next_eqn;

                // find optimizations
                vector<optimization> optim = find_optimizations(instr);

                // add_eqnidxs contains modified equations, which might have to be updated later again
                vector<int> add_eqnidxs;

                // number of optimizations to do
                int num_do_optims = max(AO.Optimize.greedyminnum, (int)optim.size()*AO.Optimize.greedymaxperc/100);

                // if best improvement is one, do all optimizations
                int best_improve=0;
                for (int i=0; i<(int)optim.size(); i++)
                        best_improve = max(best_improve, optim[i].improve);
                if (best_improve <= 1)
                        num_do_optims = MAXPOSITIVE;
                // do a number of optimizations
                while (optim.size() > 0 && num_do_optims-- > 0) {

                        // find best optimization
                        int best=0;
                        best_improve=0;
                        for (int i=0; i<(int)optim.size(); i++)
                                if (optim[i].improve > best_improve) {
                                        best=i;
                                        best_improve=optim[i].improve;
                                }

                        // add extra equations
                        for (int i=0; i<(int)add_eqnidxs.size(); i++)
                                optim[best].eqnidxs.push_back(add_eqnidxs[i]);

                        // do optimization, update next_eqn if successful
                        int next_idx = instr.size();
                        if (do_optimization(optim[best], instr, next_eqn)) {
                                next_eqn++;
                                add_eqnidxs.push_back(next_idx);
                        }

                        optim.erase(optim.begin()+best);
                }

                // partially factorize with improve >= best_improve
                next_eqn = partial_factorize(instr, next_eqn, best_improve);

                // check whether nothing has changed
                if (next_eqn == old_next_eqn) break;
        }

        // add final equation to the back (must be by definition)
        instr.push_back(next_eqn);
        instr.insert(instr.end(), instr.begin()+final_eqn_idx+1, instr.begin()+final_eqn_idx+instr[final_eqn_idx+2]);

        // removed original final equation
        instr[final_eqn_idx+3] = 0;

        // remove extra zeroes
        WORD *t = &instr[0];

        ebegin = &*instr.begin();
        eend = ebegin+instr.size();
        int de=0;

        for (WORD *e=ebegin; e!=eend; e+=de) {
                de = *(e+2);
                int n=3;
                while (*(e+n) != 0) n+=*(e+n);
                n++;
                memmove (t, e, n*sizeof(WORD));
                *(t+2) = n;
                t += n;
        }

        instr.resize(t - &instr[0]);

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: optimize_greedy(numoper=%d) : numoper=%d",
                                                thetime_str().c_str(), old_num_oper, count_operators(instr));
#endif

        return instr;
}

/*
          #] optimize_greedy : 
          #[ recycle_variables :
*/

/**  Recycle variables
 *
 *         Description
 *         ===========
 *         The current input uses many temporary variables. Many of them
 *         become obsolete at some point during the evaluation of the code,
 *         so can be recycled. This method renumbers the temporary
 *         variables, so that they are recycled. Furthermore, the input is
 *         order in depth-first order, so that the instructions can be
 *         performed consecutively.
 *
 *         Implementation details
 *         ======================
 *         First, for each subDAG, an estimate for the number of variables
 *         needed is made. This is done by the following recursive formula:
 *
 *                #vars(x) = max(#vars(ch_i(x)) + i),
 *
 *         with ch_i(x) the i-th child of x, where the childs are ordered
 *         w.r.t. #vars(ch_i). This formula is exact if the input forms a
 *         tree, and otherwise gives a reasonable estimate.
 *
 *         Then, the instructions are reordered in a depth-first order with
 *         childs ordered w.r.t. #vars. Next, the times that variables
 *         become obsolete are found. Each LHS of an instruction is
 *         renumbered to the lowest-numbered temporary variable that is
 *         available at that time.
 */
vector<WORD> recycle_variables (const vector<WORD> &all_instr) {

#ifdef DEBUG_MORE
        MesPrint ("*** [%s, w=%w] CALL: recycle_variables", thetime_str().c_str());
#endif

        // get starting positions of instructions
        vector<const WORD *> instr;
        const WORD *tbegin = &*all_instr.begin();
        const WORD *tend = tbegin+all_instr.size();
        for (const WORD *t=tbegin; t!=tend; t+=*(t+2))
                instr.push_back(t);
        int n = instr.size();

        // determine with expressions are connected, how many intermediate
        // are needed (assuming it's a expression tree instead of a DAG) and
        // sort the leaves such that you need a minimal number of variables
        vector<int> vars_needed(n);
        vector<bool> vis(n,false);
        vector<vector<int> > conn(n);

        stack<int> s;
        s.push(n);

        while (!s.empty()) {
                int i=s.top(); s.pop();
                if (i>0) {
                        i--;
                        if (vis[i]) continue;
                        vis[i]=true;
                        s.push(-(i+1));

                        // find all connections
                        for (const WORD *t=instr[i]+3; *t!=0; t+=*t)
                                if (*t!=1+ABS(*(t+*t-1)) && *(t+1)==EXTRASYMBOL) {
                                        int k = *(t+3);
                                        conn[i].push_back(k);
                                        s.push(k+1);
                                }
                }
                else {
                        i=-i-1;

                        // sort the childs w.r.t. needed variables
                        vector<pair<int,int> > need;
                        for (int j=0; j<(int)conn[i].size(); j++)
                                need.push_back(make_pair(vars_needed[conn[i][j]], conn[i][j]));

                        // keep the comma expression in proper order
                        if (*(instr[i]+1) != OPER_COMMA)
                                sort(need.rbegin(), need.rend());

                        vars_needed[i] = 1;
                        for (int j=0; j<(int)need.size(); j++) {
                                vars_needed[i] = max(vars_needed[i], need[j].first+j);
                                conn[i][j] = need[j].second;
                        }
                }
        }

        // order the instructions in depth-first order and determine the first
        // and last occurrences of variables
        vector<int> order, first(n,0), last(n,0);
        vis = vector<bool>(n,false);
        s.push(n);

        while (!s.empty()) {

                int i=s.top(); s.pop();

                if (i>0) {
                        i--;
                        if (vis[i]) continue;
                        vis[i]=true;
                        s.push(-(i+1));
                        for (int j=(int)conn[i].size()-1; j>=0; j--)
                                s.push(conn[i][j]+1);
                }
                else {
                        i=-i-1;
                        first[i] = last[i] = order.size();
                        order.push_back(i);
                        for (int j=0; j<(int)conn[i].size(); j++) {
                                int k = conn[i][j];
                                last[k] = max(last[k], first[i]);
                        }
                }
        }

        // find the renumbering to recycled variables, where at any time the
        // lowest-indexed variable that can be used is chosen
        int numvar=0;
        set<int> var;
        vector<int> renum(n);

        for (int i=0; i<(int)order.size(); i++) {
                for (int j=0; j<(int)conn[order[i]].size(); j++) {
                        int k = conn[order[i]][j];
                        if (last[k] == i) var.insert(renum[k]);
                }

                if (var.empty()) var.insert(numvar++);
                renum[order[i]] = *var.begin(); var.erase(var.begin());
        }

        // put the number of variables used in a preprocessor variable

        // generate new instructions with the renumbering
        vector<WORD> newinstr;

        for (int i=0; i<(int)order.size(); i++) {
                int x = order[i];
                int j = newinstr.size();
                newinstr.insert(newinstr.end(), instr[x], instr[x]+*(instr[x]+2));

                newinstr[j] = renum[newinstr[j]];

                for (WORD *t=&newinstr[j+3]; *t!=0; t+=*t)
                        if (*t!=1+ABS(*(t+*t-1)) && *(t+1)==EXTRASYMBOL)
                                *(t+3) = renum[*(t+3)];
        }

#ifdef DEBUG_MORE
        MesPrint ("*** [%s, w=%w] DONE: recycle_variables", thetime_str().c_str());
#endif

        return newinstr;
}

/*
          #] recycle_variables : 
          #[ optimize_expression_given_Horner :
*/

/**  Optimize expression given a Horner scheme
 *
 *         Description
 *         ===========
 *         This method picks one Horner scheme from the list of best Horner
 *         schemes, applies this scheme to the expression and then,
 *         depending on optimize.settings, does a common subexpression
 *         elimination (CSE) or performs greedy optimizations.
 *
 *         CSE is fast, while greedy might be slow. CSE followed by greedy
 *         is faster than greedy alone, but typically results in slightly
 *         worse code (not proven; just observed).
 */
void optimize_expression_given_Horner () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: optimize_expression_given_Horner", thetime_str().c_str());
#endif

        GETIDENTITY;

        // initialize timer
        LONG start_time = TimeWallClock(1);
        LONG time_limit = 100 * AO.Optimize.greedytimelimit / (AO.Optimize.horner == O_MCTS ? AO.Optimize.mctsnumkeep : 1);
        if (time_limit == 0) time_limit=MAXPOSITIVE;

        // pick a Horner scheme from the list
        LOCK(optimize_lock);
        vector<WORD> Horner_scheme = optimize_best_Horner_schemes.back();
        optimize_best_Horner_schemes.pop_back();
        UNLOCK(optimize_lock);

//        if ( ( AO.Optimize.debugflags&2 ) == 2 ) {
//                MesPrint ("Scheme: %a",Horner_scheme.size(),&(Horner_scheme[0]));
//        }

        // apply Horner scheme
        vector<WORD> tree = Horner_tree(optimize_expr, Horner_scheme);

        // generate instructions, eventually with CSE
        vector<WORD> instr;

        if (AO.Optimize.method == O_CSE || AO.Optimize.method == O_CSEGREEDY)
                instr = generate_instructions(tree, true);
        else
                instr = generate_instructions(tree, false);
        /// eventually do greedy optimizations
        if (AO.Optimize.method == O_CSEGREEDY || AO.Optimize.method == O_GREEDY) {
                instr = merge_operators(instr, false);
                instr = optimize_greedy(instr, time_limit-(TimeWallClock(1)-start_time));
                instr = merge_operators(instr, true);
                instr = optimize_greedy(instr, time_limit-(TimeWallClock(1)-start_time));
        }
        instr = merge_operators(instr, true);

        // recycle the temporary variables
        instr = recycle_variables(instr);

        // determine the quality of the code and possibly update the best code
        int num_oper = count_operators(instr);

        LOCK(optimize_lock);
        if (num_oper < optimize_best_num_oper) {
                optimize_num_vars = Horner_scheme.size();
                optimize_best_num_oper = num_oper;
                optimize_best_instr = instr;
                optimize_best_vars = vector<WORD>(AN.poly_vars, AN.poly_vars+AN.poly_num_vars);
        }
        UNLOCK(optimize_lock);

  // clean poly_vars, that are allocated by Horner_tree
        AN.poly_num_vars = 0;
        M_free(AN.poly_vars,"poly_vars");

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: optimize_expression_given_Horner", thetime_str().c_str());
#endif
}

/*
          #] optimize_expression_given_Horner : 
          #[ PF_optimize_expression_given_Horner :
*/
#ifdef WITHMPI

// Initialization.
void PF_optimize_expression_given_Horner_master_init () {
        // Nothing to do for now.
}

// Wait for an idle slave and return the process number.
int PF_optimize_expression_given_Horner_master_next() {

        // Find an idle slave.
        int next;
        PF_LongSingleReceive(PF_ANY_SOURCE, PF_OPT_HORNER_MSGTAG, &next, NULL);
        return next;

}

// The main function on the master.
void PF_optimize_expression_given_Horner_master () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_optimize_expression_given_Horner_master", thetime_str().c_str());
#endif

        // pick a Horner scheme from the list
        vector<WORD> Horner_scheme = optimize_best_Horner_schemes.back();
        optimize_best_Horner_schemes.pop_back();

        // Find an idle slave.
        int next = PF_optimize_expression_given_Horner_master_next();

        // Send a new task to the slave.
        PF_PrepareLongSinglePack();
        int len = Horner_scheme.size();
        PF_LongSinglePack(&len, 1, PF_INT);
        PF_LongSinglePack(&Horner_scheme[0], len, PF_WORD);
        PF_LongSingleSend(next, PF_OPT_HORNER_MSGTAG);

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_optimize_expression_given_Horner_master", thetime_str().c_str());
#endif

}

// Wait for all the slaves to finish their tasks.
void PF_optimize_expression_given_Horner_master_wait () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_optimize_expression_given_Horner_master_wait", thetime_str().c_str());
#endif

        // Wait for all the slaves.
        for (int i = 1; i < PF.numtasks; i++) {
                int next = PF_optimize_expression_given_Horner_master_next();
                // Send a null task.
                PF_PrepareLongSinglePack();
                int len = 0;
                PF_LongSinglePack(&len, 1, PF_INT);
                PF_LongSingleSend(next, PF_OPT_HORNER_MSGTAG);
        }

        // Combine the result from all the slaves.
        optimize_best_num_oper = INT_MAX;
        for (int i = 1; i < PF.numtasks; i++) {
                PF_LongSingleReceive(PF_ANY_SOURCE, PF_OPT_COLLECT_MSGTAG, NULL, NULL);

                int len;

                // The first integer is the number of operations.
                PF_LongSingleUnpack(&len, 1, PF_INT);

                if (len >= optimize_best_num_oper) continue;

                // Update the best result.
                optimize_best_num_oper = len;
                PF_LongSingleUnpack(&len, 1, PF_INT);
                optimize_best_instr.resize(len);
                PF_LongSingleUnpack(&optimize_best_instr[0], len, PF_WORD);
                PF_LongSingleUnpack(&len, 1, PF_INT);
                optimize_best_vars.resize(len);
                PF_LongSingleUnpack(&optimize_best_vars[0], len, PF_WORD);

                optimize_num_vars = optimize_best_vars.size();        // TODO
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_optimize_expression_given_Horner_master_wait", thetime_str().c_str());
#endif

}

// The main function on the slaves.
void PF_optimize_expression_given_Horner_slave () {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: PF_optimize_expression_given_Horner_slave", thetime_str().c_str());
#endif

        optimize_best_Horner_schemes.clear();
        optimize_best_num_oper = INT_MAX;

        int dummy = 0;
        int len;

        for (;;) {
                // Ask the master for the next task.
                PF_PrepareLongSinglePack();
                PF_LongSinglePack(&dummy, 1, PF_INT);
                PF_LongSingleSend(MASTER, PF_OPT_HORNER_MSGTAG);

                // Get a task from the master.
                PF_LongSingleReceive(MASTER, PF_OPT_HORNER_MSGTAG, NULL, NULL);

                // Length of the task.
                PF_LongSingleUnpack(&len, 1, PF_INT);

                // No task remains.
                if (len == 0) break;

                // Perform the given task.
                optimize_best_Horner_schemes.push_back(vector<WORD>());
                vector<WORD> &Horner_scheme = optimize_best_Horner_schemes.back();
                Horner_scheme.resize(len);
                PF_LongSingleUnpack(&Horner_scheme[0], len, PF_WORD);
                optimize_expression_given_Horner ();
        }

        // Send the result to the master.
        PF_PrepareLongSinglePack();
        PF_LongSinglePack(&optimize_best_num_oper, 1, PF_INT);
        if (optimize_best_num_oper != INT_MAX) {
                len = optimize_best_instr.size();
                PF_LongSinglePack(&len, 1, PF_INT);
                PF_LongSinglePack(&optimize_best_instr[0], len, PF_WORD);
                len = optimize_best_vars.size();
                PF_LongSinglePack(&len, 1, PF_INT);
                PF_LongSinglePack(&optimize_best_vars[0], len, PF_WORD);
        }
        PF_LongSingleSend(MASTER, PF_OPT_COLLECT_MSGTAG);

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: PF_optimize_expression_given_Horner_slave", thetime_str().c_str());
#endif

}

#endif
/*
          #] PF_optimize_expression_given_Horner : 
          #[ generate_output :
*/

/**  Generate output
 *
 *         Description
 *         ===========
 *         This method prepares the instructions for printing. They are
 *         stored in Form format, so that they can be printed by
 *         "PrintExtraSymbol". The results are stored in the buffer
 *         AO.OptimizeResult.
 */
VOID generate_output (const vector<WORD> &instr, int exprnr, int extraoffset, const vector<vector<WORD> > &brackets) {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: generate_output", thetime_str().c_str());
#endif

        GETIDENTITY;
        vector<WORD> output;

        // one-indexed instead of zero-indexed
        extraoffset++;
        int num = 0;
        int maxsize = (int)instr.size();
        for (int i=0; i<maxsize; i+=instr[i+2]) num++;
        int *tstart = (int *)Malloc1(num*sizeof(int),"nplaces");
        num = 0;
        for (int i=0; i<maxsize; i+=instr[i+2]) tstart[num++] = i;
        for (int j=0; j<num; j++) {
                int i = tstart[j];

                // copy arguments
                WCOPY(AT.WorkPointer, &instr[i+3], (instr[i+2]-3));

                // update maximal variable number
                AO.OptimizeResult.maxvar = MaX(AO.OptimizeResult.maxvar, instr[i]+extraoffset);

                // renumber symbols and extrasymbols to correct values
                for (WORD *t=AT.WorkPointer; *t!=0; t+=*t) {
                        if (*t == ABS(*(t+*t-1))+1) continue;
                        if (*(t+1)==SYMBOL) *(t+3) = optimize_best_vars[*(t+3)];
                        if (*(t+1)==EXTRASYMBOL) {
                                *(t+1) = SYMBOL;
                                *(t+3) = MAXVARIABLES-*(t+3)-extraoffset;
                        }
                }

                // reformat multiplication instructions, since their current
                // format is "expr.nr OPER_MUL length arguments", which differs
                // from Form's format for a product of symbols
                if (instr[i+1]==OPER_MUL) {

                        WORD *now=AT.WorkPointer+1;
                        int dt;
                        bool coeff=false;
                        int coeff_sign=1;

                        for (WORD *t=AT.WorkPointer; *t!=0; t+=dt) {
                                dt = *t;
                                if (*t == ABS(*(t+*t-1))+1) {
                                        // copy coefficient
                                        memmove(AT.WorkPointer+instr[i+2], t, *t*sizeof(WORD));
                                        coeff = true;
                                }
                                else {
                                        // move symbol
                                        int n = *(t+2);
                                        memmove(now, t+1, n*sizeof(WORD));
                                        now += n;
                                        coeff_sign *= SGN(*(t+dt-1));
                                }
                        }

                        if (coeff) {
                                // add existing coefficient
                                int n = *(AT.WorkPointer + instr[i+2]) - 1;
                                memmove(now, AT.WorkPointer+instr[i+2]+1, n*sizeof(WORD));
                                now += n;
                        }
                        else {
                                // add coefficient of one
                                *now++=1;
                                *now++=1;
                                *now++=3;
                        }

                        *(now-1) *= coeff_sign;

                        *AT.WorkPointer = now - AT.WorkPointer;
                        *now++ = 0;
                }

                // in the case of simultaneous optimization of expressions, add the
                // brackets to the final expression
                if (instr[i+1]==OPER_COMMA) {
                        WORD *start = AT.WorkPointer + instr[i+2];
                        WORD *now = start;
                        int b=0;
                        for (const WORD *t=AT.WorkPointer; *t!=0; t+=*t) {
                                if ( ( brackets[b].size() != 0 ) && ( brackets[b][0] == 0 ) ) break;
                                *now++ = *t + brackets[b].size();
                                if (brackets[b].size() != 0) {
                                        memcpy(now, &brackets[b][0], brackets[b].size()*sizeof(WORD));
                                        now += brackets[b].size();
                                }
                                memcpy(now, t+1, (*t-1)*sizeof(WORD));
                                now += *t-1;
                                b++;
                        }
                        *now++ = 0;
                        memmove(AT.WorkPointer, start, (now-start)*sizeof(WORD));
                }

                // add the number of the extra symbol; if it is the last one,
                // replace the extra symbol number with the expression number
                if (i+instr[i+2]<(int)instr.size())
                        output.push_back(instr[i]+extraoffset);
                else {
                        output.push_back(-(exprnr+1));
                }

                // add code for this symbol
                int n=0;
                while (*(AT.WorkPointer+n)!=0)
                        n += *(AT.WorkPointer+n);
                n++;
                output.insert(output.end(), AT.WorkPointer, AT.WorkPointer+n);
        }

        // trailing zero
        output.push_back(0);

        // clear buffer
        if (AO.OptimizeResult.code != NULL)
                M_free(AO.OptimizeResult.code, "optimize output");

        M_free(tstart,"nplaces");

        // copy buffer
        AO.OptimizeResult.codesize = output.size();
        AO.OptimizeResult.code = (WORD *)Malloc1(output.size()*sizeof(WORD), "optimize output");
        memcpy(AO.OptimizeResult.code, &output[0], output.size()*sizeof(WORD));

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: generate_output", thetime_str().c_str());
#endif
}

/*
          #] generate_output : 
          #[ generate_expression :
*/

/**  Generate expression
 *
 *         Description
 *         ===========
 *         This method modifies the original optimized expression by an
 *         expression with extra symbols. This is used for "#Optimize".
 */
WORD generate_expression (WORD exprnr) {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: generate_expression", thetime_str().c_str());
#endif

        GETIDENTITY;

        WORD *oldWorkPointer = AT.WorkPointer;

        CBUF *C = cbuf+AC.cbufnum;
        WORD *term = AT.WorkPointer, oldcurexpr = AR.CurExpr;
        POSITION position;
        EXPRESSIONS e = Expressions+exprnr;
        SetScratch(AR.infile,&(e->onfile));

        if ( GetTerm(BHEAD term) <= 0 ) {
                MesPrint("Expression %d has problems in scratchfile",exprnr);
                Terminate(-1);
        }
        SeekScratch(AR.outfile,&position);
        e->onfile = position;
        if ( PutOut(BHEAD term,&position,AR.outfile,0) < 0 ) {
                MesPrint("Expression %d has problems in output scratchfile",exprnr);
                Terminate(-1);
        }

        AR.CurExpr = exprnr;
        NewSort(BHEAD0);

        // scan for the original expression (marked by *t<0) and give the
        // terms to Generator
        WORD *t = AO.OptimizeResult.code;
        {
                WORD old = AR.Cnumlhs; AR.Cnumlhs = 0;
                // We can use the remaining part of the WorkSpace for every term that
                // goes through Generator. We have WorkSpace problems here if we use
                // the (modified) AT.WorkPointer every time in the loop.
                WORD* currentWorkPointer = AT.WorkPointer;
                while (*t!=0) {
                        bool is_expr = *t < 0;
                        t++;
                        while (*t!=0) {
                                if (is_expr) {
                                        memcpy(currentWorkPointer, t, *t*sizeof(WORD));
                                        Generator(BHEAD currentWorkPointer, C->numlhs);
                                }
                                t+=*t;
                        }
                        t++;
                }
                AR.Cnumlhs = old;
        }

        // final sorting
        if (EndSort(BHEAD NULL,0) < 0) {
                LowerSortLevel();
                Terminate(-1);
        }

        AT.WorkPointer = oldWorkPointer;
        AR.CurExpr = oldcurexpr;

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: generate_expression", thetime_str().c_str());
#endif

        return 0;
}

/*
          #] generate_expression : 
          #[ optimize_print_code :
*/

/**  Print optimized code
 *
 *         Description
 *         ===========
 *         This method prints the optimized code via
 *         "PrintExtraSymbol". Depending on the flag, the original
 *         expression is printed (for "Print") or not (for "#Optimize /
 *         #write "%O").
 */
VOID optimize_print_code (int print_expr) {

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: optimize_print_code", thetime_str().c_str());
#endif
        if ( ( AO.Optimize.debugflags & 1 ) != 0 ) {
/*
 *                The next code is for debugging purposes. We may want the statements
 *                in reverse order to substitute them all back.
 *                Jan used a Mathematica program to do this. Here we make that
 *                        Format Ox,debugflag=1;
 *                Creates reverse order during printing.
 *                All we have to do is put id in front of the statements. This is done
 *                in PrintExtraSymbol.
 */
                WORD *t = AO.OptimizeResult.code;
                WORD num = 0;        // First we count the number of objects.
                while (*t!=0) {
                        num++;
                        t++; while (*t!=0) t+=*t; t++;
                }
                WORD **tstart = (WORD **)Malloc1(num*sizeof(WORD *),"print objects");
                t = AO.OptimizeResult.code; num = 0; // Now we get the addresses
                while (*t!=0) {
                        tstart[num++] = t;
                        t++; while (*t!=0) t+=*t; t++;
                }
                // Flip the addresses
                int halfnum = num/2;
                for (int i=0; i<halfnum; i++) { swap(tstart[i],tstart[num-1-i]); }
                for ( int i = 0; i < num; i++ ) {
                        t = tstart[i];
                        if (*t > 0)
                                PrintExtraSymbol(*t, t+1, EXTRASYMBOL);
                        else if (print_expr)
                                PrintExtraSymbol(-*t-1, t+1, EXPRESSIONNUMBER);
                }
                CBUF *C = cbuf + AM.sbufnum;
                if (C->numrhs >= AO.OptimizeResult.minvar)
                        PrintSubtermList(AO.OptimizeResult.minvar, C->numrhs);
        }
        else {
                // print extra symbols from ConvertToPoly in optimization
                CBUF *C = cbuf + AM.sbufnum;
                if (C->numrhs >= AO.OptimizeResult.minvar)
                        PrintSubtermList(AO.OptimizeResult.minvar, C->numrhs);

                WORD *t = AO.OptimizeResult.code;
                while (*t!=0) {
                        if (*t > 0) {
                                PrintExtraSymbol(*t, t+1, EXTRASYMBOL);
                        }
                        else if (print_expr)
                                PrintExtraSymbol(-*t-1, t+1, EXPRESSIONNUMBER);
                        t++;
                        while (*t!=0) t+=*t;
                        t++;
                }
        }

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: optimize_print_code", thetime_str().c_str());
#endif
}

/*
          #] optimize_print_code : 
          #[ Optimize :
*/

/**  Optimization of expression
 *
 *         Description
 *         ===========
 *         This method takes an input expression and generates optimized
 *         code to calculate its value. The following methods are called to
 *         do so:
 *
 *         (1) get_expression : to read to expression
 *
 *         (2) get_brackets : find brackets for simultaneous optimization
 *
 *         (3) occurrence_order or find_Horner_MCTS : to determine (the)
 *         Horner scheme(s) to use; this depends on AO.optimize.horner
 *
 *         (4) optimize_expression_given_Horner : to do the optimizations
 *         for each Horner scheme; this method does either CSE or greedy
 *         optimizations dependings on AO.optimize.method
 *
 *         (5) generate_output : to format the output in Form notation and
 *         store it in a buffer
 *
 *         (6a) optimize_print_code : to print the expression (for "Print")
 *         or
 *         (6b) generate_expression : to modify the expression (for
 *         "#Optimize")
 *
 *         On ParFORM, all the processes must call this function at the same
 *         time. Then
 *
 *         (1) Because only the master can access to the expression to be
 *         optimized, the master broadcast the expression to all the slaves
 *         after reading the expression (PF_get_expression).
 *
 *         (2) get_brackets reads optimize_expr as the input and it works
 *         also on the slaves. We leave it although the bracket information
 *         is not needed on the slaves (used in (5) on the master).
 *
 *         (3) and (4) find_Horner_MCTS and optimize_expression_given_Horner
 *         are parallelized.
 *
 *         (5), (6a) and (6b) are needed only on the master.
 */
int Optimize (WORD exprnr, int do_print) {

#if defined(WITHMPI) && (defined(DEBUG) || defined(DEBUG_MORE) || defined(DEBUG_MCTS) || defined(DEBUG_GREEDY))
        // set AS.printflag negative temporary.
        struct save_printflag {
                save_printflag() {
                        oldprintflag = AS.printflag;
                        AS.printflag = -1;
                }
                ~save_printflag() {
                        AS.printflag = oldprintflag;
                }
                int oldprintflag;
        } save_printflag_;
#endif

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] CALL: Optimize", thetime_str().c_str());
        MesPrint ("*** %"); PrintRunningTime();
#endif

#ifdef WITHPTHREADS
        optimize_lock = dummylock;
#endif

        AO.OptimizeResult.minvar = (cbuf + AM.sbufnum)->numrhs + 1;

        if (get_expression(exprnr) < 0) return -1;
        vector<vector<WORD> > brackets = get_brackets();

#ifdef DEBUG
#ifdef WITHMPI
                if (PF.me == MASTER)
#endif
                MesPrint ("*** runtime after preparing the expression: %"); PrintRunningTime();
#endif

        if (optimize_expr[0]==0 ||
                        (optimize_expr[optimize_expr[0]]==0 && optimize_expr[0]==ABS(optimize_expr[optimize_expr[0]-1])+1) ||
                        (optimize_expr[optimize_expr[0]]==0 && optimize_expr[0]==8 &&
                         optimize_expr[5]==1 && optimize_expr[6]==1 && ABS(optimize_expr[7])==3)) {
                if (AO.OptimizeResult.code != NULL)
                        M_free(AO.OptimizeResult.code, "optimize output");

                // zero terms or one trivial term (number or +/-variable), so no
                // optimization, so copy expression; special case because without
                // operators the optimization crashes
                AO.OptimizeResult.code = (WORD *)Malloc1((optimize_expr[0]+3)*sizeof(WORD), "optimize output");
                AO.OptimizeResult.code[0] = -(exprnr+1);
                memcpy(AO.OptimizeResult.code+1, optimize_expr, (optimize_expr[0]+1)*sizeof(WORD));
                AO.OptimizeResult.code[optimize_expr[0]+2] = 0;
        }
        else {
                // find Horner scheme(s)
                optimize_best_Horner_schemes.clear();
                if (AO.Optimize.horner == O_OCCURRENCE) {
                        if (AO.Optimize.hornerdirection != O_BACKWARD)
                                optimize_best_Horner_schemes.push_back(occurrence_order(optimize_expr, false));
                        if (AO.Optimize.hornerdirection != O_FORWARD)
                                optimize_best_Horner_schemes.push_back(occurrence_order(optimize_expr, true));
                }
                else {
                        if (AO.Optimize.horner == O_SIMULATED_ANNEALING) {
                                optimize_best_Horner_schemes.push_back(simulated_annealing());
                        }
                        else {
                                mcts_best_schemes.clear();
                                if ( AO.inscheme ) {
                                        optimize_best_Horner_schemes.push_back(vector<WORD>(AO.schemenum));
                                        for ( int i=0; i < AO.schemenum; i++ )
                                                optimize_best_Horner_schemes[0][i] = AO.inscheme[i];
                                }
                                else {
                                        for ( int i = 0; i < AO.Optimize.mctsnumrepeat; i++ )
                                                find_Horner_MCTS();
                                        // generate results
                                        for (set<pair<int,vector<WORD> > >::iterator i=mcts_best_schemes.begin(); i!=mcts_best_schemes.end(); i++) {
                                                optimize_best_Horner_schemes.push_back(i->second);
#ifdef DEBUG_MCTS
                                                MesPrint ("{%a} -> %d",i->second.size(), &i->second[0], i->first);
#endif
                                        }
                                }
                                // clear the tree by making a new empty one.
                                mcts_root = tree_node();
                        }
                }
#ifdef DEBUG
#ifdef WITHMPI
                if (PF.me == MASTER)
#endif
                MesPrint ("*** runtime after Horner: %"); PrintRunningTime();
#endif

#ifdef WITHMPI
                if (PF.me == MASTER ) {
                        PF_optimize_expression_given_Horner_master_init();
#endif

                // find best Horner scheme and results
                optimize_best_num_oper = INT_MAX;

                int imax = (int)optimize_best_Horner_schemes.size();

                for (int i=0; i<imax; i++) {
#if defined(WITHPTHREADS)
                        if (AM.totalnumberofthreads > 1)
                                optimize_expression_given_Horner_threaded();
                        else
#elif defined(WITHMPI)
                        if (PF.numtasks > 1)
                                PF_optimize_expression_given_Horner_master();
                        else
#endif
                                optimize_expression_given_Horner();
                }

#ifdef WITHMPI
                        PF_optimize_expression_given_Horner_master_wait();
                }
                else {
                        if (PF.numtasks > 1)
                                PF_optimize_expression_given_Horner_slave();
                }
#endif

#ifdef WITHPTHREADS
                MasterWaitAll();
#endif
                // format results, then print it (for "Print") or modify
                // expression (for "#Optimize")
#ifdef WITHMPI
                if (PF.me == MASTER)
#endif
                generate_output(optimize_best_instr, exprnr, cbuf[AM.sbufnum].numrhs, brackets);
#ifdef WITHMPI
                else {
                        // non-null dummy code for slaves
                        if (AO.OptimizeResult.code == NULL) {
                                AO.OptimizeResult.code = (WORD *)Malloc1(sizeof(WORD), "optimize output");
                        }
                }
#endif
        }

#ifdef WITHMPI
        if (PF.me == MASTER) {
                PF_PreparePack();
                PF_Pack(&AO.OptimizeResult.minvar, 1, PF_WORD);
                PF_Pack(&AO.OptimizeResult.maxvar, 1, PF_WORD);
        }
        PF_Broadcast();
        if (PF.me != MASTER) {
                PF_Unpack(&AO.OptimizeResult.minvar, 1, PF_WORD);
                PF_Unpack(&AO.OptimizeResult.maxvar, 1, PF_WORD);
        }
#endif

        // set preprocessor variables
        char str[100];
        snprintf (str,100,"%d",AO.OptimizeResult.minvar);
        PutPreVar((UBYTE *)"optimminvar_",(UBYTE *)str,0,1);
        snprintf (str,100,"%d",AO.OptimizeResult.maxvar);
        PutPreVar((UBYTE *)"optimmaxvar_",(UBYTE *)str,0,1);

        if (do_print) {
#ifdef WITHMPI
                if (PF.me == MASTER)
#endif
                optimize_print_code(1);
                ClearOptimize();
        }
        else {
#ifdef WITHMPI
                if (PF.me == MASTER)
#endif
                generate_expression(exprnr);
        }

#ifdef WITHMPI
        if (PF.me == MASTER) {
#endif

        if ( AO.Optimize.printstats > 0 ) {
                char str[20];
                MesPrint("");
                count_operators(optimize_expr,true);
                int numop = count_operators(optimize_best_instr,true);
                snprintf(str,20,"%d",numop);
                PutPreVar((UBYTE *)"optimvalue_",(UBYTE *)str,0,1);
        }
        else {
                char str[20];
                int numop = count_operators(optimize_best_instr,false);
                snprintf(str,20,"%d",numop);
                PutPreVar((UBYTE *)"optimvalue_",(UBYTE *)str,0,1);
        }

        if ( ( AO.Optimize.schemeflags&1 ) == 1 ) {
                GETIDENTITY
                UBYTE *OutScr, *Out, *old1 = AO.OutputLine, *old2 = AO.OutFill;
                int i, sym;
                AO.OutputLine = AO.OutFill = (UBYTE *)AT.WorkPointer;
                FiniLine();
                OutScr = (UBYTE *)AT.WorkPointer + ( TOLONG(AT.WorkTop) - TOLONG(AT.WorkPointer) ) /2;
                TokenToLine((UBYTE *)" Scheme selected: ");
                for ( i = 0; i < optimize_num_vars; i++ ) {
                        Out = OutScr;
                        sym = optimize_best_vars[i];
                        if ( i > 0 ) TokenToLine((UBYTE *)",");
                        if ( sym < NumSymbols ) {
                                StrCopy(FindSymbol(sym),OutScr);
/*                                StrCopy(VARNAME(symbols,sym),OutScr); */
                        }
                        else {
                                Out = StrCopy((UBYTE *)AC.extrasym,Out);
                                if ( AC.extrasymbols == 0 ) {
                                        Out = NumCopy((MAXVARIABLES-sym),Out);
                                        Out = StrCopy((UBYTE *)"_",Out);
                                }
                                else if ( AC.extrasymbols == 1 ) {
                                        Out = AddArrayIndex((MAXVARIABLES-sym),Out);
                                }
                        }
                        TokenToLine(OutScr);
                }
                TokenToLine((UBYTE *)";");
                FiniLine();
                AO.OutFill = old2;
                AO.OutputLine = old1;
        }

        {
                GETIDENTITY
                UBYTE *OutScr, *Out, *outstring = 0;
                int i, sym;
                AO.OutputLine = AO.OutFill = (UBYTE *)AT.WorkPointer;
                OutScr = (UBYTE *)AT.WorkPointer + ( TOLONG(AT.WorkTop) - TOLONG(AT.WorkPointer) ) /2;
                for ( i = 0; i < optimize_num_vars; i++ ) {
                        Out = OutScr;
                        sym = optimize_best_vars[i];
                        if ( sym < NumSymbols ) {
                                StrCopy(FindSymbol(sym),OutScr);
/*                                StrCopy(VARNAME(symbols,sym),OutScr); */
                        }
                        else {
                                Out = StrCopy((UBYTE *)AC.extrasym,Out);
                                if ( AC.extrasymbols == 0 ) {
                                        Out = NumCopy((MAXVARIABLES-sym),Out);
                                        Out = StrCopy((UBYTE *)"_",Out);
                                }
                                else if ( AC.extrasymbols == 1 ) {
                                        Out = AddArrayIndex((MAXVARIABLES-sym),Out);
                                }
                        }
                        outstring = AddToString(outstring,OutScr,1);
                }
                if ( outstring == 0 ) {
                        PutPreVar((UBYTE *)"optimscheme_",(UBYTE *)"",0,1);
                }
                else {
                        PutPreVar((UBYTE *)"optimscheme_",(UBYTE *)outstring,0,1);
                        M_free(outstring,"AddToString");
                }
        }

#ifdef WITHMPI
        }

        // synchronize optimvalue_ and optimscheme_
        if ( PF.me == MASTER ) {
                UBYTE *value;
                int bytes;

                PF_PrepareLongMultiPack();

                value = GetPreVar((UBYTE *)"optimvalue_", WITHERROR);
                bytes = strlen((char *)value);
                PF_LongMultiPack(&bytes, 1, PF_INT);
                PF_LongMultiPack(value, bytes, PF_BYTE);

                value = GetPreVar((UBYTE *)"optimscheme_", WITHERROR);
                bytes = strlen((char *)value);
                PF_LongMultiPack(&bytes, 1, PF_INT);
                PF_LongMultiPack(value, bytes, PF_BYTE);
        }
        PF_LongMultiBroadcast();
        if ( PF.me != MASTER ) {
                static vector<UBYTE> prevarbuf;
                UBYTE *value;
                int bytes;

                PF_LongMultiUnpack(&bytes, 1, PF_INT);
                prevarbuf.reserve(bytes + 1);
                value = &prevarbuf[0];
                PF_LongMultiUnpack(value, bytes, PF_BYTE);
                value[bytes] = '\0';  // null terminator
                PutPreVar((UBYTE *)"optimvalue_", value, NULL, 1);

                PF_LongMultiUnpack(&bytes, 1, PF_INT);
                prevarbuf.reserve(bytes + 1);
                value = &prevarbuf[0];
                PF_LongMultiUnpack(value, bytes, PF_BYTE);
                value[bytes] = '\0';  // null terminator
                PutPreVar((UBYTE *)"optimscheme_", value, NULL, 1);
        }
#endif

        // cleanup
        M_free(optimize_expr,"LoadOptim");

#ifdef DEBUG
        MesPrint ("*** [%s, w=%w] DONE: Optimize", thetime_str().c_str());
#endif

        return 0;
}

/*
          #] Optimize : 
          #[ ClearOptimize :
*/

/**  Optimization of expression
 *
 *         Description
 *         ===========
 *         Clears the buffers that were used for optimization output.
 *         Clears the expression from the buffers (marks it to be dropped).
 *         Note: we need to use the expression by its name, because numbers
 *         may change if we drop other expressions between when we do the
 *         optimizations and clear the results (in execute.c). Also this is
 *         not 100% safe, because we could overwrite the optimized
 *         expression. But that can be done only in a Local or Global
 *         statement and hence we only have to test there that we might have
 *         to call ClearOptimize first. (in file comexpr.c)
 */
int ClearOptimize()
{
        char str[20];
        WORD numexpr, *w;
        int error = 0;
        if ( AO.OptimizeResult.code != NULL ) {
                M_free(AO.OptimizeResult.code, "optimize output");
                AO.OptimizeResult.code = NULL;
                AO.OptimizeResult.codesize = 0;
                PutPreVar((UBYTE *)"optimminvar_",(UBYTE *)("0"),0,1);
                PutPreVar((UBYTE *)"optimmaxvar_",(UBYTE *)("0"),0,1);
                PruneExtraSymbols(AO.OptimizeResult.minvar-1);
                cbuf[AM.sbufnum].numrhs = AO.OptimizeResult.minvar-1;
                AO.OptimizeResult.minvar = AO.OptimizeResult.maxvar = 0;
        }
        if ( AO.OptimizeResult.nameofexpr != NULL ) {
/*
                We have to pick up the expression by its name. Numbers may change.
                Note that this requires that when we overwrite an expression, we
                check that it is not an optimized expression. See execute.c and
                comexpr.c
*/
                if ( GetName(AC.exprnames,AO.OptimizeResult.nameofexpr,&numexpr,NOAUTO) != CEXPRESSION ) {
                        MesPrint("@Internal error while clearing optimized expression %s ",AO.OptimizeResult.nameofexpr);
                        Terminate(-1);
                }
                M_free(AO.OptimizeResult.nameofexpr, "optimize expression name");
                AO.OptimizeResult.nameofexpr = NULL;
                w = &(Expressions[numexpr].status);
                *w = SetExprCases(DROP,1,*w);
                if ( *w < 0 ) error = 1;
        }
        snprintf (str,20,"%d",cbuf[AM.sbufnum].numrhs);
        PutPreVar(AM.oldnumextrasymbols,(UBYTE *)str,0,1);
        PutPreVar((UBYTE *)"optimvalue_",(UBYTE *)("0"),0,1);
        PutPreVar((UBYTE *)"optimscheme_",(UBYTE *)("0"),0,1);
        return(error);
}

/*
          #] ClearOptimize : 
*/

tueda / form / 15241916852

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous