45 #define C(con) ((constant *) (con))
46 #define A(con) ((application *) (con))
47 #define R(con) ((reference *) (con))
48 #define D(con) (C(con)->d)
50 #define isConst(n) ((n)->getTag()==CONSTANT && C(n)->getType()==TAG_DOUBLE)
51 #define isRef(r,v) ((r)->getTag()==REFERENCE && !strcmp(R(r)->n,v))
52 #define isZero(n) (isConst(n) && D(n) == 0.0)
53 #define isOne(n) (isConst(n) && D(n) == 1.0)
54 #define isNeg(n) (isConst(n) && D(n) == -1.0)
55 #define isEuler(n) ((isConst(n) && D(n) == M_E) || isRef(n,"e"))
56 #define isval(n,v) (isConst(n) && D(n) == v)
58 #define isVar(v) ((v)->getTag()==REFERENCE)
59 #define isApp(v) ((v)->getTag()==APPLICATION)
60 #define isMul(v) (isApp(v) && !strcmp(A(v)->n,"*"))
61 #define isSqr(v) (isApp(v) && !strcmp(A(v)->n,"sqr"))
64 constant * res = new constant (TAG_DOUBLE); res->d = val; return res;
65 #define defCon(def,val) \
66 constant * def = new constant (TAG_DOUBLE); def->d = val;
67 #define defRef(def,var) \
68 reference * def = new reference (); def->n = strdup (var);
69 #define retApp1(op,f0) \
70 application * res = new application (); res->n = strdup (op); \
71 res->nargs = 1; res->args = f0; res->args->setNext (NULL); return res;
72 #define defApp1(def,op,f0) \
73 application * def = new application (); def->n = strdup (op); \
74 def->nargs = 1; def->args = f0; def->args->setNext (NULL);
75 #define defApp2(def,op,f0,f1) \
76 application * def = new application (); def->n = strdup (op); \
77 def->nargs = 2; def->args = f0; def->args->append (f1);
78 #define retApp2(op,f0,f1) \
79 application * res = new application (); res->n = strdup (op); \
80 res->nargs = 2; res->args = f0; res->args->append (f1); return res;
81 #define retApp3(op,f0,f1,f2) \
82 application * res = new application (); res->n = strdup (op); \
83 res->nargs = 3; res->args = f0; res->args->append (f1); \
84 res->args->append (f2); return res;
85 #define defApp3(def,op,f0,f1,f2) \
86 application * def = new application (); def->n = strdup (op); \
87 def->nargs = 3; def->args = f0; def->args->append (f1); \
88 def->args->append (f2);
90 #define _A(idx) app->args->get(idx)
94 #define _D0 _A(0)->differentiate (derivative)
95 #define _D1 _A(1)->differentiate (derivative)
96 #define _D2 _A(2)->differentiate (derivative)
98 #define _AF0(var) node * var = _A0;
99 #define _AF1(var) node * var = _A1;
100 #define _AF2(var) node * var = _A2;
101 #define _AD0(var) node * var = _D0;
102 #define _AD1(var) node * var = _D1;
103 #define _AD2(var) node * var = _D2;
105 #define _AA(a,idx) A(a)->args->get(idx)
106 #define _AA0(a) _AA(a,0)
107 #define _AA1(a) _AA(a,1)
109 #define _AAF0(a,var) node * var = _AA0(a);
110 #define _AAF1(a,var) node * var = _AA1(a);
118 return plus_reduce (d0, d1);
126 node * differentiate::plus_reduce (
node * f0,
node * f1) {
128 delete f0;
delete f1;
137 nr_double_t
t =
D(f0) +
D(f1);
138 delete f0;
delete f1;
151 return minus_reduce (d0, d1);
156 return minus_reduce (d0);
159 node * differentiate::minus_reduce (
node * f0) {
164 nr_double_t
t = -
D(f0);
171 node * differentiate::minus_reduce (
node * f0,
node * f1) {
173 delete f0;
delete f1;
177 return minus_reduce (f1);
182 nr_double_t
t =
D(f0) -
D(f1);
183 delete f0;
delete f1;
198 node * t1 = times_reduce (f0->recreate(), d1);
199 node * t2 = times_reduce (f1->recreate(), d0);
200 return plus_reduce (t1, t2);
203 node * differentiate::times_reduce (
node * f0,
node * f1) {
205 delete f0;
delete f1;
207 }
else if (
isOne (f0)) {
210 }
else if (
isNeg (f0)) {
212 return minus_reduce (f1);
213 }
else if (
isOne (f1)) {
216 }
else if (
isNeg (f1)) {
218 return minus_reduce (f0);
220 nr_double_t
t =
D(f0) *
D(f1);
221 delete f0;
delete f1;
238 node * t3 = minus_reduce (t2, t1);
240 return over_reduce (t3, t4);
243 node * differentiate::over_reduce (
node * f0,
node * f1) {
245 delete f0;
delete f1;
248 delete f0;
delete f1;
254 nr_double_t
t =
D(f0) /
D(f1);
255 delete f0;
delete f1;
257 }
else if (
isOne (f1)) {
260 }
else if (
isNeg (f1)) {
262 return minus_reduce (f0);
264 over_reduce_adv (f0, f1);
269 void differentiate::over_reduce_adv (
node * &f0,
node * &f1) {
274 if (!strcmp (
R(f0)->
n,
R(g1)->
n)) {
276 reference * var =
new reference (*
R(g1));
297 node * t1 = minus_reduce (f1->recreate(), one);
298 node * t2 = power_reduce (f0->recreate(), t1);
299 node * t3 = times_reduce (f1->recreate(), t2);
300 return times_reduce (t3, d0);
303 node * t1 = power_reduce (f0->recreate(), f1->recreate());
304 node *
ln = ln_reduce (f0->recreate());
305 node * t2 = times_reduce (d1, ln);
306 node * t3 = times_reduce (f1->recreate(), d0);
307 node * t4 = over_reduce (t3, f0->recreate());
308 node * t5 = plus_reduce (t2, t4);
309 return times_reduce (t1, t5);
313 node * differentiate::power_reduce (
node * f0,
node * f1) {
315 delete f0;
delete f1;
318 delete f0;
delete f1;
322 delete f0;
delete f1;
325 nr_double_t
t =
pow (
D(f0),
D(f1));
326 delete f0;
delete f1;
328 }
else if (
isOne (f1)) {
339 return over_reduce (d0, f0->
recreate ());
342 node * differentiate::ln_reduce (
node * f0) {
358 return over_reduce (t1, ln_reduce (ten));
366 return over_reduce (t1, ln_reduce (two));
373 node * t1 = times_reduce (half, d0);
375 return over_reduce (t1, t2);
378 node * differentiate::sqrt_reduce (
node * f0) {
389 node * differentiate::app_reduce (
const char * func,
node * d0,
node * f0) {
394 delete d0;
delete f0;
398 return times_reduce (d0, app);
404 return app_reduce (
"exp", d0, f0->recreate());
412 defApp2 (ask,
"<", f0->recreate(), lcon);
415 return times_reduce (d0, ite);
421 return app_reduce (
"cos", d0, f0->recreate());
427 node * t1 = minus_reduce (d0);
428 return app_reduce (
"sin", t1, f0->recreate());
436 node * t1 = power_reduce (
sec, two);
437 return times_reduce (d0, t1);
446 node * t1 = power_reduce (
cos, two);
447 node * t2 = over_reduce (
sin, t1);
448 return times_reduce (d0, t2);
456 node * t1 = minus_reduce (d0);
458 return times_reduce (t1, t2);
467 node * t1 = minus_reduce (d0);
468 node * t2 = power_reduce (
sin, two);
469 node * t3 = over_reduce (
cos, t2);
470 return times_reduce (t1, t3);
476 return app_reduce (
"sign", d0, f0->recreate());
490 node *
sqr = sqr_reduce (f0->recreate());
492 node * t1 = minus_reduce (one, sqr);
493 node * t2 = sqrt_reduce (t1);
494 return over_reduce (d0, t2);
501 node * t1 = times_reduce (two, d0);
502 return times_reduce (t1, f0->recreate());
505 node * differentiate::sqr_reduce (
node * f0) {
513 nr_double_t t =
D(f0) *
D(f0);
526 node * t1 = minus_reduce (one, sqr);
527 node * t2 = sqrt_reduce (t1);
528 node * t3 = minus_reduce (d0);
529 return over_reduce (t3, t2);
537 node * t1 = plus_reduce (one, sqr);
538 return over_reduce (d0, t1);
546 node * t1 = plus_reduce (one, sqr);
547 node * t2 = minus_reduce (d0);
548 return over_reduce (t2, t1);
556 node * t1 = minus_reduce (sqr, one);
557 node * t2 = sqrt_reduce (t1);
559 return over_reduce (d0, t3);
567 node * t1 = minus_reduce (sqr, one);
568 node * t2 = sqrt_reduce (t1);
570 node * t4 = minus_reduce (d0);
571 return over_reduce (t4, t3);
577 return app_reduce (
"cosh", d0, f0->
recreate());
583 return app_reduce (
"sinh", d0, f0->
recreate());
591 node * t1 = power_reduce (
cosh, two);
592 return over_reduce (d0, t1);
600 node * t1 = power_reduce (
sinh, two);
601 node * t2 = minus_reduce (d0);
602 return over_reduce (t2, t1);
610 node * t1 = minus_reduce (one, sqr);
611 return over_reduce (d0, t1);
619 node * t1 = minus_reduce (sqr, one);
620 node * t2 = minus_reduce (d0);
621 return over_reduce (t2, t1);
629 node * t1 = minus_reduce (sqr, one);
630 node * t2 = sqrt_reduce (t1);
631 return over_reduce (d0, t2);
639 node * t1 = plus_reduce (sqr, one);
640 node * t2 = sqrt_reduce (t1);
641 return over_reduce (d0, t2);
649 node * t1 = minus_reduce (one, sqr);
650 node * t2 = sqrt_reduce (t1);
652 node * t4 = minus_reduce (d0);
653 return over_reduce (t4, t3);
661 node * t1 = plus_reduce (one, sqr);
662 node * t2 = sqrt_reduce (t1);
664 node * t4 = minus_reduce (d0);
665 return over_reduce (t4, t3);
673 if (
D(d1) ==
D(d2)) {
674 nr_double_t t =
D(d1);
675 delete d1;
delete d2;
689 return times_reduce (d0, t2);
696 node * t1 = times_reduce (d0, two);
697 return times_reduce (t1, f0->
recreate());
708 node * t4 = plus_reduce (t2, t3);
709 return over_reduce (t4, t1);
712 node * differentiate::hypot_reduce (
node * f0,
node * f1) {
714 delete f0;
delete f1;
718 return sqrt_reduce (sqr_reduce (f1));
721 return sqrt_reduce (sqr_reduce (f0));
724 delete f0;
delete f1;
736 return times_reduce (d0, con);