"Fossies" - the Fresh Open Source Software Archive

Member "fityk-1.3.1/fityk/vm.cpp" (13 May 2016, 24233 Bytes) of package /linux/misc/fityk-1.3.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "vm.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.3.0_vs_1.3.1.

    1 // This file is part of fityk program. Copyright 2001-2013 Marcin Wojdyr
    2 // Licence: GNU General Public License ver. 2+
    3 /// virtual machine - calculates expressions using by executing bytecode
    4 
    5 #define BUILDING_LIBFITYK
    6 #include "vm.h"
    7 
    8 #include <boost/math/special_functions/gamma.hpp>
    9 #include <boost/math/special_functions/digamma.hpp>
   10 
   11 #include "common.h"
   12 #include "voigt.h"
   13 #include "numfuncs.h"
   14 #include "logic.h"
   15 
   16 #include "func.h" // %f(...)
   17 #include "model.h" // F(...)
   18 
   19 using namespace std;
   20 
   21 namespace {
   22 using namespace fityk;
   23 
   24 vector<int>::const_iterator
   25 skip_code(vector<int>::const_iterator i, int start_op, int finish_op)
   26 {
   27     //TODO: it doesn't work correctly now, because both op's and indices
   28     // are positive and mixed
   29     int counter = 1;
   30     while (counter) {
   31         ++i;
   32         if (*i == finish_op)
   33             counter--;
   34         else if (*i == start_op)
   35             counter++;
   36         else if (VMData::has_idx(*i))
   37             ++i;
   38     }
   39     return i;
   40 }
   41 
   42 template<typename T>
   43 realt get_var_with_idx(realt idx, vector<Point> const& points, T Point::*t)
   44 {
   45     if (points.empty())
   46         return 0.;
   47     else if (idx <= 0)
   48         return points[0].*t;
   49     else if (idx >= points.size() - 1)
   50         return points.back().*t;
   51     else if (is_eq(idx, iround(idx)))
   52         return points[iround(idx)].*t;
   53     else {
   54         int flo = int(floor(idx));
   55         realt fra = idx - flo;
   56         return (1-fra) * realt(points[flo].*t)
   57                + fra * realt(points[flo+1].*t);
   58     }
   59 }
   60 
   61 /// returns floating-point "index" of x in sorted vec of points
   62 realt find_idx_in_sorted(vector<Point> const& pp, realt x)
   63 {
   64     if (pp.empty())
   65         return 0;
   66     else if (x <= pp.front().x)
   67         return 0;
   68     else if (x >= pp.back().x)
   69         return pp.size() - 1;
   70     vector<Point>::const_iterator i = lower_bound(pp.begin(), pp.end(),
   71                                                   Point(x, 0));
   72     assert (i > pp.begin() && i < pp.end());
   73     if (is_eq(x, i->x))
   74         return i - pp.begin();
   75     else
   76         return i - pp.begin() - (i->x - x) / (i->x - (i-1)->x);
   77 }
   78 
   79 } // anonymous namespace
   80 
   81 namespace fityk {
   82 
   83 /// debuging utility
   84 #define OP_(x) \
   85     case OP_##x: return #x;
   86 
   87 string op2str(int op)
   88 {
   89     switch (static_cast<Op>(op)) {
   90         OP_(NUMBER) OP_(SYMBOL)
   91         OP_(X) OP_(PUT_DERIV)
   92         OP_(NEG)   OP_(EXP)  OP_(ERFC)  OP_(ERF)
   93         OP_(SIN)   OP_(COS)  OP_(TAN)  OP_(SINH) OP_(COSH)  OP_(TANH)
   94         OP_(ABS)  OP_(ROUND)
   95         OP_(ATAN) OP_(ASIN) OP_(ACOS)
   96         OP_(LOG10) OP_(LN)  OP_(SQRT)  OP_(POW)
   97         OP_(GAMMA) OP_(LGAMMA) OP_(DIGAMMA)
   98         OP_(VOIGT) OP_(DVOIGT_DX) OP_(DVOIGT_DY)
   99         OP_(XINDEX)
  100         OP_(ADD)   OP_(SUB)   OP_(MUL)   OP_(DIV)  OP_(MOD)
  101         OP_(MIN2)   OP_(MAX2) OP_(RANDNORM) OP_(RANDU)
  102         OP_(PX) OP_(PY) OP_(PS) OP_(PA)
  103         OP_(Px) OP_(Py) OP_(Ps) OP_(Pa)
  104         OP_(Pn) OP_(PM)
  105         OP_(OR) OP_(AFTER_OR) OP_(AND) OP_(AFTER_AND) OP_(NOT)
  106         OP_(TERNARY) OP_(TERNARY_MID) OP_(AFTER_TERNARY)
  107         OP_(GT) OP_(GE) OP_(LT) OP_(LE) OP_(EQ) OP_(NEQ)
  108         OP_(ASSIGN_X) OP_(ASSIGN_Y) OP_(ASSIGN_S) OP_(ASSIGN_A)
  109         OP_(FUNC) OP_(SUM_F) OP_(SUM_Z)
  110         OP_(NUMAREA) OP_(FINDX) OP_(FIND_EXTR)
  111         OP_(TILDE)
  112         OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG)
  113         OP_(OPEN_ROUND)  OP_(OPEN_SQUARE)
  114     }
  115     return S(op); // unreachable (if all OPs are listed above)
  116 }
  117 #undef OP_
  118 
  119 string vm2str(vector<int> const& code, vector<realt> const& data)
  120 {
  121     string s;
  122     for (vector<int>::const_iterator i = code.begin(); i < code.end(); ++i) {
  123         s += op2str(*i);
  124         if (*i == OP_NUMBER) {
  125             ++i;
  126             assert (*i >= 0 && *i < size(data));
  127             s += "[" + S(*i) + "](" + S(data[*i]) + ")";
  128         } else if (VMData::has_idx(*i)) {
  129             ++i;
  130             s += "[" + S(*i) + "]";
  131         }
  132         s += " ";
  133     }
  134     return s;
  135 }
  136 
  137 
  138 void VMData::append_number(realt d)
  139 {
  140     append_code(OP_NUMBER);
  141     int number_pos = numbers_.size();
  142     append_code(number_pos);
  143     numbers_.push_back(d);
  144 }
  145 
  146 void VMData::replace_symbols(const vector<realt>& vv)
  147 {
  148     for (vector<int>::iterator op = code_.begin(); op < code_.end(); ++op) {
  149         if (*op == OP_SYMBOL) {
  150             *op = OP_NUMBER;
  151             ++op;
  152             realt value = vv[*op];
  153             vector<realt>::iterator x =
  154                     find(numbers_.begin(), numbers_.end(), value);
  155             if (x != numbers_.end())
  156                 *op = x - numbers_.begin();
  157             else {
  158                 *op = numbers_.size();
  159                 numbers_.push_back(value);
  160             }
  161         } else if (has_idx(*op))
  162             ++op;
  163     }
  164 }
  165 
  166 /// switches between non-negative and negative indices (a -> -1-a),
  167 /// the point of having negative indices is to avoid conflicts with opcodes.
  168 /// The same transformation is used in OpTree. 
  169 void VMData::flip_indices()
  170 {
  171     for (vector<int>::iterator i = code_.begin(); i < code_.end(); ++i)
  172         if (has_idx(*i)) {
  173             ++i;
  174             *i = -1 - *i;
  175         }
  176 }
  177 
  178 bool VMData::has_op(int op) const
  179 {
  180     for (vector<int>::const_iterator i = code_.begin(); i < code_.end(); ++i) {
  181         if (*i == op)
  182             return true;
  183         if (has_idx(*i))
  184             ++i;
  185     }
  186     return false;
  187 }
  188 
  189 inline realt as_bool(realt d) { return fabs(d) < 0.5 ? 0 : 1; }
  190 
  191 #define STACK_OFFSET_CHANGE(ch) stackPtr+=(ch)
  192 
  193 inline
  194 void run_const_op(const Full* F, const std::vector<realt>& numbers,
  195                   vector<int>::const_iterator& i,
  196                   realt*& stackPtr,
  197                   const int n,
  198                   const vector<Point>& old_points,
  199                   const vector<Point>& new_points)
  200 {
  201     switch (*i) {
  202         //unary-operators
  203         case OP_NEG:
  204             *stackPtr = - *stackPtr;
  205             break;
  206         case OP_SQRT:
  207             *stackPtr = sqrt(*stackPtr);
  208             break;
  209         case OP_GAMMA:
  210             *stackPtr = boost::math::tgamma(*stackPtr);
  211             break;
  212         case OP_LGAMMA:
  213             *stackPtr = boost::math::lgamma(*stackPtr);
  214             break;
  215         case OP_DIGAMMA:
  216             *stackPtr = boost::math::digamma(*stackPtr);
  217             break;
  218         case OP_EXP:
  219             *stackPtr = exp(*stackPtr);
  220             break;
  221         case OP_ERFC:
  222             *stackPtr = erfc(*stackPtr);
  223             break;
  224         case OP_ERF:
  225             *stackPtr = erf(*stackPtr);
  226             break;
  227         case OP_LOG10:
  228             *stackPtr = log10(*stackPtr);
  229             break;
  230         case OP_LN:
  231             *stackPtr = log(*stackPtr);
  232             break;
  233         case OP_SIN:
  234             *stackPtr = sin(*stackPtr);
  235             break;
  236         case OP_COS:
  237             *stackPtr = cos(*stackPtr);
  238             break;
  239         case OP_TAN:
  240             *stackPtr = tan(*stackPtr);
  241             break;
  242         case OP_SINH:
  243             *stackPtr = sinh(*stackPtr);
  244             break;
  245         case OP_COSH:
  246             *stackPtr = cosh(*stackPtr);
  247             break;
  248         case OP_TANH:
  249             *stackPtr = tanh(*stackPtr);
  250             break;
  251         case OP_ATAN:
  252             *stackPtr = atan(*stackPtr);
  253             break;
  254         case OP_ASIN:
  255             *stackPtr = asin(*stackPtr);
  256             break;
  257         case OP_ACOS:
  258             *stackPtr = acos(*stackPtr);
  259             break;
  260         case OP_ABS:
  261             *stackPtr = fabs(*stackPtr);
  262             break;
  263         case OP_ROUND:
  264             *stackPtr = floor(*stackPtr + 0.5);
  265             break;
  266 
  267         case OP_XINDEX:
  268             *stackPtr = find_idx_in_sorted(old_points, *stackPtr);
  269             break;
  270 
  271 #ifndef STANDALONE_DATATRANS
  272         case OP_FUNC:
  273             i++;
  274             *stackPtr = F->mgr.get_function(*i)->calculate_value(*stackPtr);
  275             break;
  276         case OP_SUM_F:
  277             i++;
  278             *stackPtr = F->dk.get_model(*i)->value(*stackPtr);
  279             break;
  280         case OP_SUM_Z:
  281             i++;
  282             *stackPtr = F->dk.get_model(*i)->zero_shift(*stackPtr);
  283             break;
  284         case OP_NUMAREA:
  285             i += 2;
  286             STACK_OFFSET_CHANGE(-2);
  287             if (*(i-1) == OP_FUNC) {
  288                 *stackPtr = F->mgr.get_function(*i)->numarea(*stackPtr,
  289                                     *(stackPtr+1), iround(*(stackPtr+2)));
  290             } else if (*(i-1) == OP_SUM_F) {
  291                 *stackPtr = F->dk.get_model(*i)->numarea(*stackPtr,
  292                                     *(stackPtr+1), iround(*(stackPtr+2)));
  293             } else // OP_SUM_Z
  294                 throw ExecuteError("Z.numarea() is not implemented. "
  295                                    "Does anyone need it?");
  296             break;
  297 
  298         case OP_FINDX:
  299             i += 2;
  300             STACK_OFFSET_CHANGE(-2);
  301             if (*(i-1) == OP_FUNC) {
  302                 *stackPtr = find_x_with_value(F->mgr.get_function(*i),
  303                                   *stackPtr, *(stackPtr+1), *(stackPtr+2));
  304             } else if (*(i-1) == OP_SUM_F) {
  305                 *stackPtr = find_x_with_value(F->dk.get_model(*i),
  306                                   *stackPtr, *(stackPtr+1), *(stackPtr+2));
  307             } else // OP_SUM_Z
  308                 throw ExecuteError("Z.findx() is not implemented. "
  309                                    "Does anyone need it?");
  310             break;
  311 
  312         case OP_FIND_EXTR:
  313             i += 2;
  314             STACK_OFFSET_CHANGE(-1);
  315             if (*(i-1) == OP_FUNC) {
  316                 *stackPtr = find_extremum(F->mgr.get_function(*i),
  317                                           *stackPtr, *(stackPtr+1));
  318             } else if (*(i-1) == OP_SUM_F) {
  319                 *stackPtr = find_extremum(F->dk.get_model(*i),
  320                                           *stackPtr, *(stackPtr+1));
  321             } else // OP_SUM_Z
  322                 throw ExecuteError("Z.extremum() is not implemented. "
  323                                    "Does anyone need it?");
  324             break;
  325 #endif //not STANDALONE_DATATRANS
  326 
  327         //binary-operators
  328         case OP_MIN2:
  329             STACK_OFFSET_CHANGE(-1);
  330             *stackPtr = min(*stackPtr, *(stackPtr+1));
  331             break;
  332         case OP_MAX2:
  333             STACK_OFFSET_CHANGE(-1);
  334             *stackPtr = max(*stackPtr, *(stackPtr+1));
  335             break;
  336         case OP_RANDU:
  337             STACK_OFFSET_CHANGE(-1);
  338             *stackPtr = rand_uniform(*stackPtr, *(stackPtr+1));
  339             break;
  340         case OP_RANDNORM:
  341             STACK_OFFSET_CHANGE(-1);
  342             *stackPtr += rand_gauss() * *(stackPtr+1);
  343             break;
  344         case OP_ADD:
  345             STACK_OFFSET_CHANGE(-1);
  346             *stackPtr += *(stackPtr+1);
  347             break;
  348         case OP_SUB:
  349             STACK_OFFSET_CHANGE(-1);
  350             *stackPtr -= *(stackPtr+1);
  351             break;
  352         case OP_MUL:
  353             STACK_OFFSET_CHANGE(-1);
  354             *stackPtr *= *(stackPtr+1);
  355             break;
  356         case OP_DIV:
  357             STACK_OFFSET_CHANGE(-1);
  358             *stackPtr /= *(stackPtr+1);
  359             break;
  360         case OP_MOD:
  361             STACK_OFFSET_CHANGE(-1);
  362             *stackPtr -= floor(*stackPtr / *(stackPtr+1)) * *(stackPtr+1);
  363             break;
  364         case OP_POW:
  365             STACK_OFFSET_CHANGE(-1);
  366             *stackPtr = pow(*stackPtr, *(stackPtr+1));
  367             break;
  368         case OP_VOIGT:
  369             STACK_OFFSET_CHANGE(-1);
  370             *stackPtr = humlik(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  371             break;
  372         case OP_DVOIGT_DX:
  373             STACK_OFFSET_CHANGE(-1);
  374             *stackPtr = humdev_dkdx(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  375             break;
  376         case OP_DVOIGT_DY:
  377             STACK_OFFSET_CHANGE(-1);
  378             *stackPtr = humdev_dkdy(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  379             break;
  380 
  381         // comparisions
  382         case OP_LT:
  383             STACK_OFFSET_CHANGE(-1);
  384             *stackPtr = is_lt(*stackPtr, *(stackPtr+1));
  385             break;
  386         case OP_GT:
  387             STACK_OFFSET_CHANGE(-1);
  388             *stackPtr = is_gt(*stackPtr, *(stackPtr+1));
  389             break;
  390         case OP_LE:
  391             STACK_OFFSET_CHANGE(-1);
  392             *stackPtr = is_le(*stackPtr, *(stackPtr+1));
  393             break;
  394         case OP_GE:
  395             STACK_OFFSET_CHANGE(-1);
  396             *stackPtr = is_ge(*stackPtr, *(stackPtr+1));
  397             break;
  398         case OP_EQ:
  399             STACK_OFFSET_CHANGE(-1);
  400             *stackPtr = is_eq(*stackPtr, *(stackPtr+1));
  401             break;
  402         case OP_NEQ:
  403             STACK_OFFSET_CHANGE(-1);
  404             *stackPtr = is_neq(*stackPtr, *(stackPtr+1));
  405             break;
  406 
  407         // putting-number-to-stack-operators
  408         case OP_NUMBER:
  409             STACK_OFFSET_CHANGE(+1);
  410             i++;
  411             *stackPtr = numbers[*i];
  412             break;
  413 
  414         case OP_Pn:
  415             STACK_OFFSET_CHANGE(+1);
  416             *stackPtr = static_cast<realt>(n);
  417             break;
  418         case OP_PM:
  419             STACK_OFFSET_CHANGE(+1);
  420             *stackPtr = static_cast<realt>(old_points.size());
  421             break;
  422 
  423         case OP_Px:
  424             *stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::x);
  425             break;
  426         case OP_Py:
  427             *stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::y);
  428             break;
  429         case OP_Ps:
  430             *stackPtr = get_var_with_idx(*stackPtr, old_points,
  431                                          &Point::sigma);
  432             break;
  433         case OP_Pa:
  434             *stackPtr = as_bool(get_var_with_idx(*stackPtr, old_points,
  435                                                  &Point::is_active));
  436             break;
  437         case OP_PX:
  438             *stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::x);
  439             break;
  440         case OP_PY:
  441             *stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::y);
  442             break;
  443         case OP_PS:
  444             *stackPtr = get_var_with_idx(*stackPtr, new_points,
  445                                          &Point::sigma);
  446             break;
  447         case OP_PA:
  448             *stackPtr = as_bool(get_var_with_idx(*stackPtr, new_points,
  449                                                  &Point::is_active));
  450             break;
  451 
  452         // logical; can skip part of VM code !
  453         case OP_NOT:
  454             *stackPtr = is_eq(*stackPtr, 0.);
  455             break;
  456         case OP_AND:
  457             if (is_neq(*stackPtr, 0))    //return second
  458                 stackPtr--;
  459             else              // return first
  460                 i = skip_code(i, OP_AND, OP_AFTER_AND);
  461             break;
  462 
  463         case OP_OR:
  464             if (is_neq(*stackPtr, 0))    //return first
  465                 i = skip_code(i, OP_OR, OP_AFTER_OR);
  466             else              // return second
  467                 stackPtr--;
  468             break;
  469 
  470         case OP_TERNARY:
  471             if (! *stackPtr)
  472                 i = skip_code(i, OP_TERNARY, OP_TERNARY_MID);
  473             stackPtr--;
  474             break;
  475         case OP_TERNARY_MID:
  476             //if we are here, condition was true. Skip.
  477             i = skip_code(i, OP_TERNARY_MID, OP_AFTER_TERNARY);
  478             break;
  479 
  480         case OP_AFTER_AND: //do nothing
  481         case OP_AFTER_OR:
  482         case OP_AFTER_TERNARY:
  483         case OP_TILDE:
  484             break;
  485 
  486         case OP_DATASET:
  487             throw ExecuteError("@n can not be used in this context");
  488             //break; // unreachable
  489 
  490         default:
  491             //cerr << "Unknown operator in VM code: " << *i << endl;
  492             assert(0);
  493     }
  494 }
  495 
  496 inline
  497 void run_mutab_op(const Full* F, const std::vector<realt>& numbers,
  498                   vector<int>::const_iterator& i,
  499                   realt*& stackPtr,
  500                   const int n,
  501                   const vector<Point>& old_points,
  502                   vector<Point>& new_points)
  503 {
  504     //cerr << "op " << dt_op(*i) << endl;
  505     switch (*i) {
  506         //assignment-operators
  507         case OP_ASSIGN_X:
  508             new_points[n].x = *stackPtr;
  509             STACK_OFFSET_CHANGE(-1);
  510             break;
  511         case OP_ASSIGN_Y:
  512             new_points[n].y = *stackPtr;
  513             STACK_OFFSET_CHANGE(-1);
  514             break;
  515         case OP_ASSIGN_S:
  516             new_points[n].sigma = *stackPtr;
  517             STACK_OFFSET_CHANGE(-1);
  518             break;
  519         case OP_ASSIGN_A:
  520             new_points[n].is_active = is_neq(*stackPtr, 0.);
  521             STACK_OFFSET_CHANGE(-1);
  522             break;
  523         default:
  524             run_const_op(F, numbers, i, stackPtr, n, old_points, new_points);
  525     }
  526 }
  527 
  528 inline
  529 void run_func_op(const vector<realt>& numbers, vector<int>::const_iterator &i,
  530                  realt*& stackPtr)
  531 {
  532     switch (*i) {
  533         //unary operators
  534         case OP_NEG:
  535             *stackPtr = - *stackPtr;
  536             break;
  537         case OP_SQRT:
  538             *stackPtr = sqrt(*stackPtr);
  539             break;
  540         case OP_EXP:
  541             *stackPtr = exp(*stackPtr);
  542             break;
  543         case OP_ERFC:
  544             *stackPtr = erfc(*stackPtr);
  545             break;
  546         case OP_ERF:
  547             *stackPtr = erf(*stackPtr);
  548             break;
  549         case OP_LOG10:
  550             *stackPtr = log10(*stackPtr);
  551             break;
  552         case OP_LN:
  553             *stackPtr = log(*stackPtr);
  554             break;
  555         case OP_SINH:
  556             *stackPtr = sinh(*stackPtr);
  557             break;
  558         case OP_COSH:
  559             *stackPtr = cosh(*stackPtr);
  560             break;
  561         case OP_TANH:
  562             *stackPtr = tanh(*stackPtr);
  563             break;
  564         case OP_SIN:
  565             *stackPtr = sin(*stackPtr);
  566             break;
  567         case OP_COS:
  568             *stackPtr = cos(*stackPtr);
  569             break;
  570         case OP_TAN:
  571             *stackPtr = tan(*stackPtr);
  572             break;
  573         case OP_ATAN:
  574             *stackPtr = atan(*stackPtr);
  575             break;
  576         case OP_ASIN:
  577             *stackPtr = asin(*stackPtr);
  578             break;
  579         case OP_ACOS:
  580             *stackPtr = acos(*stackPtr);
  581             break;
  582         case OP_LGAMMA:
  583             *stackPtr = boost::math::lgamma(*stackPtr);
  584             break;
  585         case OP_DIGAMMA:
  586             *stackPtr = boost::math::digamma(*stackPtr);
  587             break;
  588         case OP_ABS:
  589             *stackPtr = fabs(*stackPtr);
  590             break;
  591 
  592         //binary operators
  593         case OP_ADD:
  594             STACK_OFFSET_CHANGE(-1);
  595             *stackPtr += *(stackPtr+1);
  596             break;
  597         case OP_SUB:
  598             STACK_OFFSET_CHANGE(-1);
  599             *stackPtr -= *(stackPtr+1);
  600             break;
  601         case OP_MUL:
  602             STACK_OFFSET_CHANGE(-1);
  603             *stackPtr *= *(stackPtr+1);
  604             break;
  605         case OP_DIV:
  606             STACK_OFFSET_CHANGE(-1);
  607             *stackPtr /= *(stackPtr+1);
  608             break;
  609         case OP_POW:
  610             STACK_OFFSET_CHANGE(-1);
  611             *stackPtr = pow(*stackPtr, *(stackPtr+1));
  612             break;
  613         case OP_VOIGT:
  614             STACK_OFFSET_CHANGE(-1);
  615             *stackPtr = humlik(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  616             break;
  617         case OP_DVOIGT_DX:
  618             STACK_OFFSET_CHANGE(-1);
  619             *stackPtr = humdev_dkdx(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  620             break;
  621         case OP_DVOIGT_DY:
  622             STACK_OFFSET_CHANGE(-1);
  623             *stackPtr = humdev_dkdy(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
  624             break;
  625 
  626         // putting-number-to-stack-operators
  627         // stack overflow not checked
  628         case OP_NUMBER:
  629             STACK_OFFSET_CHANGE(+1);
  630             i++; // OP_NUMBER opcode is always followed by index
  631             *stackPtr = numbers[*i];
  632             break;
  633 
  634         case OP_TILDE:
  635             // used for default values in Runner::make_func_from_template()
  636             break;
  637 
  638         default:
  639             throw ExecuteError("op " + op2str(*i) +
  640                                " is not allowed for variables and functions");
  641     }
  642 }
  643 
  644 
  645 void ExprCalculator::transform_data(vector<Point>& points)
  646 {
  647     if (points.empty())
  648         return;
  649 
  650     realt stack[16];
  651     realt* stackPtr = stack - 1; // will be ++'ed first
  652     vector<Point> new_points = points;
  653 
  654     // do the time-consuming overflow checking only for the first point
  655     v_foreach (int, i, vm_.code()) {
  656         run_mutab_op(F_, vm_.numbers(), i, stackPtr, 0, points, new_points);
  657         if (stackPtr - stack >= 16)
  658             throw ExecuteError("stack overflow");
  659     }
  660     assert(stackPtr == stack - 1); // ASSIGN_ op must be at the end
  661 
  662     // the same for the rest of points, but without checks
  663     for (int n = 1; n != size(points); ++n)
  664         v_foreach (int, i, vm_.code())
  665             run_mutab_op(F_, vm_.numbers(), i, stackPtr, n, points, new_points);
  666 
  667     points.swap(new_points);
  668 }
  669 
  670 realt ExprCalculator::calculate(int n, const vector<Point>& points) const
  671 {
  672     realt stack[16];
  673     realt* stackPtr = stack - 1; // will be ++'ed first
  674     v_foreach (int, i, vm_.code()) {
  675         run_const_op(F_, vm_.numbers(), i, stackPtr, n, points, points);
  676         if (stackPtr - stack >= 16)
  677             throw ExecuteError("stack overflow");
  678     }
  679     //cerr << "stackPtr: " << stackPtr - stack << endl;
  680     assert(stackPtr == stack); // no ASSIGN_ at the end
  681     return stack[0];
  682 }
  683 
  684 realt ExprCalculator::calculate_custom(const vector<realt>& custom_val) const
  685 {
  686     realt stack[16];
  687     realt* stackPtr = stack - 1; // will be ++'ed first
  688     const vector<Point> dummy;
  689     v_foreach (int, i, vm_.code()) {
  690         if (*i == OP_SYMBOL) {
  691             STACK_OFFSET_CHANGE(+1);
  692             ++i;
  693             if (is_index(*i, custom_val))
  694                 *stackPtr = custom_val[*i];
  695             else
  696                 throw ExecuteError("[internal] variable mismatch");
  697         } else
  698             run_const_op(F_, vm_.numbers(), i, stackPtr, 0, dummy, dummy);
  699         if (stackPtr - stack >= 16)
  700             throw ExecuteError("stack overflow");
  701     }
  702     assert(stackPtr == stack); // no ASSIGN_ at the end
  703     return stack[0];
  704 }
  705 
  706 /// executes VM code, sets derivatives and returns value
  707 realt run_code_for_variable(const VMData& vm,
  708                             const vector<Variable*> &variables,
  709                             vector<realt> &derivatives)
  710 {
  711     realt stack[16];
  712     realt* stackPtr = stack - 1; // will be ++'ed first
  713     v_foreach (int, i, vm.code()) {
  714         if (*i == OP_SYMBOL) {
  715             STACK_OFFSET_CHANGE(+1);
  716             ++i; // skip the next one
  717             *stackPtr = variables[*i]->value();
  718         } else if (*i == OP_PUT_DERIV) {
  719             ++i;
  720             // the OP_PUT_DERIV opcode is followed by a number n,
  721             // the derivative is calculated with respect to n'th variable
  722             assert(*i < (int) derivatives.size());
  723             derivatives[*i] = *stackPtr;
  724             STACK_OFFSET_CHANGE(-1);
  725         } else
  726             run_func_op(vm.numbers(), i, stackPtr);
  727     }
  728     assert(stackPtr == stack);
  729     return stack[0];
  730 }
  731 
  732 
  733 realt run_code_for_custom_func(const VMData& vm, realt x,
  734                                vector<realt> &derivatives)
  735 {
  736     realt stack[16];
  737     realt* stackPtr = stack - 1; // will be ++'ed first
  738     v_foreach (int, i, vm.code()) {
  739         if (*i == OP_X) {
  740             STACK_OFFSET_CHANGE(+1);
  741             *stackPtr = x;
  742         } else if (*i == OP_PUT_DERIV) {
  743             ++i;
  744             // the OP_PUT_DERIV opcode is followed by a number n,
  745             // the derivative is calculated with respect to n'th variable
  746             assert(*i < (int) derivatives.size());
  747             derivatives[*i] = *stackPtr;
  748             STACK_OFFSET_CHANGE(-1);
  749         } else
  750             run_func_op(vm.numbers(), i, stackPtr);
  751     }
  752     assert(stackPtr == stack);
  753     return stack[0];
  754 }
  755 
  756 realt run_code_for_custom_func_value(const VMData& vm, realt x,
  757                                      int code_offset)
  758 {
  759     realt stack[16];
  760     realt* stackPtr = stack - 1; // will be ++'ed first
  761     for (vector<int>::const_iterator i = vm.code().begin() + code_offset;
  762                                                  i != vm.code().end(); ++i) {
  763         if (*i == OP_X) {
  764             STACK_OFFSET_CHANGE(+1);
  765             *stackPtr = x;
  766         } else
  767             run_func_op(vm.numbers(), i, stackPtr);
  768     }
  769     assert(stackPtr == stack);
  770     return stack[0];
  771 }
  772 
  773 } // namespace fityk