00001
00002
00003
00004 #include <boost/math/special_functions/gamma.hpp>
00005 #include <boost/math/special_functions/digamma.hpp>
00006
00007 #include <string>
00008 #include <vector>
00009 #include <cassert>
00010 #include <cstdlib>
00011 #include <cmath>
00012
00013 #include "common.h"
00014 #include "ui.h"
00015 #include "ast.h"
00016 #include "var.h"
00017 #include "logic.h"
00018 #include "numfuncs.h"
00019 #include "voigt.h"
00020 #include "lexer.h"
00021 #include "eparser.h"
00022
00023 using namespace std;
00024
00025
00026 string OpTree::str(const vector<string> *vars)
00027 {
00028 if (op < 0) {
00029 int v_nr = -op-1;
00030 return vars == NULL || vars->empty() ? "var"+S(v_nr) : (*vars)[v_nr];
00031 }
00032 switch (op) {
00033 case OP_NUMBER: return S(val);
00034 case OP_NEG: return "-" + c1->str_b(c1->op >= OP_POW, vars);
00035 case OP_EXP: return "exp(" + c1->str(vars) + ")";
00036 case OP_ERFC: return "erfc(" + c1->str(vars) + ")";
00037 case OP_ERF: return "erf(" + c1->str(vars) + ")";
00038 case OP_SINH: return "sinh(" + c1->str(vars) + ")";
00039 case OP_COSH: return "cosh(" + c1->str(vars) + ")";
00040 case OP_TANH: return "tanh("+ c1->str(vars) + ")";
00041 case OP_SIN: return "sin(" + c1->str(vars) + ")";
00042 case OP_COS: return "cos(" + c1->str(vars) + ")";
00043 case OP_TAN: return "tan(" + c1->str(vars) + ")";
00044 case OP_ASIN: return "asin("+ c1->str(vars) + ")";
00045 case OP_ACOS: return "acos("+ c1->str(vars) + ")";
00046 case OP_ATAN: return "atan("+ c1->str(vars) + ")";
00047 case OP_LGAMMA: return "lgamma("+ c1->str(vars) + ")";
00048 case OP_DIGAMMA: return "digamma("+ c1->str(vars) + ")";
00049 case OP_ABS: return "abs(" + c1->str(vars) + ")";
00050 case OP_LOG10:return "log10("+c1->str(vars) + ")";
00051 case OP_LN: return "ln(" + c1->str(vars) + ")";
00052 case OP_SQRT: return "sqrt("+ c1->str(vars) + ")";
00053 case OP_VOIGT: return "voigt("+ c1->str(vars) +","+ c2->str(vars) +")";
00054 case OP_DVOIGT_DX: return "dvoigt_dx("+ c1->str(vars)
00055 + "," + c2->str(vars) + ")";
00056 case OP_DVOIGT_DY: return "dvoigt_dy("+ c1->str(vars)
00057 + "," + c2->str(vars) + ")";
00058 case OP_POW: return c1->str_b(c1->op >= OP_POW, vars)
00059 + "^" + c2->str_b(c2->op >= OP_POW, vars);
00060 case OP_ADD: return c1->str(vars) + "+" + c2->str(vars);
00061 case OP_SUB: return c1->str(vars) + "-"
00062 + c2->str_b(c2->op >= OP_ADD, vars);
00063 case OP_MUL: return c1->str_b(c1->op >= OP_ADD, vars)
00064 + "*" + c2->str_b(c2->op >= OP_ADD, vars);
00065 case OP_DIV: return c1->str_b(c1->op >= OP_ADD, vars)
00066 + "/" + c2->str_b(c2->op >= OP_MUL, vars);
00067 default: assert(0); return "";
00068 }
00069 }
00070
00071
00072 OpTree* OpTree::clone() const
00073 {
00074 OpTree *t = new OpTree(*this);
00075 if (c1) t->c1 = c1->clone();
00076 if (c2) t->c2 = c2->clone();
00077 return t;
00078 }
00079
00080
00081
00082 OpTree* simplify_terms(OpTree *a);
00083 OpTree* do_multiply(OpTree *a, OpTree *b);
00084
00085 OpTree* do_neg(OpTree *a)
00086 {
00087 if (a->op == 0) {
00088 realt val = - a->val;
00089 delete a;
00090 return new OpTree(val);
00091 }
00092 else if (a->op == OP_NEG) {
00093 OpTree *t = a->c1->clone();
00094 delete a;
00095 return t;
00096 }
00097 else
00098 return new OpTree(OP_NEG, a);
00099 }
00100
00101 OpTree* do_add(int op, OpTree *a, OpTree *b)
00102 {
00103 if (a->op == 0 && b->op == 0) {
00104 realt val = (op == OP_ADD ? a->val + b->val : a->val - b->val);
00105 delete a;
00106 delete b;
00107 return new OpTree(val);
00108 }
00109 else if (a->op == 0 && is_eq(a->val, 0.)) {
00110 delete a;
00111 if (op == OP_ADD)
00112 return b;
00113 else
00114 return do_neg(b);
00115 }
00116 else if (b->op == 0 && is_eq(b->val, 0.)) {
00117 delete b;
00118 return a;
00119 }
00120 else if (b->op == OP_NEG) {
00121 OpTree *t = b->remove_c1();
00122 delete b;
00123 return do_add(op == OP_ADD ? OP_SUB : OP_ADD, a, t);
00124 }
00125 else if ((b->op == OP_MUL || b->op == OP_DIV)
00126 && b->c1->op == 0 && b->c1->val < 0) {
00127 b->c1->val = - b->c1->val;
00128 return do_add(op == OP_ADD ? OP_SUB : OP_ADD, a, b);
00129 }
00130 else if (*a == *b) {
00131 delete b;
00132 if (op == OP_ADD)
00133 return do_multiply(new OpTree(2.), a);
00134 else {
00135 delete a;
00136 return new OpTree(0.);
00137 }
00138 }
00139 else
00140 return new OpTree(op, a, b);
00141 }
00142
00143 OpTree* do_add(OpTree *a, OpTree *b)
00144 {
00145 return do_add(OP_ADD, a, b);
00146 }
00147
00148 OpTree* do_sub(OpTree *a, OpTree *b)
00149 {
00150 return do_add(OP_SUB, a, b);
00151 }
00152
00153 OpTree* do_multiply(OpTree *a, OpTree *b)
00154 {
00155 if (a->op == 0 && b->op == 0)
00156 {
00157 realt val = a->val * b->val;
00158 delete a;
00159 delete b;
00160 return new OpTree(val);
00161 }
00162 else if ((a->op == 0 && is_eq(a->val, 0.))
00163 || (b->op == 0 && is_eq(b->val, 0.)))
00164 {
00165 delete a;
00166 delete b;
00167 return new OpTree(0.);
00168 }
00169 else if (a->op == 0 && is_eq(a->val, 1.)) {
00170 delete a;
00171 return b;
00172 }
00173 else if (b->op == 0 && is_eq(b->val, 1.)) {
00174 delete b;
00175 return a;
00176 }
00177 else if (a->op == 0 && is_eq(a->val, -1.)) {
00178 delete a;
00179 return do_neg(b);
00180 }
00181 else if (b->op == 0 && is_eq(b->val, -1.)) {
00182 delete b;
00183 return do_neg(a);
00184 }
00185
00186 else if (a->op == 0 && b->op == OP_DIV && b->c1->op == 0) {
00187 b->c1->val *= a->val;
00188 delete a;
00189 return b;
00190 }
00191 else {
00192 return new OpTree(OP_MUL, a, b);
00193 }
00194 }
00195
00196 OpTree* do_divide(OpTree *a, OpTree *b)
00197 {
00198
00199 if (a->op == 0 && b->op == 0)
00200 {
00201 realt val = a->val / b->val;
00202 delete a;
00203 delete b;
00204 return new OpTree(val);
00205 }
00206 else if (a->op == 0 && is_eq(a->val, 0.)) {
00207 delete a;
00208 delete b;
00209 return new OpTree(0.);
00210 }
00211 else if (b->op == 0 && is_eq(b->val, 1.)) {
00212 delete b;
00213 return a;
00214 }
00215 else {
00216 return new OpTree(OP_DIV, a, b);
00217 }
00218 }
00219
00220 OpTree *do_sqr(OpTree *a)
00221 {
00222 return do_multiply(a, a->clone());
00223
00224 }
00225
00226 OpTree *do_oneover(OpTree *a)
00227 {
00228 return do_divide(new OpTree(1.), a);
00229 }
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 OpTree* do_exp(OpTree *a)
00246 {
00247 if (a->op == 0) {
00248 realt val = exp(a->val);
00249 delete a;
00250 return new OpTree(val);
00251 }
00252 else
00253 return new OpTree(OP_EXP, simplify_terms(a));
00254 }
00255
00256 OpTree* do_erf(OpTree *a)
00257 {
00258 if (a->op == 0) {
00259 realt val = erf(a->val);
00260 delete a;
00261 return new OpTree(val);
00262 }
00263 else
00264 return new OpTree(OP_ERF, simplify_terms(a));
00265 }
00266
00267 OpTree* do_erfc(OpTree *a)
00268 {
00269 if (a->op == 0) {
00270 realt val = erfc(a->val);
00271 delete a;
00272 return new OpTree(val);
00273 }
00274 else
00275 return new OpTree(OP_ERFC, simplify_terms(a));
00276 }
00277
00278 OpTree* do_sqrt(OpTree *a)
00279 {
00280 if (a->op == 0) {
00281 realt val = sqrt(a->val);
00282 delete a;
00283 return new OpTree(val);
00284 }
00285 else
00286 return new OpTree(OP_SQRT, a);
00287 }
00288
00289 OpTree* do_log10(OpTree *a)
00290 {
00291 if (a->op == 0) {
00292 realt val = log10(a->val);
00293 delete a;
00294 return new OpTree(val);
00295 }
00296 else
00297 return new OpTree(OP_LOG10, simplify_terms(a));
00298 }
00299
00300 OpTree* do_ln(OpTree *a)
00301 {
00302 if (a->op == 0) {
00303 realt val = log(a->val);
00304 delete a;
00305 return new OpTree(val);
00306 }
00307 else
00308 return new OpTree(OP_LN, simplify_terms(a));
00309 }
00310
00311 OpTree* do_sinh(OpTree *a)
00312 {
00313 if (a->op == 0) {
00314 realt val = sinh(a->val);
00315 delete a;
00316 return new OpTree(val);
00317 }
00318 else
00319 return new OpTree(OP_SINH, simplify_terms(a));
00320 }
00321
00322 OpTree* do_cosh(OpTree *a)
00323 {
00324 if (a->op == 0) {
00325 realt val = cosh(a->val);
00326 delete a;
00327 return new OpTree(val);
00328 }
00329 else
00330 return new OpTree(OP_COSH, simplify_terms(a));
00331 }
00332
00333 OpTree* do_tanh(OpTree *a)
00334 {
00335 if (a->op == 0) {
00336 realt val = tanh(a->val);
00337 delete a;
00338 return new OpTree(val);
00339 }
00340 else
00341 return new OpTree(OP_TANH, simplify_terms(a));
00342 }
00343
00344 OpTree* do_sin(OpTree *a)
00345 {
00346 if (a->op == 0) {
00347 realt val = sin(a->val);
00348 delete a;
00349 return new OpTree(val);
00350 }
00351 else
00352 return new OpTree(OP_SIN, simplify_terms(a));
00353 }
00354
00355 OpTree* do_cos(OpTree *a)
00356 {
00357 if (a->op == 0) {
00358 realt val = cos(a->val);
00359 delete a;
00360 return new OpTree(val);
00361 }
00362 else
00363 return new OpTree(OP_COS, simplify_terms(a));
00364 }
00365
00366 OpTree* do_tan(OpTree *a)
00367 {
00368 if (a->op == 0) {
00369 realt val = tan(a->val);
00370 delete a;
00371 return new OpTree(val);
00372 }
00373 else
00374 return new OpTree(OP_TAN, simplify_terms(a));
00375 }
00376
00377 OpTree* do_atan(OpTree *a)
00378 {
00379 if (a->op == 0) {
00380 realt val = atan(a->val);
00381 delete a;
00382 return new OpTree(val);
00383 }
00384 else
00385 return new OpTree(OP_ATAN, simplify_terms(a));
00386 }
00387
00388 OpTree* do_asin(OpTree *a)
00389 {
00390 if (a->op == 0) {
00391 realt val = asin(a->val);
00392 delete a;
00393 return new OpTree(val);
00394 }
00395 else
00396 return new OpTree(OP_ASIN, simplify_terms(a));
00397 }
00398
00399 OpTree* do_acos(OpTree *a)
00400 {
00401 if (a->op == 0) {
00402 realt val = acos(a->val);
00403 delete a;
00404 return new OpTree(val);
00405 }
00406 else
00407 return new OpTree(OP_ACOS, simplify_terms(a));
00408 }
00409
00410 OpTree* do_lgamma(OpTree *a)
00411 {
00412 if (a->op == 0) {
00413 realt val = boost::math::lgamma(a->val);
00414 delete a;
00415 return new OpTree(val);
00416 }
00417 else
00418 return new OpTree(OP_LGAMMA, simplify_terms(a));
00419 }
00420
00421 OpTree* do_digamma(OpTree *a)
00422 {
00423 if (a->op == 0) {
00424 realt val = boost::math::digamma(a->val);
00425 delete a;
00426 return new OpTree(val);
00427 }
00428 else
00429 return new OpTree(OP_DIGAMMA, simplify_terms(a));
00430 }
00431
00432 OpTree* do_abs(OpTree *a)
00433 {
00434 if (a->op == 0) {
00435 realt val = fabs(a->val);
00436 delete a;
00437 return new OpTree(val);
00438 }
00439 else
00440 return new OpTree(OP_ABS, simplify_terms(a));
00441 }
00442
00443 OpTree* do_pow(OpTree *a, OpTree *b)
00444 {
00445 if (a->op == 0 && b->op == 0) {
00446 realt val = pow(a->val, b->val);
00447 delete a;
00448 delete b;
00449 return new OpTree(val);
00450 }
00451 else if (a->op == 0 && is_eq(a->val, 0.)) {
00452 delete a;
00453 delete b;
00454 return new OpTree(0.);
00455 }
00456 else if ((b->op == 0 && is_eq(b->val, 0.))
00457 || (a->op == 0 && is_eq(a->val, 1.))) {
00458 delete a;
00459 delete b;
00460 return new OpTree(1.);
00461 }
00462 else if (b->op == 0 && is_eq(b->val, 1.)) {
00463 delete b;
00464 return a;
00465 }
00466 else if (b->op == 0 && is_eq(b->val, -1.)) {
00467 delete b;
00468 return do_oneover(a);
00469 }
00470 else {
00471 return new OpTree(OP_POW, a, simplify_terms(b));
00472 }
00473 }
00474
00475 OpTree* do_voigt(OpTree *a, OpTree *b)
00476 {
00477 if (a->op == 0 && b->op == 0) {
00478 realt val = humlik(a->val, b->val) / sqrt(M_PI);
00479 delete a;
00480 return new OpTree(val);
00481 }
00482 else
00483 return new OpTree(OP_VOIGT, simplify_terms(a), simplify_terms(b));
00484 }
00485
00486 OpTree* do_dvoigt_dx(OpTree *a, OpTree *b)
00487 {
00488 if (a->op == 0 && b->op == 0) {
00489 realt val = humdev_dkdx(a->val, b->val) / sqrt(M_PI);
00490 delete a;
00491 return new OpTree(val);
00492 }
00493 else
00494 return new OpTree(OP_DVOIGT_DX, simplify_terms(a), simplify_terms(b));
00495 }
00496
00497 OpTree* do_dvoigt_dy(OpTree *a, OpTree *b)
00498 {
00499 if (a->op == 0 && b->op == 0) {
00500 realt val = humdev_dkdy(a->val, b->val) / sqrt(M_PI);
00501 delete a;
00502 return new OpTree(val);
00503 }
00504 else
00505 return new OpTree(OP_DVOIGT_DY, simplify_terms(a), simplify_terms(b));
00506 }
00507
00508
00509
00510 struct MultFactor
00511 {
00512
00513 OpTree *t, *e;
00514 MultFactor(OpTree *t_, OpTree *e_) : t(t_), e(e_) {}
00515 void clear() { delete t; delete e; t=e=0; }
00516 };
00517
00518
00519
00520
00521
00522 void get_factors(OpTree *a, OpTree *expo,
00523 realt &constant, vector<MultFactor>& v)
00524 {
00525 if (a->op == OP_ADD || a->op == OP_SUB)
00526 a = simplify_terms(a);
00527 if (a->op == 0 && expo->op == 0)
00528 constant *= pow(a->val, expo->val);
00529 else if (a->op == OP_MUL) {
00530 get_factors(a->c1, expo, constant, v);
00531 get_factors(a->c2, expo, constant, v);
00532 }
00533 else if (a->op == OP_DIV) {
00534 get_factors(a->c1, expo, constant, v);
00535 OpTree *expo2 = do_neg(expo->clone());
00536 get_factors(a->c2, expo2, constant, v);
00537 delete expo2;
00538 }
00539 else if (a->op == OP_NEG) {
00540 get_factors(a->c1, expo, constant, v);
00541 get_factors(new OpTree(-1.), expo, constant, v);
00542 }
00543 else if (a->op == OP_SQRT) {
00544 OpTree *expo2 = do_multiply(new OpTree(0.5), expo->clone());
00545 get_factors(a->c1, expo2, constant, v);
00546 delete expo2;
00547 }
00548 else if (a->op == OP_POW) {
00549 OpTree *expo2 = do_multiply(a->remove_c2(), expo->clone());
00550 get_factors(a->c1, expo2, constant, v);
00551 delete expo2;
00552 }
00553 else {
00554 bool found = false;
00555 for (vector<MultFactor>::iterator i = v.begin(); i != v.end(); ++i)
00556 if (*i->t == *a) {
00557 i->e = do_add(i->e, expo->clone());
00558 found = true;
00559 break;
00560 }
00561 if (!found) {
00562 v.push_back(MultFactor(a, expo->clone()));
00563 return;
00564 }
00565 }
00566
00567 a->c1 = a->c2 = 0;
00568 delete a;
00569 }
00570
00571
00572 OpTree* simplify_factors(OpTree *a)
00573 {
00574 #ifdef DEBUG_SIMPLIFY
00575 cout << "simplify_factors() [<] " << a->str() << endl;
00576 #endif
00577 vector<MultFactor> v;
00578 OpTree expo(1.);
00579 realt constant = 1;
00580 get_factors(a, &expo, constant, v);
00581 #ifdef DEBUG_SIMPLIFY
00582 cout << "simplify_factors(): [.] {" << constant << "} ";
00583 for (vector<MultFactor>::iterator i = v.begin(); i != v.end(); ++i)
00584 cout << "{" << i->t->str() << "|" << i->e->str() << "} ";
00585 cout << endl;
00586 #endif
00587
00588
00589 for (vector<MultFactor>::iterator i = v.begin(); i != v.end(); ++i)
00590 if (i->t && i->t->op == OP_TAN) {
00591 for (vector<MultFactor>::iterator j = v.begin(); j != v.end(); ++j){
00592 if (j->t && j->t->op == OP_COS && *j->e == *i->e) {
00593 i->t->change_op(OP_SIN);
00594 j->clear();
00595 }
00596 if (j->t && j->t->op == OP_SIN
00597 && ((j->e->op==0 && i->e->op==0 && j->e->val==-i->e->val)
00598 || (j->e->op==OP_NEG && *j->e->c1 == *i->e)
00599 || (i->e->op==OP_NEG && *i->e->c1 == *j->e))) {
00600 i->t->change_op(OP_COS);
00601 j->clear();
00602 }
00603 }
00604 }
00605
00606 for (vector<MultFactor>::iterator i = v.begin(); i != v.end(); ++i)
00607 if (i->t && i->t->op == OP_SIN) {
00608 for (vector<MultFactor>::iterator j = v.begin(); j != v.end(); ++j){
00609 if (j->t && j->t->op == OP_COS
00610 && ((j->e->op==0 && i->e->op==0 && j->e->val==-i->e->val)
00611 || (j->e->op==OP_NEG && *j->e->c1 == *i->e)
00612 || (i->e->op==OP_NEG && *i->e->c1 == *j->e))) {
00613 i->t->change_op(OP_TAN);
00614 j->clear();
00615 }
00616 }
00617 }
00618
00619
00620
00621 OpTree *tu = 0, *tb = 0;
00622 for (vector<MultFactor>::iterator i = v.begin(); i != v.end(); ++i)
00623 if (i->t) {
00624 if ((i->e->op == 0 && i->e->val < 0) || i->e->op == OP_NEG) {
00625 OpTree *p = do_pow(i->t, do_neg(i->e));
00626 tb = (tb == 0 ? p : do_multiply(tb, p));
00627 }
00628 else {
00629 OpTree *p = do_pow(i->t, i->e);
00630 tu = (tu == 0 ? p : do_multiply(tu, p));
00631 }
00632 }
00633 OpTree *constant_t = new OpTree(constant);
00634 OpTree *ret = 0;
00635 if (tu) {
00636 if (tb)
00637 ret = do_multiply(constant_t, do_divide(tu, tb));
00638 else
00639 ret = do_multiply(constant_t, tu);
00640 }
00641 else {
00642 if (tb)
00643 ret = do_divide(constant_t, tb);
00644 else
00645 ret = constant_t;
00646 }
00647 #ifdef DEBUG_SIMPLIFY
00648 cout << "simplify_factors() [>] " << ret->str() << endl;
00649 #endif
00650 return ret;
00651 }
00652
00653
00654
00655 struct MultTerm
00656 {
00657 OpTree *t;
00658 realt k;
00659 MultTerm(OpTree *t_, realt k_) : t(t_), k(k_) {}
00660 void clear() { delete t; t=0; }
00661 };
00662
00663 void get_terms(OpTree *a, realt multiplier, vector<MultTerm> &v)
00664 {
00665 if (a->op == OP_MUL || a->op == OP_DIV || a->op == OP_SQRT
00666 || a->op == OP_POW)
00667 a = simplify_factors(a);
00668
00669 if (a->op == OP_ADD) {
00670 get_terms(a->c1, multiplier, v);
00671 get_terms(a->c2, multiplier, v);
00672 a->c1 = a->c2 = 0;
00673 delete a;
00674 }
00675 else if (a->op == OP_SUB) {
00676 get_terms(a->c1, multiplier, v);
00677 get_terms(a->c2, -multiplier, v);
00678 a->c1 = a->c2 = 0;
00679 delete a;
00680 }
00681 else if (a->op == OP_NEG) {
00682 get_terms(a->c1, -multiplier, v);
00683 a->c1 = a->c2 = 0;
00684 delete a;
00685 }
00686 else if (a->op == OP_MUL && a->c1->op == 0) {
00687 get_terms(a->c2, multiplier*(a->c1->val), v);
00688 a->c2 = 0;
00689 delete a;
00690 }
00691
00692 else if (a->op == OP_DIV && a->c1->op == 0 && a->c1->val != 1.) {
00693 get_terms(do_oneover(a->c2), multiplier*(a->c1->val), v);
00694 a->c2 = 0;
00695 delete a;
00696 }
00697 else {
00698 for (vector<MultTerm>::iterator i = v.begin(); i != v.end(); ++i) {
00699 if (a->op == 0 && i->t && i->t->op == 0) {
00700 i->k += multiplier * a->val;
00701 delete a;
00702 return;
00703 }
00704 if (i->t && *i->t == *a) {
00705 i->k += multiplier;
00706 delete a;
00707 return;
00708 }
00709 }
00710
00711 if (a->op == 0) {
00712 v.push_back(MultTerm(new OpTree(1.), multiplier * a->val));
00713 delete a;
00714 }
00715 else {
00716 v.push_back(MultTerm(a, multiplier));
00717 }
00718 }
00719 }
00720
00721 OpTree* simplify_terms(OpTree *a)
00722 {
00723
00724
00725
00726 if (a->op == OP_MUL || a->op == OP_DIV || a->op == OP_SQRT
00727 || a->op == OP_POW)
00728 return simplify_factors(a);
00729 else if (!(a->op == OP_ADD || a->op == OP_SUB || a->op == OP_NEG))
00730 return a;
00731 #ifdef DEBUG_SIMPLIFY
00732 cout << "simplify_terms() [<] " << a->str() << endl;
00733 #endif
00734 vector<MultTerm> v;
00735 get_terms(a, 1., v);
00736 #ifdef DEBUG_SIMPLIFY
00737 cout << "simplify_terms() [.] ";
00738 for (vector<MultTerm>::iterator i = v.begin(); i != v.end(); ++i)
00739 cout << "{" << i->t->str() << "|" << i->k << "} ";
00740 cout << endl;
00741 #endif
00742
00743
00744 realt to_add = 0.;
00745 for (vector<MultTerm>::iterator i = v.begin(); i != v.end(); ++i)
00746 if (i->t && i->t->op == OP_POW && i->t->c1->op == OP_SIN
00747 && i->t->c2->op == 0 && is_eq(i->t->c2->val, 2.))
00748 for (vector<MultTerm>::iterator j = v.begin(); j != v.end(); ++j)
00749 if (j->t && j->t->op == OP_POW && j->t->c1->op == OP_COS
00750 && j->t->c2->op == 0 && is_eq(j->t->c2->val, 2.)) {
00751 realt k = j->k;
00752 i->k -= k;
00753 j->clear();
00754 to_add += k;
00755 }
00756 if (to_add)
00757 get_terms(new OpTree(1.), to_add, v);
00758
00759
00760 OpTree *t = 0;
00761 for (vector<MultTerm>::iterator i = v.begin(); i != v.end(); ++i)
00762 if (i->t && !is_eq(i->k, 0)) {
00763 if (!t)
00764 t = do_multiply(new OpTree(i->k), i->t);
00765 else if (i->k > 0)
00766 t = do_add(t, do_multiply(new OpTree(i->k), i->t));
00767 else
00768 t = do_sub(t, do_multiply(new OpTree(-i->k), i->t));
00769 }
00770 if (!t)
00771 t = new OpTree(0.);
00772 #ifdef DEBUG_SIMPLIFY
00773 cout << "simplify_terms() [>] " << t->str() << endl;
00774 #endif
00775 return t;
00776 }
00777
00778
00779
00780
00781
00782
00783 string simplify_formula(string const &formula)
00784 {
00785 Lexer lex(formula.c_str());
00786 ExpressionParser ep(NULL);
00787
00788
00789 try {
00790 ep.parse_expr(lex, -1, NULL, NULL, ExpressionParser::kAstMode);
00791 }
00792 catch (fityk::SyntaxError&) {
00793 return formula;
00794 }
00795
00796 vector<OpTree*> trees = prepare_ast_with_der(ep.vm(), 1);
00797 vector<string> vars(1, "x");
00798 string simplified = trees.back()->str(&vars);
00799 purge_all_elements(trees);
00800 return simplified;
00801 }
00802
00803
00804 void get_derivatives_str(const char* formula, string& result)
00805 {
00806 Lexer lex(formula);
00807 ExpressionParser ep(NULL);
00808 vector<string> vars;
00809 ep.parse_expr(lex, -1, NULL, &vars);
00810 vector<OpTree*> trees = prepare_ast_with_der(ep.vm(), vars.size());
00811
00812 result += "f(" + join_vector(vars,", ") + ") = " + trees.back()->str(&vars);
00813 for (size_t i = 0; i != vars.size(); ++i)
00814 result += "\ndf / d " + vars[i] + " = " + trees[i]->str(&vars);
00815 purge_all_elements(trees);
00816 }
00817
00818
00819
00820 static
00821 vector<OpTree*> calculate_deriv(vector<int>::const_iterator &i, int len,
00822 const VMData& vm)
00823 {
00824 vector<OpTree*> results(len + 1, (OpTree*) NULL);
00825 --i;
00826 assert(i >= vm.code().begin());
00827
00828 if (*i < 0) {
00829 --i;
00830 assert(i >= vm.code().begin());
00831 assert(VMData::has_idx(*i));
00832 }
00833
00834 switch (*i) {
00835
00836 case OP_NUMBER: {
00837 for (int k = 0; k < len; ++k)
00838 results[k] = new OpTree(0.);
00839 int negative_idx = *(i+1);
00840 int idx = -1 - negative_idx;
00841 assert(is_index(idx, vm.numbers()));
00842 realt value = vm.numbers()[idx];
00843 results[len] = new OpTree(value);
00844 break;
00845 }
00846
00847 case OP_SYMBOL: {
00848 int negative_idx = *(i+1);
00849 int idx = -1 - negative_idx;
00850 assert(idx <= len);
00851 for (int k = 0; k < len; ++k) {
00852 realt der_value = (k == idx ? 1. : 0.);
00853 results[k] = new OpTree(der_value);
00854 }
00855 results[len] = new OpTree(NULL, idx);
00856 break;
00857 }
00858
00859 case OP_X: {
00860 int n = len - 1;
00861
00862 for (int k = 0; k < len; ++k) {
00863 realt der_value = (k == n ? 1. : 0.);
00864 results[k] = new OpTree(der_value);
00865 }
00866 results[len] = new OpTree(NULL, n);
00867 break;
00868 }
00869
00870
00871 case OP_SQRT:
00872 case OP_EXP:
00873 case OP_ERFC:
00874 case OP_ERF:
00875 case OP_LOG10:
00876 case OP_LN:
00877 case OP_SINH:
00878 case OP_COSH:
00879 case OP_TANH:
00880 case OP_SIN:
00881 case OP_COS:
00882 case OP_TAN:
00883 case OP_ASIN:
00884 case OP_ACOS:
00885 case OP_ATAN:
00886 case OP_LGAMMA:
00887 case OP_ABS:
00888 {
00889 int op = *i;
00890 vector<OpTree*> arg = calculate_deriv(i, len, vm);
00891 OpTree* (* do_op)(OpTree *) = NULL;
00892 OpTree* der = 0;
00893 OpTree* larg = arg.back()->clone();
00894 if (op == OP_SQRT) {
00895 der = do_divide(new OpTree(0.5), do_sqrt(larg));
00896 do_op = do_sqrt;
00897 }
00898 else if (op == OP_EXP) {
00899 der = do_exp(larg);
00900 do_op = do_exp;
00901 }
00902 else if (op == OP_ERFC) {
00903
00904 der = do_multiply(do_exp(do_neg(do_sqr(larg))),
00905 new OpTree(-2/sqrt(M_PI)));
00906 do_op = do_erfc;
00907 }
00908 else if (op == OP_ERF) {
00909
00910 der = do_multiply(do_exp(do_neg(do_sqr(larg))),
00911 new OpTree(2/sqrt(M_PI)));
00912 do_op = do_erf;
00913 }
00914 else if (op == OP_LOG10) {
00915 OpTree *ln_10 = do_ln(new OpTree(10.));
00916 der = do_oneover(do_multiply(larg, ln_10));
00917 do_op = do_log10;
00918 }
00919 else if (op == OP_LN) {
00920 der = do_oneover(larg);
00921 do_op = do_ln;
00922 }
00923 else if (op == OP_SINH) {
00924 der = do_cosh(larg);
00925 do_op = do_sinh;
00926 }
00927 else if (op == OP_COSH) {
00928 der = do_sinh(larg);
00929 do_op = do_cosh;
00930 }
00931 else if (op == OP_TANH) {
00932 der = do_oneover(do_sqr(do_cosh(larg)));
00933 do_op = do_tanh;
00934 }
00935 else if (op == OP_SIN) {
00936 der = do_cos(larg);
00937 do_op = do_sin;
00938 }
00939 else if (op == OP_COS) {
00940 der = do_neg(do_sin(larg));
00941 do_op = do_cos;
00942 }
00943 else if (op == OP_TAN) {
00944 der = do_oneover(do_sqr(do_cos(larg)));
00945 do_op = do_tan;
00946 }
00947 else if (op == OP_ASIN) {
00948 OpTree *root_arg = do_sub(new OpTree(1.), do_sqr(larg));
00949 der = do_oneover(do_sqrt(root_arg));
00950 do_op = do_asin;
00951 }
00952 else if (op == OP_ACOS) {
00953 OpTree *root_arg = do_sub(new OpTree(1.), do_sqr(larg));
00954 der = do_divide(new OpTree(-1.), do_sqrt(root_arg));
00955 do_op = do_acos;
00956 }
00957 else if (op == OP_ATAN) {
00958 der = do_oneover(do_add(new OpTree(1.), do_sqr(larg)));
00959 do_op = do_atan;
00960 }
00961 else if (op == OP_LGAMMA) {
00962 der = do_digamma(larg);
00963 do_op = do_lgamma;
00964 }
00965 else if (op == OP_ABS) {
00966 der = do_divide(do_abs(larg), larg->clone());
00967 do_op = do_abs;
00968 }
00969 else
00970 assert(0);
00971 for (int k = 0; k < len; ++k)
00972 results[k] = do_multiply(der->clone(), arg[k]);
00973 delete der;
00974 results[len] = (*do_op)(arg[len]);
00975 break;
00976 }
00977
00978 case OP_VOIGT:
00979 {
00980 int op = *i;
00981 vector<OpTree*> right = calculate_deriv(i, len, vm);
00982 vector<OpTree*> left = calculate_deriv(i, len, vm);
00983 assert (op == OP_VOIGT);
00984 OpTree *d1 = do_dvoigt_dx(left[len]->clone(), right[len]->clone());
00985 OpTree *d2 = do_dvoigt_dy(left[len]->clone(), right[len]->clone());
00986 for (int k = 0; k < len; ++k) {
00987 results[k] = do_add(do_multiply(d1->clone(), left[k]),
00988 do_multiply(d2->clone(), right[k]));
00989 }
00990 results[len] = do_voigt(left[len], right[len]);
00991 delete d1;
00992 delete d2;
00993 break;
00994 }
00995
00996 case OP_POW: {
00997 vector<OpTree*> right = calculate_deriv(i, len, vm);
00998 vector<OpTree*> left = calculate_deriv(i, len, vm);
00999
01000
01001 if (left[len]->op == 0 && right[len]->op == 0) {
01002 for (int k = 0; k < len; ++k) {
01003 assert(left[k]->op == 0 && right[k]->op == 0);
01004 delete left[k];
01005 delete right[k];
01006 results[k] = new OpTree(0.);
01007 }
01008 results[len] = do_pow(left[len], right[len]);
01009 }
01010
01011
01012 else {
01013 for (int k = 0; k < len; ++k) {
01014 OpTree *a = left[len],
01015 *b = right[len],
01016 *ap = left[k],
01017 *bp = right[k];
01018
01019
01020
01021 OpTree *pow_a_b = do_pow(a->clone(), b->clone());
01022 OpTree *term1 = do_divide(do_multiply(b->clone(), ap),
01023 a->clone());
01024 OpTree *term2 = do_multiply(do_ln(a->clone()), bp);
01025 results[k] = do_multiply(pow_a_b, do_add(term1, term2));
01026 }
01027 results[len] = do_pow(left[len], right[len]);
01028 }
01029 break;
01030 }
01031
01032 case OP_NEG: {
01033 vector<OpTree*> arg = calculate_deriv(i, len, vm);
01034 for (int k = 0; k < len+1; ++k)
01035 results[k] = do_neg(arg[k]);
01036 break;
01037 }
01038
01039 case OP_MUL: {
01040 vector<OpTree*> right = calculate_deriv(i, len, vm);
01041 vector<OpTree*> left = calculate_deriv(i, len, vm);
01042 for (int k = 0; k < len; ++k) {
01043
01044 OpTree *a = left[len],
01045 *b = right[len],
01046 *ap = left[k],
01047 *bp = right[k];
01048 results[k] = do_add(do_multiply(a->clone(), bp),
01049 do_multiply(ap, b->clone()));
01050 }
01051 results[len] = do_multiply(left[len], right[len]);
01052 break;
01053 }
01054
01055 case OP_DIV: {
01056 vector<OpTree*> right = calculate_deriv(i, len, vm);
01057 vector<OpTree*> left = calculate_deriv(i, len, vm);
01058 for (int k = 0; k < len; ++k) {
01059
01060 OpTree *a = left[len],
01061 *b = right[len],
01062 *ap = left[k],
01063 *bp = right[k];
01064 OpTree *upper = do_sub(do_multiply(ap, b->clone()),
01065 do_multiply(bp, a->clone()));
01066 results[k] = do_divide(upper, do_sqr(b->clone()));
01067 }
01068 results[len] = do_divide(left[len], right[len]);
01069 break;
01070 }
01071
01072 case OP_ADD: {
01073 vector<OpTree*> right = calculate_deriv(i, len, vm);
01074 vector<OpTree*> left = calculate_deriv(i, len, vm);
01075 for (int k = 0; k < len+1; ++k)
01076 results[k] = do_add(left[k], right[k]);
01077 break;
01078 }
01079
01080 case OP_SUB: {
01081 vector<OpTree*> right = calculate_deriv(i, len, vm);
01082 vector<OpTree*> left = calculate_deriv(i, len, vm);
01083 for (int k = 0; k < len+1; ++k)
01084 results[k] = do_sub(left[k], right[k]);
01085 break;
01086 }
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110 default:
01111 purge_all_elements(results);
01112 throw ExecuteError("`" + op2str(*i) + "' is not allowed for "
01113 "variables/functions");
01114 }
01115 assert(results[0] != NULL);
01116 assert(results[len] != NULL);
01117
01118 for (int k = 0; k < len+1; ++k)
01119 results[k] = simplify_terms(results[k]);
01120
01121
01122 return results;
01123 }
01124
01125 vector<OpTree*> prepare_ast_with_der(const VMData& vm, int len)
01126 {
01127 assert(!vm.code().empty());
01128 const_cast<VMData&>(vm).flip_indices();
01129 vector<int>::const_iterator iter = vm.code().end();
01130 vector<OpTree*> r = calculate_deriv(iter, len, vm);
01131 assert (iter == vm.code().begin());
01132 const_cast<VMData&>(vm).flip_indices();
01133 return r;
01134 }
01135
01136 void add_bytecode_from_tree(const OpTree* tree, const vector<int> &symbol_map,
01137 VMData& vm)
01138 {
01139 int op = tree->op;
01140 if (op < 0) {
01141 size_t n = -op-1;
01142 if (n == symbol_map.size()) {
01143 vm.append_code(OP_X);
01144 }
01145 else {
01146 assert(is_index(n, symbol_map));
01147 vm.append_code(OP_SYMBOL);
01148 vm.append_code(symbol_map[n]);
01149 }
01150 }
01151 else if (op == 0) {
01152 vm.append_number(tree->val);
01153 }
01154 else if (op >= OP_ONE_ARG && op < OP_TWO_ARG) {
01155 add_bytecode_from_tree(tree->c1, symbol_map, vm);
01156 vm.append_code(op);
01157 }
01158 else if (op >= OP_TWO_ARG) {
01159 add_bytecode_from_tree(tree->c1, symbol_map, vm);
01160 add_bytecode_from_tree(tree->c2, symbol_map, vm);
01161 vm.append_code(op);
01162 }
01163 }
01164