Branch data Line data Source code
1 : : /*
2 : : * differentiate.cpp - the Qucs equation differentiator implementations
3 : : *
4 : : * Copyright (C) 2007, 2008 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : : #include <ctype.h>
33 : : #include <cmath>
34 : :
35 : : #include "logging.h"
36 : : #include "complex.h"
37 : : #include "object.h"
38 : : #include "consts.h"
39 : : #include "equation.h"
40 : : #include "differentiate.h"
41 : :
42 : : using namespace qucs;
43 : : using namespace qucs::eqn;
44 : :
45 : : // Short helper macros.
46 : : #define C(con) ((constant *) (con))
47 : : #define A(con) ((application *) (con))
48 : : #define R(con) ((reference *) (con))
49 : : #define D(con) (C(con)->d)
50 : :
51 : : #define isConst(n) ((n)->getTag()==CONSTANT && C(n)->getType()==TAG_DOUBLE)
52 : : #define isRef(r,v) ((r)->getTag()==REFERENCE && !strcmp(R(r)->n,v))
53 : : #define isZero(n) (isConst(n) && D(n) == 0.0)
54 : : #define isOne(n) (isConst(n) && D(n) == 1.0)
55 : : #define isNeg(n) (isConst(n) && D(n) == -1.0)
56 : : #define isEuler(n) ((isConst(n) && D(n) == M_E) || isRef(n,"e"))
57 : : #define isval(n,v) (isConst(n) && D(n) == v)
58 : :
59 : : #define isVar(v) ((v)->getTag()==REFERENCE)
60 : : #define isApp(v) ((v)->getTag()==APPLICATION)
61 : : #define isMul(v) (isApp(v) && !strcmp(A(v)->n,"*"))
62 : : #define isSqr(v) (isApp(v) && !strcmp(A(v)->n,"sqr"))
63 : :
64 : : #define retCon(val) \
65 : : constant * res = new constant (TAG_DOUBLE); res->d = val; return res;
66 : : #define defCon(def,val) \
67 : : constant * def = new constant (TAG_DOUBLE); def->d = val;
68 : : #define defRef(def,var) \
69 : : reference * def = new reference (); def->n = strdup (var);
70 : : #define retApp1(op,f0) \
71 : : application * res = new application (); res->n = strdup (op); \
72 : : res->nargs = 1; res->args = f0; res->args->setNext (NULL); return res;
73 : : #define defApp1(def,op,f0) \
74 : : application * def = new application (); def->n = strdup (op); \
75 : : def->nargs = 1; def->args = f0; def->args->setNext (NULL);
76 : : #define defApp2(def,op,f0,f1) \
77 : : application * def = new application (); def->n = strdup (op); \
78 : : def->nargs = 2; def->args = f0; def->args->append (f1);
79 : : #define retApp2(op,f0,f1) \
80 : : application * res = new application (); res->n = strdup (op); \
81 : : res->nargs = 2; res->args = f0; res->args->append (f1); return res;
82 : : #define retApp3(op,f0,f1,f2) \
83 : : application * res = new application (); res->n = strdup (op); \
84 : : res->nargs = 3; res->args = f0; res->args->append (f1); \
85 : : res->args->append (f2); return res;
86 : : #define defApp3(def,op,f0,f1,f2) \
87 : : application * def = new application (); def->n = strdup (op); \
88 : : def->nargs = 3; def->args = f0; def->args->append (f1); \
89 : : def->args->append (f2);
90 : :
91 : : #define _A(idx) app->args->get(idx)
92 : : #define _A0 _A(0)
93 : : #define _A1 _A(1)
94 : : #define _A2 _A(2)
95 : : #define _D0 _A(0)->differentiate (derivative)
96 : : #define _D1 _A(1)->differentiate (derivative)
97 : : #define _D2 _A(2)->differentiate (derivative)
98 : :
99 : : #define _AF0(var) node * var = _A0;
100 : : #define _AF1(var) node * var = _A1;
101 : : #define _AF2(var) node * var = _A2;
102 : : #define _AD0(var) node * var = _D0;
103 : : #define _AD1(var) node * var = _D1;
104 : : #define _AD2(var) node * var = _D2;
105 : :
106 : : #define _AA(a,idx) A(a)->args->get(idx)
107 : : #define _AA0(a) _AA(a,0)
108 : : #define _AA1(a) _AA(a,1)
109 : :
110 : : #define _AAF0(a,var) node * var = _AA0(a);
111 : : #define _AAF1(a,var) node * var = _AA1(a);
112 : :
113 : 0 : node * differentiate::plus_binary (application * app, char * derivative) {
114 : 0 : _AD0 (d0);
115 : 0 : _AD1 (d1);
116 [ # # ][ # # ]: 0 : if (!isConst (d0) && !isConst (d1)) {
[ # # ][ # # ]
[ # # ]
117 [ # # ]: 0 : retApp2 ("+", d0, d1);
118 : : }
119 : 0 : return plus_reduce (d0, d1);
120 : : }
121 : :
122 : 0 : node * differentiate::plus_unary (application * app, char * derivative) {
123 : 0 : _AD0 (d0);
124 : 0 : return d0;
125 : : }
126 : :
127 : 0 : node * differentiate::plus_reduce (node * f0, node * f1) {
128 [ # # ][ # # ]: 0 : if (isZero (f0) && isZero (f1)) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
129 [ # # ][ # # ]: 0 : delete f0; delete f1;
130 [ # # ]: 0 : retCon (0);
131 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
132 [ # # ]: 0 : delete f0;
133 : 0 : return f1;
134 [ # # ][ # # ]: 0 : } else if (isZero (f1)) {
[ # # ][ # # ]
135 [ # # ]: 0 : delete f1;
136 : 0 : return f0;
137 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
138 : 0 : nr_double_t t = D(f0) + D(f1);
139 [ # # ][ # # ]: 0 : delete f0; delete f1;
140 [ # # ]: 0 : retCon (t);
141 : : } else {
142 [ # # ]: 0 : retApp2 ("+", f0, f1);
143 : : }
144 : : }
145 : :
146 : 0 : node * differentiate::minus_binary (application * app, char * derivative) {
147 : 0 : _AD0 (d0);
148 : 0 : _AD1 (d1);
149 [ # # ][ # # ]: 0 : if (!isConst (d0) && !isConst (d1)) {
[ # # ][ # # ]
[ # # ]
150 [ # # ]: 0 : retApp2 ("-", d0, d1);
151 : : }
152 : 0 : return minus_reduce (d0, d1);
153 : : }
154 : :
155 : 0 : node * differentiate::minus_unary (application * app, char * derivative) {
156 : 0 : _AD0 (d0);
157 : 0 : return minus_reduce (d0);
158 : : }
159 : :
160 : 0 : node * differentiate::minus_reduce (node * f0) {
161 [ # # ][ # # ]: 0 : if (isZero (f0)) {
[ # # ][ # # ]
162 [ # # ]: 0 : delete f0;
163 [ # # ]: 0 : retCon (0);
164 [ # # ][ # # ]: 0 : } else if (isConst (f0)) {
[ # # ]
165 : 0 : nr_double_t t = -D(f0);
166 [ # # ]: 0 : delete f0;
167 [ # # ]: 0 : retCon (t);
168 : : }
169 [ # # ]: 0 : retApp1 ("-", f0);
170 : : }
171 : :
172 : 0 : node * differentiate::minus_reduce (node * f0, node * f1) {
173 [ # # ][ # # ]: 0 : if (isZero (f0) && isZero (f1)) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
174 [ # # ][ # # ]: 0 : delete f0; delete f1;
175 [ # # ]: 0 : retCon (0);
176 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
177 [ # # ]: 0 : delete f0;
178 : 0 : return minus_reduce (f1);
179 [ # # ][ # # ]: 0 : } else if (isZero (f1)) {
[ # # ][ # # ]
180 [ # # ]: 0 : delete f1;
181 : 0 : return f0;
182 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
183 : 0 : nr_double_t t = D(f0) - D(f1);
184 [ # # ][ # # ]: 0 : delete f0; delete f1;
185 [ # # ]: 0 : retCon (t);
186 : : } else {
187 [ # # ]: 0 : retApp2 ("-", f0, f1);
188 : : }
189 : : }
190 : :
191 : 0 : node * differentiate::times (application * app, char * derivative) {
192 : 0 : _AF0 (f0);
193 : 0 : _AF1 (f1);
194 [ # # ][ # # ]: 0 : if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
195 [ # # ]: 0 : retCon (0);
196 : : }
197 : 0 : _AD0 (d0);
198 : 0 : _AD1 (d1);
199 : 0 : node * t1 = times_reduce (f0->recreate(), d1);
200 : 0 : node * t2 = times_reduce (f1->recreate(), d0);
201 : 0 : return plus_reduce (t1, t2);
202 : : }
203 : :
204 : 0 : node * differentiate::times_reduce (node * f0, node * f1) {
205 [ # # ][ # # ]: 0 : if (isZero (f0) || isZero (f1)) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
206 [ # # ][ # # ]: 0 : delete f0; delete f1;
207 [ # # ]: 0 : retCon (0);
208 [ # # ][ # # ]: 0 : } else if (isOne (f0)) {
[ # # ][ # # ]
209 [ # # ]: 0 : delete f0;
210 : 0 : return f1;
211 [ # # ][ # # ]: 0 : } else if (isNeg (f0)) {
[ # # ][ # # ]
212 [ # # ]: 0 : delete f0;
213 : 0 : return minus_reduce (f1);
214 [ # # ][ # # ]: 0 : } else if (isOne (f1)) {
[ # # ][ # # ]
215 [ # # ]: 0 : delete f1;
216 : 0 : return f0;
217 [ # # ][ # # ]: 0 : } else if (isNeg (f1)) {
[ # # ][ # # ]
218 [ # # ]: 0 : delete f1;
219 : 0 : return minus_reduce (f0);
220 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
221 : 0 : nr_double_t t = D(f0) * D(f1);
222 [ # # ][ # # ]: 0 : delete f0; delete f1;
223 [ # # ]: 0 : retCon (t);
224 : : } else {
225 [ # # ]: 0 : retApp2 ("*", f0, f1);
226 : : }
227 : : }
228 : :
229 : 0 : node * differentiate::over (application * app, char * derivative) {
230 : 0 : _AF0 (f0);
231 : 0 : _AF1 (f1);
232 [ # # ][ # # ]: 0 : if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
233 [ # # ]: 0 : retCon (0);
234 : : }
235 : 0 : _AD0 (d0);
236 : 0 : _AD1 (d1);
237 : 0 : node * t1 = times_reduce (f0->recreate(), d1);
238 : 0 : node * t2 = times_reduce (f1->recreate(), d0);
239 : 0 : node * t3 = minus_reduce (t2, t1);
240 : 0 : node * t4 = sqr_reduce (f1->recreate());
241 : 0 : return over_reduce (t3, t4);
242 : : }
243 : :
244 : 0 : node * differentiate::over_reduce (node * f0, node * f1) {
245 [ # # ][ # # ]: 0 : if (isOne (f0) && isOne (f1)) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
246 [ # # ][ # # ]: 0 : delete f0; delete f1;
247 [ # # ]: 0 : retCon (1);
248 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
249 [ # # ][ # # ]: 0 : delete f0; delete f1;
250 [ # # ]: 0 : retCon (0);
251 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
252 [ # # ][ # # ]: 0 : if (isZero (f1)) {
[ # # ][ # # ]
253 [ # # ]: 0 : retApp2 ("/", f0, f1);
254 : : }
255 : 0 : nr_double_t t = D(f0) / D(f1);
256 [ # # ][ # # ]: 0 : delete f0; delete f1;
257 [ # # ]: 0 : retCon (t);
258 [ # # ][ # # ]: 0 : } else if (isOne (f1)) {
[ # # ][ # # ]
259 [ # # ]: 0 : delete f1;
260 : 0 : return f0;
261 [ # # ][ # # ]: 0 : } else if (isNeg (f1)) {
[ # # ][ # # ]
262 [ # # ]: 0 : delete f1;
263 : 0 : return minus_reduce (f0);
264 : : } else {
265 : 0 : over_reduce_adv (f0, f1);
266 [ # # ]: 0 : retApp2 ("/", f0, f1);
267 : : }
268 : : }
269 : :
270 : 0 : void differentiate::over_reduce_adv (node * &f0, node * &f1) {
271 [ # # ]: 0 : if (isVar (f0)) {
272 [ # # ][ # # ]: 0 : if (isSqr (f1)) {
[ # # ]
273 : 0 : _AAF0 (f1,g1);
274 [ # # ]: 0 : if (isVar (g1)) {
275 [ # # ]: 0 : if (!strcmp (R(f0)->n, R(g1)->n)) {
276 [ # # ]: 0 : defCon (one, 1);
277 [ # # ]: 0 : reference * var = new reference (*R(g1));
278 [ # # ]: 0 : delete f0;
279 [ # # ]: 0 : delete f1;
280 : 0 : f0 = one;
281 : 0 : f1 = var;
282 : : }
283 : : }
284 : : }
285 : : }
286 : 0 : }
287 : :
288 : 0 : node * differentiate::power (application * app, char * derivative) {
289 : 0 : _AF0 (f0);
290 : 0 : _AF1 (f1);
291 [ # # ][ # # ]: 0 : if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
292 [ # # ]: 0 : retCon (0);
293 : : }
294 : 0 : _AD0 (d0);
295 : 0 : _AD1 (d1);
296 [ # # ][ # # ]: 0 : if (isZero (d1)) {
[ # # ][ # # ]
297 [ # # ]: 0 : defCon (one, 1);
298 : 0 : node * t1 = minus_reduce (f1->recreate(), one);
299 : 0 : node * t2 = power_reduce (f0->recreate(), t1);
300 : 0 : node * t3 = times_reduce (f1->recreate(), t2);
301 : 0 : return times_reduce (t3, d0);
302 : : }
303 : : else {
304 : 0 : node * t1 = power_reduce (f0->recreate(), f1->recreate());
305 : 0 : node * ln = ln_reduce (f0->recreate());
306 : 0 : node * t2 = times_reduce (d1, ln);
307 : 0 : node * t3 = times_reduce (f1->recreate(), d0);
308 : 0 : node * t4 = over_reduce (t3, f0->recreate());
309 : 0 : node * t5 = plus_reduce (t2, t4);
310 : 0 : return times_reduce (t1, t5);
311 : : }
312 : : }
313 : :
314 : 0 : node * differentiate::power_reduce (node * f0, node * f1) {
315 [ # # ][ # # ]: 0 : if (isOne (f0)) {
[ # # ][ # # ]
316 [ # # ][ # # ]: 0 : delete f0; delete f1;
317 [ # # ]: 0 : retCon (1);
318 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
319 [ # # ][ # # ]: 0 : delete f0; delete f1;
320 [ # # ]: 0 : retCon (0);
321 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
322 [ # # ][ # # ]: 0 : if (isZero (f1)) {
[ # # ][ # # ]
323 [ # # ][ # # ]: 0 : delete f0; delete f1;
324 [ # # ]: 0 : retCon (1);
325 : : }
326 : 0 : nr_double_t t = std::pow (D(f0), D(f1));
327 [ # # ][ # # ]: 0 : delete f0; delete f1;
328 [ # # ]: 0 : retCon (t);
329 [ # # ][ # # ]: 0 : } else if (isOne (f1)) {
[ # # ][ # # ]
330 [ # # ]: 0 : delete f1;
331 : 0 : return f0;
332 : : } else {
333 [ # # ]: 0 : retApp2 ("^", f0, f1);
334 : : }
335 : : }
336 : :
337 : 0 : node * differentiate::ln (application * app, char * derivative) {
338 : 0 : _AF0 (f0);
339 : 0 : _AD0 (d0);
340 : 0 : return over_reduce (d0, f0->recreate ());
341 : : }
342 : :
343 : 0 : node * differentiate::ln_reduce (node * f0) {
344 [ # # ][ # # ]: 0 : if (isOne (f0)) {
[ # # ][ # # ]
345 [ # # ]: 0 : delete f0;
346 [ # # ]: 0 : retCon (0);
347 [ # # ][ # # ]: 0 : } else if (isEuler (f0)) {
[ # # ][ # # ]
[ # # ][ # # ]
348 [ # # ]: 0 : delete f0;
349 [ # # ]: 0 : retCon (1);
350 : : }
351 [ # # ]: 0 : retApp1 ("ln", f0);
352 : : }
353 : :
354 : 0 : node * differentiate::log10 (application * app, char * derivative) {
355 : 0 : _AF0 (f0);
356 : 0 : _AD0 (d0);
357 : 0 : node * t1 = over_reduce (d0, f0->recreate ());
358 [ # # ]: 0 : defCon (ten, 10);
359 : 0 : return over_reduce (t1, ln_reduce (ten));
360 : : }
361 : :
362 : 0 : node * differentiate::log2 (application * app, char * derivative) {
363 : 0 : _AF0 (f0);
364 : 0 : _AD0 (d0);
365 : 0 : node * t1 = over_reduce (d0, f0->recreate ());
366 [ # # ]: 0 : defCon (two, 2);
367 : 0 : return over_reduce (t1, ln_reduce (two));
368 : : }
369 : :
370 : 0 : node * differentiate::sqrt (application * app, char * derivative) {
371 : 0 : _AF0 (f0);
372 : 0 : _AD0 (d0);
373 [ # # ]: 0 : defCon (half, 0.5);
374 : 0 : node * t1 = times_reduce (half, d0);
375 : 0 : node * t2 = sqrt_reduce (f0->recreate());
376 : 0 : return over_reduce (t1, t2);
377 : : }
378 : :
379 : 0 : node * differentiate::sqrt_reduce (node * f0) {
380 [ # # ][ # # ]: 0 : if (isOne (f0)) {
[ # # ][ # # ]
381 [ # # ]: 0 : delete f0;
382 [ # # ]: 0 : retCon (1);
383 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
384 [ # # ]: 0 : delete f0;
385 [ # # ]: 0 : retCon (0);
386 : : }
387 [ # # ]: 0 : retApp1 ("sqrt", f0);
388 : : }
389 : :
390 : 0 : node * differentiate::app_reduce (const char * func, node * d0, node * f0) {
391 [ # # ][ # # ]: 0 : if (isOne (d0)) {
[ # # ][ # # ]
392 [ # # ]: 0 : delete d0;
393 [ # # ]: 0 : retApp1 (func, f0);
394 [ # # ][ # # ]: 0 : } else if (isZero (d0)) {
[ # # ][ # # ]
395 [ # # ][ # # ]: 0 : delete d0; delete f0;
396 [ # # ]: 0 : retCon (0);
397 : : }
398 [ # # ]: 0 : defApp1 (app, func, f0);
399 : 0 : return times_reduce (d0, app);
400 : : }
401 : :
402 : 0 : node * differentiate::exp (application * app, char * derivative) {
403 : 0 : _AF0 (f0);
404 : 0 : _AD0 (d0);
405 : 0 : return app_reduce ("exp", d0, f0->recreate());
406 : : }
407 : :
408 : 0 : node * differentiate::limexp (application * app, char * derivative) {
409 : 0 : _AF0 (f0);
410 : 0 : _AD0 (d0);
411 [ # # ]: 0 : defCon (lexp, ::exp (M_LIMEXP));
412 [ # # ]: 0 : defCon (lcon, M_LIMEXP);
413 [ # # ]: 0 : defApp2 (ask, "<", f0->recreate(), lcon);
414 [ # # ]: 0 : defApp1 (exp, "exp", f0->recreate());
415 [ # # ]: 0 : defApp3 (ite, "?:", ask, exp, lexp);
416 : 0 : return times_reduce (d0, ite);
417 : : }
418 : :
419 : 0 : node * differentiate::sin (application * app, char * derivative) {
420 : 0 : _AF0 (f0);
421 : 0 : _AD0 (d0);
422 : 0 : return app_reduce ("cos", d0, f0->recreate());
423 : : }
424 : :
425 : 0 : node * differentiate::cos (application * app, char * derivative) {
426 : 0 : _AF0 (f0);
427 : 0 : _AD0 (d0);
428 : 0 : node * t1 = minus_reduce (d0);
429 : 0 : return app_reduce ("sin", t1, f0->recreate());
430 : : }
431 : :
432 : 0 : node * differentiate::tan (application * app, char * derivative) {
433 : 0 : _AF0 (f0);
434 : 0 : _AD0 (d0);
435 [ # # ]: 0 : defApp1 (sec, "sec", f0->recreate());
436 [ # # ]: 0 : defCon (two, 2);
437 : 0 : node * t1 = power_reduce (sec, two);
438 : 0 : return times_reduce (d0, t1);
439 : : }
440 : :
441 : 0 : node * differentiate::sec (application * app, char * derivative) {
442 : 0 : _AF0 (f0);
443 : 0 : _AD0 (d0);
444 [ # # ]: 0 : defApp1 (sin, "sin", f0->recreate());
445 [ # # ]: 0 : defApp1 (cos, "cos", f0->recreate());
446 [ # # ]: 0 : defCon (two, 2);
447 : 0 : node * t1 = power_reduce (cos, two);
448 : 0 : node * t2 = over_reduce (sin, t1);
449 : 0 : return times_reduce (d0, t2);
450 : : }
451 : :
452 : 0 : node * differentiate::cot (application * app, char * derivative) {
453 : 0 : _AF0 (f0);
454 : 0 : _AD0 (d0);
455 [ # # ]: 0 : defApp1 (cosec, "cosec", f0->recreate());
456 [ # # ]: 0 : defCon (two, 2);
457 : 0 : node * t1 = minus_reduce (d0);
458 : 0 : node * t2 = power_reduce (cosec, two);
459 : 0 : return times_reduce (t1, t2);
460 : : }
461 : :
462 : 0 : node * differentiate::cosec (application * app, char * derivative) {
463 : 0 : _AF0 (f0);
464 : 0 : _AD0 (d0);
465 [ # # ]: 0 : defApp1 (sin, "sin", f0->recreate());
466 [ # # ]: 0 : defApp1 (cos, "cos", f0->recreate());
467 [ # # ]: 0 : defCon (two, 2);
468 : 0 : node * t1 = minus_reduce (d0);
469 : 0 : node * t2 = power_reduce (sin, two);
470 : 0 : node * t3 = over_reduce (cos, t2);
471 : 0 : return times_reduce (t1, t3);
472 : : }
473 : :
474 : 0 : node * differentiate::abs (application * app, char * derivative) {
475 : 0 : _AF0 (f0);
476 : 0 : _AD0 (d0);
477 : 0 : return app_reduce ("sign", d0, f0->recreate());
478 : : }
479 : :
480 : 0 : node * differentiate::step (application *, char *) {
481 [ # # ]: 0 : retCon (0);
482 : : }
483 : :
484 : 0 : node * differentiate::sign (application *, char *) {
485 [ # # ]: 0 : retCon (0);
486 : : }
487 : :
488 : 0 : node * differentiate::arcsin (application * app, char * derivative) {
489 : 0 : _AF0 (f0);
490 : 0 : _AD0 (d0);
491 : 0 : node * sqr = sqr_reduce (f0->recreate());
492 [ # # ]: 0 : defCon (one, 1);
493 : 0 : node * t1 = minus_reduce (one, sqr);
494 : 0 : node * t2 = sqrt_reduce (t1);
495 : 0 : return over_reduce (d0, t2);
496 : : }
497 : :
498 : 0 : node * differentiate::square (application * app, char * derivative) {
499 : 0 : _AF0 (f0);
500 : 0 : _AD0 (d0);
501 [ # # ]: 0 : defCon (two, 2);
502 : 0 : node * t1 = times_reduce (two, d0);
503 : 0 : return times_reduce (t1, f0->recreate());
504 : : }
505 : :
506 : 0 : node * differentiate::sqr_reduce (node * f0) {
507 [ # # ][ # # ]: 0 : if (isOne (f0)) {
[ # # ][ # # ]
508 [ # # ]: 0 : delete f0;
509 [ # # ]: 0 : retCon (1);
510 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
511 [ # # ]: 0 : delete f0;
512 [ # # ]: 0 : retCon (0);
513 [ # # ][ # # ]: 0 : } else if (isConst (f0)) {
[ # # ]
514 : 0 : nr_double_t t = D(f0) * D(f0);
515 [ # # ]: 0 : delete f0;
516 [ # # ]: 0 : retCon (t);
517 : : } else {
518 [ # # ]: 0 : retApp1 ("sqr", f0);
519 : : }
520 : : }
521 : :
522 : 0 : node * differentiate::arccos (application * app, char * derivative) {
523 : 0 : _AF0 (f0);
524 : 0 : _AD0 (d0);
525 : 0 : node * sqr = sqr_reduce (f0->recreate());
526 [ # # ]: 0 : defCon (one, 1);
527 : 0 : node * t1 = minus_reduce (one, sqr);
528 : 0 : node * t2 = sqrt_reduce (t1);
529 : 0 : node * t3 = minus_reduce (d0);
530 : 0 : return over_reduce (t3, t2);
531 : : }
532 : :
533 : 0 : node * differentiate::arctan (application * app, char * derivative) {
534 : 0 : _AF0 (f0);
535 : 0 : _AD0 (d0);
536 : 0 : node * sqr = sqr_reduce (f0->recreate());
537 [ # # ]: 0 : defCon (one, 1);
538 : 0 : node * t1 = plus_reduce (one, sqr);
539 : 0 : return over_reduce (d0, t1);
540 : : }
541 : :
542 : 0 : node * differentiate::arccot (application * app, char * derivative) {
543 : 0 : _AF0 (f0);
544 : 0 : _AD0 (d0);
545 : 0 : node * sqr = sqr_reduce (f0->recreate());
546 [ # # ]: 0 : defCon (one, 1);
547 : 0 : node * t1 = plus_reduce (one, sqr);
548 : 0 : node * t2 = minus_reduce (d0);
549 : 0 : return over_reduce (t2, t1);
550 : : }
551 : :
552 : 0 : node * differentiate::arcsec (application * app, char * derivative) {
553 : 0 : _AF0 (f0);
554 : 0 : _AD0 (d0);
555 : 0 : node * sqr = sqr_reduce (f0->recreate());
556 [ # # ]: 0 : defCon (one, 1);
557 : 0 : node * t1 = minus_reduce (sqr, one);
558 : 0 : node * t2 = sqrt_reduce (t1);
559 : 0 : node * t3 = times_reduce (f0->recreate(), t2);
560 : 0 : return over_reduce (d0, t3);
561 : : }
562 : :
563 : 0 : node * differentiate::arccosec (application * app, char * derivative) {
564 : 0 : _AF0 (f0);
565 : 0 : _AD0 (d0);
566 : 0 : node * sqr = sqr_reduce (f0->recreate());
567 [ # # ]: 0 : defCon (one, 1);
568 : 0 : node * t1 = minus_reduce (sqr, one);
569 : 0 : node * t2 = sqrt_reduce (t1);
570 : 0 : node * t3 = times_reduce (f0->recreate(), t2);
571 : 0 : node * t4 = minus_reduce (d0);
572 : 0 : return over_reduce (t4, t3);
573 : : }
574 : :
575 : 0 : node * differentiate::sinh (application * app, char * derivative) {
576 : 0 : _AF0 (f0);
577 : 0 : _AD0 (d0);
578 : 0 : return app_reduce ("cosh", d0, f0->recreate());
579 : : }
580 : :
581 : 0 : node * differentiate::cosh (application * app, char * derivative) {
582 : 0 : _AF0 (f0);
583 : 0 : _AD0 (d0);
584 : 0 : return app_reduce ("sinh", d0, f0->recreate());
585 : : }
586 : :
587 : 0 : node * differentiate::tanh (application * app, char * derivative) {
588 : 0 : _AF0 (f0);
589 : 0 : _AD0 (d0);
590 [ # # ]: 0 : defApp1 (cosh, "cosh", f0->recreate());
591 [ # # ]: 0 : defCon (two, 2);
592 : 0 : node * t1 = power_reduce (cosh, two);
593 : 0 : return over_reduce (d0, t1);
594 : : }
595 : :
596 : 0 : node * differentiate::coth (application * app, char * derivative) {
597 : 0 : _AF0 (f0);
598 : 0 : _AD0 (d0);
599 [ # # ]: 0 : defApp1 (sinh, "sinh", f0->recreate());
600 [ # # ]: 0 : defCon (two, 2);
601 : 0 : node * t1 = power_reduce (sinh, two);
602 : 0 : node * t2 = minus_reduce (d0);
603 : 0 : return over_reduce (t2, t1);
604 : : }
605 : :
606 : 0 : node * differentiate::artanh (application * app, char * derivative) {
607 : 0 : _AF0 (f0);
608 : 0 : _AD0 (d0);
609 : 0 : node * sqr = sqr_reduce (f0->recreate());
610 [ # # ]: 0 : defCon (one, 1);
611 : 0 : node * t1 = minus_reduce (one, sqr);
612 : 0 : return over_reduce (d0, t1);
613 : : }
614 : :
615 : 0 : node * differentiate::arcoth (application * app, char * derivative) {
616 : 0 : _AF0 (f0);
617 : 0 : _AD0 (d0);
618 : 0 : node * sqr = sqr_reduce (f0->recreate());
619 [ # # ]: 0 : defCon (one, 1);
620 : 0 : node * t1 = minus_reduce (sqr, one);
621 : 0 : node * t2 = minus_reduce (d0);
622 : 0 : return over_reduce (t2, t1);
623 : : }
624 : :
625 : 0 : node * differentiate::arcosh (application * app, char * derivative) {
626 : 0 : _AF0 (f0);
627 : 0 : _AD0 (d0);
628 : 0 : node * sqr = sqr_reduce (f0->recreate());
629 [ # # ]: 0 : defCon (one, 1);
630 : 0 : node * t1 = minus_reduce (sqr, one);
631 : 0 : node * t2 = sqrt_reduce (t1);
632 : 0 : return over_reduce (d0, t2);
633 : : }
634 : :
635 : 0 : node * differentiate::arsinh (application * app, char * derivative) {
636 : 0 : _AF0 (f0);
637 : 0 : _AD0 (d0);
638 : 0 : node * sqr = sqr_reduce (f0->recreate());
639 [ # # ]: 0 : defCon (one, 1);
640 : 0 : node * t1 = plus_reduce (sqr, one);
641 : 0 : node * t2 = sqrt_reduce (t1);
642 : 0 : return over_reduce (d0, t2);
643 : : }
644 : :
645 : 0 : node * differentiate::arsech (application * app, char * derivative) {
646 : 0 : _AF0 (f0);
647 : 0 : _AD0 (d0);
648 : 0 : node * sqr = sqr_reduce (f0->recreate());
649 [ # # ]: 0 : defCon (one, 1);
650 : 0 : node * t1 = minus_reduce (one, sqr);
651 : 0 : node * t2 = sqrt_reduce (t1);
652 : 0 : node * t3 = times_reduce (f0->recreate(), t2);
653 : 0 : node * t4 = minus_reduce (d0);
654 : 0 : return over_reduce (t4, t3);
655 : : }
656 : :
657 : 0 : node * differentiate::arcosech (application * app, char * derivative) {
658 : 0 : _AF0 (f0);
659 : 0 : _AD0 (d0);
660 : 0 : node * sqr = sqr_reduce (f0->recreate());
661 [ # # ]: 0 : defCon (one, 1);
662 : 0 : node * t1 = plus_reduce (one, sqr);
663 : 0 : node * t2 = sqrt_reduce (t1);
664 : 0 : node * t3 = times_reduce (f0->recreate(), t2);
665 : 0 : node * t4 = minus_reduce (d0);
666 : 0 : return over_reduce (t4, t3);
667 : : }
668 : :
669 : 0 : node * differentiate::ifthenelse (application * app, char * derivative) {
670 : 0 : _AF0 (f0);
671 : 0 : _AD1 (d1);
672 : 0 : _AD2 (d2);
673 [ # # ][ # # ]: 0 : if (isConst (d1) && isConst (d2)) {
[ # # ][ # # ]
[ # # ]
674 [ # # ]: 0 : if (D(d1) == D(d2)) {
675 : 0 : nr_double_t t = D(d1);
676 [ # # ][ # # ]: 0 : delete d1; delete d2;
677 [ # # ]: 0 : retCon (t);
678 : : }
679 : : }
680 [ # # ]: 0 : retApp3 ("?:", f0->recreate(), d1, d2);
681 : : }
682 : :
683 : 0 : node * differentiate::sinc (application * app, char * derivative) {
684 : 0 : _AF0 (f0);
685 : 0 : _AD0 (d0);
686 [ # # ]: 0 : defApp1 (sinc, "sinc", f0->recreate());
687 [ # # ]: 0 : defApp1 (cos, "cos", f0->recreate());
688 : 0 : node * t1 = minus_reduce (cos, sinc);
689 : 0 : node * t2 = over_reduce (t1, f0->recreate());
690 : 0 : return times_reduce (d0, t2);
691 : : }
692 : :
693 : 0 : node * differentiate::norm (application * app, char * derivative) {
694 : 0 : _AF0 (f0);
695 : 0 : _AD0 (d0);
696 [ # # ]: 0 : defCon (two, 2);
697 : 0 : node * t1 = times_reduce (d0, two);
698 : 0 : return times_reduce (t1, f0->recreate());
699 : : }
700 : :
701 : 0 : node * differentiate::xhypot (application * app, char * derivative) {
702 : 0 : _AF0 (f0);
703 : 0 : _AF1 (f1);
704 : 0 : _AD0 (d0);
705 : 0 : _AD1 (d1);
706 : 0 : node * t1 = hypot_reduce (f0->recreate(), f1->recreate());
707 : 0 : node * t2 = times_reduce (d0, f0->recreate());
708 : 0 : node * t3 = times_reduce (d1, f1->recreate());
709 : 0 : node * t4 = plus_reduce (t2, t3);
710 : 0 : return over_reduce (t4, t1);
711 : : }
712 : :
713 : 0 : node * differentiate::hypot_reduce (node * f0, node * f1) {
714 [ # # ][ # # ]: 0 : if (isZero (f0) && isZero (f1)) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
715 [ # # ][ # # ]: 0 : delete f0; delete f1;
716 [ # # ]: 0 : retCon (0);
717 [ # # ][ # # ]: 0 : } else if (isZero (f0)) {
[ # # ][ # # ]
718 [ # # ]: 0 : delete f0;
719 : 0 : return sqrt_reduce (sqr_reduce (f1));
720 [ # # ][ # # ]: 0 : } else if (isZero (f1)) {
[ # # ][ # # ]
721 [ # # ]: 0 : delete f1;
722 : 0 : return sqrt_reduce (sqr_reduce (f0));
723 [ # # ][ # # ]: 0 : } else if (isConst (f0) && isConst (f1)) {
[ # # ][ # # ]
[ # # ]
724 : 0 : nr_double_t t = ::xhypot (D(f0), D(f1));
725 [ # # ][ # # ]: 0 : delete f0; delete f1;
726 [ # # ]: 0 : retCon (t);
727 : : } else {
728 [ # # ]: 0 : retApp2 ("hypot", f0, f1);
729 : : }
730 : : }
731 : :
732 : : #include "constants.h"
733 : :
734 : 0 : node * differentiate::vt (application * app, char * derivative) {
735 : 0 : _AD0 (d0);
736 [ # # ]: 0 : defCon (con, kBoverQ);
737 : 0 : return times_reduce (d0, con);
738 : : }
739 : :
740 : : // List of differentiators.
741 : : struct differentiation_t eqn::differentiations[] = {
742 : : { "+", differentiate::plus_binary, 2 },
743 : : { "+", differentiate::plus_unary, 1 },
744 : : { "-", differentiate::minus_binary, 2 },
745 : : { "-", differentiate::minus_unary, 1 },
746 : : { "*", differentiate::times, 2 },
747 : : { "/", differentiate::over, 2 },
748 : : { "^", differentiate::power, 2 },
749 : :
750 : : { "?:", differentiate::ifthenelse, 3 },
751 : :
752 : : { "ln", differentiate::ln, 1 },
753 : : { "log10", differentiate::log10, 1 },
754 : : { "log2", differentiate::log2, 1 },
755 : : { "sqrt", differentiate::sqrt, 1 },
756 : : { "exp", differentiate::exp, 1 },
757 : : { "sinc", differentiate::sinc, 1 },
758 : : { "norm", differentiate::norm, 1 },
759 : : { "sin", differentiate::sin, 1 },
760 : : { "cos", differentiate::cos, 1 },
761 : : { "tan", differentiate::tan, 1 },
762 : : { "sec", differentiate::sec, 1 },
763 : : { "cot", differentiate::cot, 1 },
764 : : { "cosec", differentiate::cosec, 1 },
765 : : { "abs", differentiate::abs, 1 },
766 : : { "step", differentiate::step, 1 },
767 : : { "sign", differentiate::sign, 1 },
768 : : { "arcsin", differentiate::arcsin, 1 },
769 : : { "arccos", differentiate::arccos, 1 },
770 : : { "arctan", differentiate::arctan, 1 },
771 : : { "arccot", differentiate::arccot, 1 },
772 : : { "arcsec", differentiate::arcsec, 1 },
773 : : { "arccosec", differentiate::arccosec, 1 },
774 : : { "sqr", differentiate::square, 1 },
775 : : { "sinh", differentiate::sinh, 1 },
776 : : { "cosh", differentiate::cosh, 1 },
777 : : { "tanh", differentiate::tanh, 1 },
778 : : { "coth", differentiate::coth, 1 },
779 : : { "arsinh", differentiate::arsinh, 1 },
780 : : { "arcosh", differentiate::arcosh, 1 },
781 : : { "artanh", differentiate::artanh, 1 },
782 : : { "arcoth", differentiate::arcoth, 1 },
783 : : { "arsech", differentiate::arsech, 1 },
784 : : { "arcosech", differentiate::arcosech, 1 },
785 : : { "hypot", differentiate::xhypot, 2 },
786 : : { "limexp", differentiate::limexp, 1 },
787 : : { "vt", differentiate::vt, 1 },
788 : :
789 : : { NULL, NULL, 0 /* end of list */ }
790 : : };
|
