fityk: /home/wojdyr/fityk/src/vm.cpp Source File

00001 // This file is part of fityk program. Copyright (C) Marcin Wojdyr
00002 // Licence: GNU General Public License ver. 2+
00003 /// virtual machine - calculates expressions using by executing bytecode
00004 
00005 #include "vm.h"
00006 
00007 #include <boost/math/special_functions/gamma.hpp>
00008 #include <boost/math/special_functions/digamma.hpp>
00009 
00010 #include "common.h"
00011 #include "voigt.h"
00012 #include "numfuncs.h"
00013 #include "logic.h"
00014 
00015 #include "func.h" // %f(...)
00016 #include "model.h" // F(...)
00017 
00018 using namespace std;
00019 
00020 namespace {
00021 
00022 vector<int>::const_iterator
00023 skip_code(vector<int>::const_iterator i, int start_op, int finish_op)
00024 {
00025     //TODO: it doesn't work correctly now, because both op's and indices
00026     // are positive and mixed
00027     int counter = 1;
00028     while (counter) {
00029         ++i;
00030         if (*i == finish_op)
00031             counter--;
00032         else if (*i == start_op)
00033             counter++;
00034         else if (VMData::has_idx(*i))
00035             ++i;
00036     }
00037     return i;
00038 }
00039 
00040 template<typename T>
00041 realt get_var_with_idx(realt idx, vector<Point> const& points, T Point::*t)
00042 {
00043     if (points.empty())
00044         return 0.;
00045     else if (idx <= 0)
00046         return points[0].*t;
00047     else if (idx >= points.size() - 1)
00048         return points.back().*t;
00049     else if (is_eq(idx, iround(idx)))
00050         return points[iround(idx)].*t;
00051     else {
00052         int flo = int(floor(idx));
00053         realt fra = idx - flo;
00054         return (1-fra) * realt(points[flo].*t)
00055                + fra * realt(points[flo+1].*t);
00056     }
00057 }
00058 
00059 /// returns floating-point "index" of x in sorted vec of points
00060 realt find_idx_in_sorted(vector<Point> const& pp, realt x)
00061 {
00062     if (pp.empty())
00063         return 0;
00064     else if (x <= pp.front().x)
00065         return 0;
00066     else if (x >= pp.back().x)
00067         return pp.size() - 1;
00068     vector<Point>::const_iterator i = lower_bound(pp.begin(), pp.end(),
00069                                                   Point(x, 0));
00070     assert (i > pp.begin() && i < pp.end());
00071     if (is_eq(x, i->x))
00072         return i - pp.begin();
00073     else
00074         return i - pp.begin() - (i->x - x) / (i->x - (i-1)->x);
00075 }
00076 
00077 } // anonymous namespace
00078 
00079 
00080 /// debuging utility
00081 #define OP_(x) \
00082     case OP_##x: return #x;
00083 
00084 string op2str(int op)
00085 {
00086     switch (static_cast<Op>(op)) {
00087         OP_(NUMBER) OP_(SYMBOL)
00088         OP_(X) OP_(PUT_DERIV)
00089         OP_(NEG)   OP_(EXP)  OP_(ERFC)  OP_(ERF)
00090         OP_(SIN)   OP_(COS)  OP_(TAN)  OP_(SINH) OP_(COSH)  OP_(TANH)
00091         OP_(ABS)  OP_(ROUND)
00092         OP_(ATAN) OP_(ASIN) OP_(ACOS)
00093         OP_(LOG10) OP_(LN)  OP_(SQRT)  OP_(POW)
00094         OP_(GAMMA) OP_(LGAMMA) OP_(DIGAMMA)
00095         OP_(VOIGT) OP_(DVOIGT_DX) OP_(DVOIGT_DY)
00096         OP_(XINDEX)
00097         OP_(ADD)   OP_(SUB)   OP_(MUL)   OP_(DIV)  OP_(MOD)
00098         OP_(MIN2)   OP_(MAX2) OP_(RANDNORM) OP_(RANDU)
00099         OP_(PX) OP_(PY) OP_(PS) OP_(PA)
00100         OP_(Px) OP_(Py) OP_(Ps) OP_(Pa)
00101         OP_(Pn) OP_(PM)
00102         OP_(OR) OP_(AFTER_OR) OP_(AND) OP_(AFTER_AND) OP_(NOT)
00103         OP_(TERNARY) OP_(TERNARY_MID) OP_(AFTER_TERNARY)
00104         OP_(GT) OP_(GE) OP_(LT) OP_(LE) OP_(EQ) OP_(NEQ)
00105         OP_(ASSIGN_X) OP_(ASSIGN_Y) OP_(ASSIGN_S) OP_(ASSIGN_A)
00106         OP_(FUNC) OP_(SUM_F) OP_(SUM_Z)
00107         OP_(NUMAREA) OP_(FINDX) OP_(FIND_EXTR)
00108         OP_(TILDE)
00109         OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG)
00110         OP_(OPEN_ROUND)  OP_(OPEN_SQUARE)
00111     }
00112     return S(op);
00113 };
00114 #undef OP_
00115 
00116 string vm2str(vector<int> const& code, vector<realt> const& data)
00117 {
00118     string s;
00119     v_foreach (int, i, code) {
00120         s += op2str(*i);
00121         if (*i == OP_NUMBER) {
00122             ++i;
00123             assert (*i >= 0 && *i < size(data));
00124             s += "[" + S(*i) + "](" + S(data[*i]) + ")";
00125         }
00126         else if (VMData::has_idx(*i)) {
00127             ++i;
00128             s += "[" + S(*i) + "]";
00129         }
00130         s += " ";
00131     }
00132     return s;
00133 }
00134 
00135 
00136 void VMData::append_number(realt d)
00137 {
00138     append_code(OP_NUMBER);
00139     int number_pos = numbers_.size();
00140     append_code(number_pos);
00141     numbers_.push_back(d);
00142 }
00143 
00144 void VMData::replace_symbols(const vector<realt>& vv)
00145 {
00146     vm_foreach (int, op, code_) {
00147         if (*op == OP_SYMBOL) {
00148             *op = OP_NUMBER;
00149             ++op;
00150             realt value = vv[*op];
00151             vector<realt>::iterator x =
00152                     find(numbers_.begin(), numbers_.end(), value);
00153             if (x != numbers_.end())
00154                 *op = x - numbers_.begin();
00155             else {
00156                 *op = numbers_.size();
00157                 numbers_.push_back(value);
00158             }
00159         }
00160         else if (has_idx(*op))
00161             ++op;
00162     }
00163 }
00164 
00165 /// switches between non-negative and negative indices (a -> -1-a),
00166 /// the point of having negative indices is to avoid conflicts with opcodes.
00167 /// The same transformation is used in OpTree. 
00168 void VMData::flip_indices()
00169 {
00170     vm_foreach (int, i, code_)
00171         if (has_idx(*i)) {
00172             ++i;
00173             *i = -1 - *i;
00174         }
00175 }
00176 
00177 bool VMData::has_op(int op) const
00178 {
00179     v_foreach (int, i, code_) {
00180         if (*i == op)
00181             return true;
00182         if (has_idx(*i))
00183             ++i;
00184     }
00185     return false;
00186 }
00187 
00188 #define STACK_OFFSET_CHANGE(ch) stackPtr+=(ch)
00189 
00190 inline
00191 void run_const_op(const Ftk* F, const std::vector<realt>& numbers,
00192                   vector<int>::const_iterator& i,
00193                   realt*& stackPtr,
00194                   const int n,
00195                   const vector<Point>& old_points,
00196                   const vector<Point>& new_points)
00197 {
00198     switch (*i) {
00199         //unary-operators
00200         case OP_NEG:
00201             *stackPtr = - *stackPtr;
00202             break;
00203         case OP_SQRT:
00204             *stackPtr = sqrt(*stackPtr);
00205             break;
00206         case OP_GAMMA:
00207             *stackPtr = boost::math::tgamma(*stackPtr);
00208             break;
00209         case OP_LGAMMA:
00210             *stackPtr = boost::math::lgamma(*stackPtr);
00211             break;
00212         case OP_DIGAMMA:
00213             *stackPtr = boost::math::digamma(*stackPtr);
00214             break;
00215         case OP_EXP:
00216             *stackPtr = exp(*stackPtr);
00217             break;
00218         case OP_ERFC:
00219             *stackPtr = erfc(*stackPtr);
00220             break;
00221         case OP_ERF:
00222             *stackPtr = erf(*stackPtr);
00223             break;
00224         case OP_LOG10:
00225             *stackPtr = log10(*stackPtr);
00226             break;
00227         case OP_LN:
00228             *stackPtr = log(*stackPtr);
00229             break;
00230         case OP_SIN:
00231             *stackPtr = sin(*stackPtr);
00232             break;
00233         case OP_COS:
00234             *stackPtr = cos(*stackPtr);
00235             break;
00236         case OP_TAN:
00237             *stackPtr = tan(*stackPtr);
00238             break;
00239         case OP_SINH:
00240             *stackPtr = sinh(*stackPtr);
00241             break;
00242         case OP_COSH:
00243             *stackPtr = cosh(*stackPtr);
00244             break;
00245         case OP_TANH:
00246             *stackPtr = tanh(*stackPtr);
00247             break;
00248         case OP_ATAN:
00249             *stackPtr = atan(*stackPtr);
00250             break;
00251         case OP_ASIN:
00252             *stackPtr = asin(*stackPtr);
00253             break;
00254         case OP_ACOS:
00255             *stackPtr = acos(*stackPtr);
00256             break;
00257         case OP_ABS:
00258             *stackPtr = fabs(*stackPtr);
00259             break;
00260         case OP_ROUND:
00261             *stackPtr = floor(*stackPtr + 0.5);
00262             break;
00263 
00264         case OP_XINDEX:
00265             *stackPtr = find_idx_in_sorted(old_points, *stackPtr);
00266             break;
00267 
00268 #ifndef STANDALONE_DATATRANS
00269         case OP_FUNC:
00270             i++;
00271             *stackPtr = F->get_function(*i)->calculate_value(*stackPtr);
00272             break;
00273         case OP_SUM_F:
00274             i++;
00275             *stackPtr = F->get_model(*i)->value(*stackPtr);
00276             break;
00277         case OP_SUM_Z:
00278             i++;
00279             *stackPtr = F->get_model(*i)->zero_shift(*stackPtr);
00280             break;
00281         case OP_NUMAREA:
00282             i += 2;
00283             STACK_OFFSET_CHANGE(-2);
00284             if (*(i-1) == OP_FUNC) {
00285                 *stackPtr = F->get_function(*i)->numarea(*stackPtr,
00286                                     *(stackPtr+1), iround(*(stackPtr+2)));
00287             }
00288             else if (*(i-1) == OP_SUM_F) {
00289                 *stackPtr = F->get_model(*i)->numarea(*stackPtr,
00290                                     *(stackPtr+1), iround(*(stackPtr+2)));
00291             }
00292             else // OP_SUM_Z
00293                 throw ExecuteError("Z.numarea() is not implemented. "
00294                                    "Does anyone need it?");
00295             break;
00296 
00297         case OP_FINDX:
00298             i += 2;
00299             STACK_OFFSET_CHANGE(-2);
00300             if (*(i-1) == OP_FUNC) {
00301                 *stackPtr = F->get_function(*i)->find_x_with_value(
00302                                   *stackPtr, *(stackPtr+1), *(stackPtr+2));
00303             }
00304             else if (*(i-1) == OP_SUM_F) {
00305                 throw ExecuteError("F.findx() is not implemented. "
00306                                    "Does anyone need it?");
00307             }
00308             else // OP_SUM_Z
00309                 throw ExecuteError("Z.findx() is not implemented. "
00310                                    "Does anyone need it?");
00311             break;
00312 
00313         case OP_FIND_EXTR:
00314             i += 2;
00315             STACK_OFFSET_CHANGE(-1);
00316             if (*(i-1) == OP_FUNC) {
00317                 *stackPtr = F->get_function(*i)->find_extremum(*stackPtr,
00318                                                             *(stackPtr+1));
00319             }
00320             else if (*(i-1) == OP_SUM_F) {
00321                 throw ExecuteError("F.extremum() is not implemented. "
00322                                    "Does anyone need it?");
00323             }
00324             else // OP_SUM_Z
00325                 throw ExecuteError("Z.extremum() is not implemented. "
00326                                    "Does anyone need it?");
00327             break;
00328 #endif //not STANDALONE_DATATRANS
00329 
00330         //binary-operators
00331         case OP_MIN2:
00332             STACK_OFFSET_CHANGE(-1);
00333             *stackPtr = min(*stackPtr, *(stackPtr+1));
00334             break;
00335         case OP_MAX2:
00336             STACK_OFFSET_CHANGE(-1);
00337             *stackPtr = max(*stackPtr, *(stackPtr+1));
00338             break;
00339         case OP_RANDU:
00340             STACK_OFFSET_CHANGE(-1);
00341             *stackPtr = rand_uniform(*stackPtr, *(stackPtr+1));
00342             break;
00343         case OP_RANDNORM:
00344             STACK_OFFSET_CHANGE(-1);
00345             *stackPtr += rand_gauss() * *(stackPtr+1);
00346             break;
00347         case OP_ADD:
00348             STACK_OFFSET_CHANGE(-1);
00349             *stackPtr += *(stackPtr+1);
00350             break;
00351         case OP_SUB:
00352             STACK_OFFSET_CHANGE(-1);
00353             *stackPtr -= *(stackPtr+1);
00354             break;
00355         case OP_MUL:
00356             STACK_OFFSET_CHANGE(-1);
00357             *stackPtr *= *(stackPtr+1);
00358             break;
00359         case OP_DIV:
00360             STACK_OFFSET_CHANGE(-1);
00361             *stackPtr /= *(stackPtr+1);
00362             break;
00363         case OP_MOD:
00364             STACK_OFFSET_CHANGE(-1);
00365             *stackPtr -= floor(*stackPtr / *(stackPtr+1)) * *(stackPtr+1);
00366             break;
00367         case OP_POW:
00368             STACK_OFFSET_CHANGE(-1);
00369             *stackPtr = pow(*stackPtr, *(stackPtr+1));
00370             break;
00371         case OP_VOIGT:
00372             STACK_OFFSET_CHANGE(-1);
00373             *stackPtr = humlik(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00374             break;
00375         case OP_DVOIGT_DX:
00376             STACK_OFFSET_CHANGE(-1);
00377             *stackPtr = humdev_dkdx(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00378             break;
00379         case OP_DVOIGT_DY:
00380             STACK_OFFSET_CHANGE(-1);
00381             *stackPtr = humdev_dkdy(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00382             break;
00383 
00384         // comparisions
00385         case OP_LT:
00386             STACK_OFFSET_CHANGE(-1);
00387             *stackPtr = is_lt(*stackPtr, *(stackPtr+1));
00388             break;
00389         case OP_GT:
00390             STACK_OFFSET_CHANGE(-1);
00391             *stackPtr = is_gt(*stackPtr, *(stackPtr+1));
00392             break;
00393         case OP_LE:
00394             STACK_OFFSET_CHANGE(-1);
00395             *stackPtr = is_le(*stackPtr, *(stackPtr+1));
00396             break;
00397         case OP_GE:
00398             STACK_OFFSET_CHANGE(-1);
00399             *stackPtr = is_ge(*stackPtr, *(stackPtr+1));
00400             break;
00401         case OP_EQ:
00402             STACK_OFFSET_CHANGE(-1);
00403             *stackPtr = is_eq(*stackPtr, *(stackPtr+1));
00404             break;
00405         case OP_NEQ:
00406             STACK_OFFSET_CHANGE(-1);
00407             *stackPtr = is_neq(*stackPtr, *(stackPtr+1));
00408             break;
00409 
00410         // putting-number-to-stack-operators
00411         case OP_NUMBER:
00412             STACK_OFFSET_CHANGE(+1);
00413             i++;
00414             *stackPtr = numbers[*i];
00415             break;
00416 
00417         case OP_Pn:
00418             STACK_OFFSET_CHANGE(+1);
00419             *stackPtr = static_cast<realt>(n);
00420             break;
00421         case OP_PM:
00422             STACK_OFFSET_CHANGE(+1);
00423             *stackPtr = static_cast<realt>(old_points.size());
00424             break;
00425 
00426         case OP_Px:
00427             *stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::x);
00428             break;
00429         case OP_Py:
00430             *stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::y);
00431             break;
00432         case OP_Ps:
00433             *stackPtr = get_var_with_idx(*stackPtr, old_points,
00434                                          &Point::sigma);
00435             break;
00436         case OP_Pa:
00437             *stackPtr = bool(iround(get_var_with_idx(*stackPtr, old_points,
00438                                                      &Point::is_active)));
00439             break;
00440         case OP_PX:
00441             *stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::x);
00442             break;
00443         case OP_PY:
00444             *stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::y);
00445             break;
00446         case OP_PS:
00447             *stackPtr = get_var_with_idx(*stackPtr, new_points,
00448                                          &Point::sigma);
00449             break;
00450         case OP_PA:
00451             *stackPtr = bool(iround(get_var_with_idx(*stackPtr, new_points,
00452                                                      &Point::is_active)));
00453             break;
00454 
00455         // logical; can skip part of VM code !
00456         case OP_NOT:
00457             *stackPtr = is_eq(*stackPtr, 0.);
00458             break;
00459         case OP_AND:
00460             if (is_neq(*stackPtr, 0))    //return second
00461                 stackPtr--;
00462             else              // return first
00463                 i = skip_code(i, OP_AND, OP_AFTER_AND);
00464             break;
00465 
00466         case OP_OR:
00467             if (is_neq(*stackPtr, 0))    //return first
00468                 i = skip_code(i, OP_OR, OP_AFTER_OR);
00469             else              // return second
00470                 stackPtr--;
00471             break;
00472 
00473         case OP_TERNARY:
00474             if (! *stackPtr)
00475                 i = skip_code(i, OP_TERNARY, OP_TERNARY_MID);
00476             stackPtr--;
00477             break;
00478         case OP_TERNARY_MID:
00479             //if we are here, condition was true. Skip.
00480             i = skip_code(i, OP_TERNARY_MID, OP_AFTER_TERNARY);
00481             break;
00482 
00483         case OP_AFTER_AND: //do nothing
00484         case OP_AFTER_OR:
00485         case OP_AFTER_TERNARY:
00486         case OP_TILDE:
00487             break;
00488 
00489         case OP_DATASET:
00490             throw ExecuteError("@n can not be used in this context");
00491             break;
00492 
00493         default:
00494             //cerr << "Unknown operator in VM code: " << *i << endl;
00495             assert(0);
00496     }
00497 }
00498 
00499 inline
00500 void run_mutab_op(const Ftk* F, const std::vector<realt>& numbers,
00501                   vector<int>::const_iterator& i,
00502                   realt*& stackPtr,
00503                   const int n,
00504                   const vector<Point>& old_points,
00505                   vector<Point>& new_points)
00506 {
00507     //cerr << "op " << dt_op(*i) << endl;
00508     switch (*i) {
00509         //assignment-operators
00510         case OP_ASSIGN_X:
00511             new_points[n].x = *stackPtr;
00512             STACK_OFFSET_CHANGE(-1);
00513             break;
00514         case OP_ASSIGN_Y:
00515             new_points[n].y = *stackPtr;
00516             STACK_OFFSET_CHANGE(-1);
00517             break;
00518         case OP_ASSIGN_S:
00519             new_points[n].sigma = *stackPtr;
00520             STACK_OFFSET_CHANGE(-1);
00521             break;
00522         case OP_ASSIGN_A:
00523             new_points[n].is_active = is_neq(*stackPtr, 0.);
00524             STACK_OFFSET_CHANGE(-1);
00525             break;
00526         default:
00527             run_const_op(F, numbers, i, stackPtr, n, old_points, new_points);
00528     }
00529 }
00530 
00531 inline
00532 void run_func_op(const vector<realt>& numbers, vector<int>::const_iterator &i,
00533                  realt*& stackPtr)
00534 {
00535     switch (*i) {
00536         //unary operators
00537         case OP_NEG:
00538             *stackPtr = - *stackPtr;
00539             break;
00540         case OP_SQRT:
00541             *stackPtr = sqrt(*stackPtr);
00542             break;
00543         case OP_EXP:
00544             *stackPtr = exp(*stackPtr);
00545             break;
00546         case OP_ERFC:
00547             *stackPtr = erfc(*stackPtr);
00548             break;
00549         case OP_ERF:
00550             *stackPtr = erf(*stackPtr);
00551             break;
00552         case OP_LOG10:
00553             *stackPtr = log10(*stackPtr);
00554             break;
00555         case OP_LN:
00556             *stackPtr = log(*stackPtr);
00557             break;
00558         case OP_SINH:
00559             *stackPtr = sinh(*stackPtr);
00560             break;
00561         case OP_COSH:
00562             *stackPtr = cosh(*stackPtr);
00563             break;
00564         case OP_TANH:
00565             *stackPtr = tanh(*stackPtr);
00566             break;
00567         case OP_SIN:
00568             *stackPtr = sin(*stackPtr);
00569             break;
00570         case OP_COS:
00571             *stackPtr = cos(*stackPtr);
00572             break;
00573         case OP_TAN:
00574             *stackPtr = tan(*stackPtr);
00575             break;
00576         case OP_ATAN:
00577             *stackPtr = atan(*stackPtr);
00578             break;
00579         case OP_ASIN:
00580             *stackPtr = asin(*stackPtr);
00581             break;
00582         case OP_ACOS:
00583             *stackPtr = acos(*stackPtr);
00584             break;
00585         case OP_LGAMMA:
00586             *stackPtr = boost::math::lgamma(*stackPtr);
00587             break;
00588         case OP_DIGAMMA:
00589             *stackPtr = boost::math::digamma(*stackPtr);
00590             break;
00591         case OP_ABS:
00592             *stackPtr = fabs(*stackPtr);
00593             break;
00594 
00595         //binary operators
00596         case OP_ADD:
00597             STACK_OFFSET_CHANGE(-1);
00598             *stackPtr += *(stackPtr+1);
00599             break;
00600         case OP_SUB:
00601             STACK_OFFSET_CHANGE(-1);
00602             *stackPtr -= *(stackPtr+1);
00603             break;
00604         case OP_MUL:
00605             STACK_OFFSET_CHANGE(-1);
00606             *stackPtr *= *(stackPtr+1);
00607             break;
00608         case OP_DIV:
00609             STACK_OFFSET_CHANGE(-1);
00610             *stackPtr /= *(stackPtr+1);
00611             break;
00612         case OP_POW:
00613             STACK_OFFSET_CHANGE(-1);
00614             *stackPtr = pow(*stackPtr, *(stackPtr+1));
00615             break;
00616         case OP_VOIGT:
00617             STACK_OFFSET_CHANGE(-1);
00618             *stackPtr = humlik(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00619             break;
00620         case OP_DVOIGT_DX:
00621             STACK_OFFSET_CHANGE(-1);
00622             *stackPtr = humdev_dkdx(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00623             break;
00624         case OP_DVOIGT_DY:
00625             STACK_OFFSET_CHANGE(-1);
00626             *stackPtr = humdev_dkdy(*stackPtr, *(stackPtr+1)) / sqrt(M_PI);
00627             break;
00628 
00629         // putting-number-to-stack-operators
00630         // stack overflow not checked
00631         case OP_NUMBER:
00632             STACK_OFFSET_CHANGE(+1);
00633             i++; // OP_NUMBER opcode is always followed by index
00634             *stackPtr = numbers[*i];
00635             break;
00636 
00637         case OP_TILDE:
00638             // used for default values in Runner::make_func_from_template()
00639             break;
00640 
00641         default:
00642             throw ExecuteError("op " + op2str(*i) +
00643                                " is not allowed for variables and functions");
00644     }
00645 }
00646 
00647 
00648 void ExprCalculator::transform_data(vector<Point>& points)
00649 {
00650     if (points.empty())
00651         return;
00652 
00653     realt stack[16];
00654     realt* stackPtr = stack - 1; // will be ++'ed first
00655     vector<Point> new_points = points;
00656 
00657     // do the time-consuming overflow checking only for the first point
00658     v_foreach (int, i, vm_.code()) {
00659         run_mutab_op(F_, vm_.numbers(), i, stackPtr, 0, points, new_points);
00660         if (stackPtr - stack >= 16)
00661             throw ExecuteError("stack overflow");
00662     }
00663     assert(stackPtr == stack - 1); // ASSIGN_ op must be at the end
00664 
00665     // the same for the rest of points, but without checks
00666     for (int n = 1; n != size(points); ++n)
00667         v_foreach (int, i, vm_.code())
00668             run_mutab_op(F_, vm_.numbers(), i, stackPtr, n, points, new_points);
00669 
00670     points.swap(new_points);
00671 }
00672 
00673 realt ExprCalculator::calculate(int n, const vector<Point>& points) const
00674 {
00675     realt stack[16];
00676     realt* stackPtr = stack - 1; // will be ++'ed first
00677     v_foreach (int, i, vm_.code()) {
00678         run_const_op(F_, vm_.numbers(), i, stackPtr, n, points, points);
00679         if (stackPtr - stack >= 16)
00680             throw ExecuteError("stack overflow");
00681     }
00682     //cerr << "stackPtr: " << stackPtr - stack << endl;
00683     assert(stackPtr == stack); // no ASSIGN_ at the end
00684     return stack[0];
00685 }
00686 
00687 realt ExprCalculator::calculate_custom(const vector<realt>& custom_val) const
00688 {
00689     realt stack[16];
00690     realt* stackPtr = stack - 1; // will be ++'ed first
00691     const vector<Point> dummy;
00692     v_foreach (int, i, vm_.code()) {
00693         if (*i == OP_SYMBOL) {
00694             STACK_OFFSET_CHANGE(+1);
00695             i++;
00696             if (is_index(*i, custom_val))
00697                 *stackPtr = custom_val[*i];
00698             else
00699                 throw ExecuteError("[internal] variable mismatch");
00700         }
00701         else
00702             run_const_op(F_, vm_.numbers(), i, stackPtr, 0, dummy, dummy);
00703         if (stackPtr - stack >= 16)
00704             throw ExecuteError("stack overflow");
00705     }
00706     assert(stackPtr == stack); // no ASSIGN_ at the end
00707     return stack[0];
00708 }
00709 
00710 /// executes VM code, sets derivatives and returns value
00711 realt run_code_for_variable(const VMData& vm,
00712                             const vector<Variable*> &variables,
00713                             vector<realt> &derivatives)
00714 {
00715     realt stack[16];
00716     realt* stackPtr = stack - 1; // will be ++'ed first
00717     v_foreach (int, i, vm.code()) {
00718         if (*i == OP_SYMBOL) {
00719             STACK_OFFSET_CHANGE(+1);
00720             ++i; // skip the next one
00721             *stackPtr = variables[*i]->get_value();
00722         }
00723         else if (*i == OP_PUT_DERIV) {
00724             ++i;
00725             // the OP_PUT_DERIV opcode is followed by a number n,
00726             // the derivative is calculated with respect to n'th variable
00727             assert(*i < (int) derivatives.size());
00728             derivatives[*i] = *stackPtr;
00729             STACK_OFFSET_CHANGE(-1);
00730         }
00731         else
00732             run_func_op(vm.numbers(), i, stackPtr);
00733     }
00734     assert(stackPtr == stack);
00735     return stack[0];
00736 }
00737 
00738 
00739 realt run_code_for_custom_func(const VMData& vm, realt x,
00740                                vector<realt> &derivatives)
00741 {
00742     realt stack[16];
00743     realt* stackPtr = stack - 1; // will be ++'ed first
00744     v_foreach (int, i, vm.code()) {
00745         if (*i == OP_X) {
00746             STACK_OFFSET_CHANGE(+1);
00747             *stackPtr = x;
00748         }
00749         else if (*i == OP_PUT_DERIV) {
00750             ++i;
00751             // the OP_PUT_DERIV opcode is followed by a number n,
00752             // the derivative is calculated with respect to n'th variable
00753             assert(*i < (int) derivatives.size());
00754             derivatives[*i] = *stackPtr;
00755             STACK_OFFSET_CHANGE(-1);
00756         }
00757         else
00758             run_func_op(vm.numbers(), i, stackPtr);
00759     }
00760     assert(stackPtr == stack);
00761     return stack[0];
00762 }
00763 
00764 realt run_code_for_custom_func_value(const VMData& vm, realt x,
00765                                      int code_offset)
00766 {
00767     realt stack[16];
00768     realt* stackPtr = stack - 1; // will be ++'ed first
00769     for (vector<int>::const_iterator i = vm.code().begin() + code_offset;
00770                                                  i != vm.code().end(); ++i) {
00771         if (*i == OP_X) {
00772             STACK_OFFSET_CHANGE(+1);
00773             *stackPtr = x;
00774         }
00775         else
00776             run_func_op(vm.numbers(), i, stackPtr);
00777     }
00778     assert(stackPtr == stack);
00779     return stack[0];
00780 }
00781