Qucs-core  0.0.19
differentiate.cpp
Go to the documentation of this file.
00001 /*
00002  * differentiate.cpp - the Qucs equation differentiator implementations
00003  *
00004  * Copyright (C) 2007, 2008 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include <ctype.h>
00033 #include <cmath>
00034 
00035 #include "logging.h"
00036 #include "complex.h"
00037 #include "object.h"
00038 #include "consts.h"
00039 #include "equation.h"
00040 #include "differentiate.h"
00041 
00042 using namespace qucs;
00043 using namespace qucs::eqn;
00044 
00045 // Short helper macros.
00046 #define C(con) ((constant *) (con))
00047 #define A(con) ((application *) (con))
00048 #define R(con) ((reference *) (con))
00049 #define D(con) (C(con)->d)
00050 
00051 #define isConst(n) ((n)->getTag()==CONSTANT && C(n)->getType()==TAG_DOUBLE)
00052 #define isRef(r,v) ((r)->getTag()==REFERENCE && !strcmp(R(r)->n,v))
00053 #define isZero(n)  (isConst(n) && D(n) == 0.0)
00054 #define isOne(n)   (isConst(n) && D(n) == 1.0)
00055 #define isNeg(n)   (isConst(n) && D(n) == -1.0)
00056 #define isEuler(n) ((isConst(n) && D(n) == euler) || isRef(n,"e"))
00057 #define isval(n,v) (isConst(n) && D(n) == v)
00058 
00059 #define isVar(v)   ((v)->getTag()==REFERENCE)
00060 #define isApp(v)   ((v)->getTag()==APPLICATION)
00061 #define isMul(v)   (isApp(v) && !strcmp(A(v)->n,"*"))
00062 #define isSqr(v)   (isApp(v) && !strcmp(A(v)->n,"sqr"))
00063 
00064 #define retCon(val) \
00065   constant * res = new constant (TAG_DOUBLE); res->d = val; return res;
00066 #define defCon(def,val) \
00067   constant * def = new constant (TAG_DOUBLE); def->d = val;
00068 #define defRef(def,var) \
00069   reference * def = new reference (); def->n = strdup (var);
00070 #define retApp1(op,f0) \
00071   application * res = new application (); res->n = strdup (op); \
00072   res->nargs = 1; res->args = f0; res->args->setNext (NULL); return res;
00073 #define defApp1(def,op,f0) \
00074   application * def = new application (); def->n = strdup (op); \
00075   def->nargs = 1; def->args = f0; def->args->setNext (NULL);
00076 #define defApp2(def,op,f0,f1) \
00077   application * def = new application (); def->n = strdup (op); \
00078   def->nargs = 2; def->args = f0; def->args->append (f1);
00079 #define retApp2(op,f0,f1) \
00080   application * res = new application (); res->n = strdup (op); \
00081   res->nargs = 2; res->args = f0; res->args->append (f1); return res;
00082 #define retApp3(op,f0,f1,f2) \
00083   application * res = new application (); res->n = strdup (op); \
00084   res->nargs = 3; res->args = f0; res->args->append (f1); \
00085   res->args->append (f2); return res;
00086 #define defApp3(def,op,f0,f1,f2) \
00087   application * def = new application (); def->n = strdup (op); \
00088   def->nargs = 3; def->args = f0; def->args->append (f1); \
00089   def->args->append (f2);
00090 
00091 #define _A(idx) app->args->get(idx)
00092 #define _A0 _A(0)
00093 #define _A1 _A(1)
00094 #define _A2 _A(2)
00095 #define _D0 _A(0)->differentiate (derivative)
00096 #define _D1 _A(1)->differentiate (derivative)
00097 #define _D2 _A(2)->differentiate (derivative)
00098 
00099 #define _AF0(var) node * var = _A0;
00100 #define _AF1(var) node * var = _A1;
00101 #define _AF2(var) node * var = _A2;
00102 #define _AD0(var) node * var = _D0;
00103 #define _AD1(var) node * var = _D1;
00104 #define _AD2(var) node * var = _D2;
00105 
00106 #define _AA(a,idx) A(a)->args->get(idx)
00107 #define _AA0(a) _AA(a,0)
00108 #define _AA1(a) _AA(a,1)
00109 
00110 #define _AAF0(a,var) node * var = _AA0(a);
00111 #define _AAF1(a,var) node * var = _AA1(a);
00112 
00113 node * differentiate::plus_binary (application * app, char * derivative) {
00114   _AD0 (d0);
00115   _AD1 (d1);
00116   if (!isConst (d0) && !isConst (d1)) {
00117     retApp2 ("+", d0, d1);
00118   }
00119   return plus_reduce (d0, d1);
00120 }
00121 
00122 node * differentiate::plus_unary (application * app, char * derivative) {
00123   _AD0 (d0);
00124   return d0;
00125 }
00126 
00127 node * differentiate::plus_reduce (node * f0, node * f1) {
00128   if (isZero (f0) && isZero (f1)) {
00129     delete f0; delete f1;
00130     retCon (0);
00131   } else if (isZero (f0)) {
00132     delete f0;
00133     return f1;
00134   } else if (isZero (f1)) {
00135     delete f1;
00136     return f0;
00137   } else if (isConst (f0) && isConst (f1)) {
00138     nr_double_t t = D(f0) + D(f1);
00139     delete f0; delete f1;
00140     retCon (t);
00141   } else {
00142     retApp2 ("+", f0, f1);
00143   }
00144 }
00145 
00146 node * differentiate::minus_binary (application * app, char * derivative) {
00147   _AD0 (d0);
00148   _AD1 (d1);
00149   if (!isConst (d0) && !isConst (d1)) {
00150     retApp2 ("-", d0, d1);
00151   }
00152   return minus_reduce (d0, d1);
00153 }
00154 
00155 node * differentiate::minus_unary (application * app, char * derivative) {
00156   _AD0 (d0);
00157   return minus_reduce (d0);
00158 }
00159 
00160 node * differentiate::minus_reduce (node * f0) {
00161   if (isZero (f0)) {
00162     delete f0;
00163     retCon (0);
00164   } else if (isConst (f0)) {
00165     nr_double_t t = -D(f0);
00166     delete f0;
00167     retCon (t);
00168   }
00169   retApp1 ("-", f0);
00170 }
00171 
00172 node * differentiate::minus_reduce (node * f0, node * f1) {
00173   if (isZero (f0) && isZero (f1)) {
00174     delete f0; delete f1;
00175     retCon (0);
00176   } else if (isZero (f0)) {
00177     delete f0;
00178     return minus_reduce (f1);
00179   } else if (isZero (f1)) {
00180     delete f1;
00181     return f0;
00182   } else if (isConst (f0) && isConst (f1)) {
00183     nr_double_t t = D(f0) - D(f1);
00184     delete f0; delete f1;
00185     retCon (t);
00186   } else {
00187     retApp2 ("-", f0, f1);
00188   }
00189 }
00190 
00191 node * differentiate::times (application * app, char * derivative) {
00192   _AF0 (f0);
00193   _AF1 (f1);
00194   if (isConst (f0) && isConst (f1)) {
00195     retCon (0);
00196   }
00197   _AD0 (d0);
00198   _AD1 (d1);
00199   node * t1 = times_reduce (f0->recreate(), d1);
00200   node * t2 = times_reduce (f1->recreate(), d0);
00201   return plus_reduce (t1, t2);
00202 }
00203 
00204 node * differentiate::times_reduce (node * f0, node * f1) {
00205   if (isZero (f0) || isZero (f1)) {
00206     delete f0; delete f1;
00207     retCon (0);
00208   } else if (isOne (f0)) {
00209     delete f0;
00210     return f1;
00211   } else if (isNeg (f0)) {
00212     delete f0;
00213     return minus_reduce (f1);
00214   } else if (isOne (f1)) {
00215     delete f1;
00216     return f0;
00217   } else if (isNeg (f1)) {
00218     delete f1;
00219     return minus_reduce (f0);
00220   } else if (isConst (f0) && isConst (f1)) {
00221     nr_double_t t = D(f0) * D(f1);
00222     delete f0; delete f1;
00223     retCon (t);
00224   } else {
00225     retApp2 ("*", f0, f1);
00226   }
00227 }
00228 
00229 node * differentiate::over (application * app, char * derivative) {
00230   _AF0 (f0);
00231   _AF1 (f1);
00232   if (isConst (f0) && isConst (f1)) {
00233     retCon (0);
00234   }
00235   _AD0 (d0);
00236   _AD1 (d1);
00237   node * t1 = times_reduce (f0->recreate(), d1);
00238   node * t2 = times_reduce (f1->recreate(), d0);
00239   node * t3 = minus_reduce (t2, t1);
00240   node * t4 = sqr_reduce (f1->recreate());
00241   return over_reduce (t3, t4);
00242 }
00243 
00244 node * differentiate::over_reduce (node * f0, node * f1) {
00245   if (isOne (f0) && isOne (f1)) {
00246     delete f0; delete f1;
00247     retCon (1);
00248   } else if (isZero (f0)) {
00249     delete f0; delete f1;
00250     retCon (0);
00251   } else if (isConst (f0) && isConst (f1)) {
00252     if (isZero (f1)) {
00253       retApp2 ("/", f0, f1);
00254     }
00255     nr_double_t t = D(f0) / D(f1);
00256     delete f0; delete f1;
00257     retCon (t);
00258   } else if (isOne (f1)) {
00259     delete f1;
00260     return f0;
00261   } else if (isNeg (f1)) {
00262     delete f1;
00263     return minus_reduce (f0);
00264   } else {
00265     over_reduce_adv (f0, f1);
00266     retApp2 ("/", f0, f1);
00267   }
00268 }
00269 
00270 void differentiate::over_reduce_adv (node * &f0, node * &f1) {
00271   if (isVar (f0)) {
00272     if (isSqr (f1)) {
00273       _AAF0 (f1,g1);
00274       if (isVar (g1)) {
00275         if (!strcmp (R(f0)->n, R(g1)->n)) {
00276           defCon (one, 1);
00277           reference * var = new reference (*R(g1));
00278           delete f0;
00279           delete f1;
00280           f0 = one;
00281           f1 = var;
00282         }
00283       }
00284     }
00285   }
00286 }
00287 
00288 node * differentiate::power (application * app, char * derivative) {
00289   _AF0 (f0);
00290   _AF1 (f1);
00291   if (isConst (f0) && isConst (f1)) {
00292     retCon (0);
00293   }
00294   _AD0 (d0);
00295   _AD1 (d1);
00296   if (isZero (d1)) {
00297     defCon (one, 1);
00298     node * t1 = minus_reduce (f1->recreate(), one);
00299     node * t2 = power_reduce (f0->recreate(), t1);
00300     node * t3 = times_reduce (f1->recreate(), t2);
00301     return times_reduce (t3, d0);
00302   }
00303   else {
00304     node * t1 = power_reduce (f0->recreate(), f1->recreate());
00305     node * ln = ln_reduce (f0->recreate());
00306     node * t2 = times_reduce (d1, ln);
00307     node * t3 = times_reduce (f1->recreate(), d0);
00308     node * t4 = over_reduce (t3, f0->recreate());
00309     node * t5 = plus_reduce (t2, t4);
00310     return times_reduce (t1, t5);
00311   }
00312 }
00313 
00314 node * differentiate::power_reduce (node * f0, node * f1) {
00315   if (isOne (f0)) {
00316     delete f0; delete f1;
00317     retCon (1);
00318   } else if (isZero (f0)) {
00319     delete f0; delete f1;
00320     retCon (0);
00321   } else if (isConst (f0) && isConst (f1)) {
00322     if (isZero (f1)) {
00323       delete f0; delete f1;
00324       retCon (1);
00325     }
00326     nr_double_t t = std::pow (D(f0), D(f1));
00327     delete f0; delete f1;
00328     retCon (t);
00329   } else if (isOne (f1)) {
00330     delete f1;
00331     return f0;
00332   } else {
00333     retApp2 ("^", f0, f1);
00334   }
00335 }
00336 
00337 node * differentiate::ln (application * app, char * derivative) {
00338   _AF0 (f0);
00339   _AD0 (d0);
00340   return over_reduce (d0, f0->recreate ());
00341 }
00342 
00343 node * differentiate::ln_reduce (node * f0) {
00344   if (isOne (f0)) {
00345     delete f0;
00346     retCon (0);
00347   } else if (isEuler (f0)) {
00348     delete f0;
00349     retCon (1);
00350   }
00351   retApp1 ("ln", f0);
00352 }
00353 
00354 node * differentiate::log10 (application * app, char * derivative) {
00355   _AF0 (f0);
00356   _AD0 (d0);
00357   node * t1 = over_reduce (d0, f0->recreate ());
00358   defCon (ten, 10);
00359   return over_reduce (t1, ln_reduce (ten));
00360 }
00361 
00362 node * differentiate::log2 (application * app, char * derivative) {
00363   _AF0 (f0);
00364   _AD0 (d0);
00365   node * t1 = over_reduce (d0, f0->recreate ());
00366   defCon (two, 2);
00367   return over_reduce (t1, ln_reduce (two));
00368 }
00369 
00370 node * differentiate::sqrt (application * app, char * derivative) {
00371   _AF0 (f0);
00372   _AD0 (d0);
00373   defCon (half, 0.5);
00374   node * t1 = times_reduce (half, d0);
00375   node * t2 = sqrt_reduce (f0->recreate());
00376   return over_reduce (t1, t2);
00377 }
00378 
00379 node * differentiate::sqrt_reduce (node * f0) {
00380   if (isOne (f0)) {
00381     delete f0;
00382     retCon (1);
00383   } else if (isZero (f0)) {
00384     delete f0;
00385     retCon (0);
00386   }
00387   retApp1 ("sqrt", f0);
00388 }
00389 
00390 node * differentiate::app_reduce (const char * func, node * d0, node * f0) {
00391   if (isOne (d0)) {
00392     delete d0;
00393     retApp1 (func, f0);
00394   } else if (isZero (d0)) {
00395     delete d0; delete f0;
00396     retCon (0);
00397   }
00398   defApp1 (app, func, f0);
00399   return times_reduce (d0, app);
00400 }
00401 
00402 node * differentiate::exp (application * app, char * derivative) {
00403   _AF0 (f0);
00404   _AD0 (d0);
00405   return app_reduce ("exp", d0, f0->recreate());
00406 }
00407 
00408 node * differentiate::limexp (application * app, char * derivative) {
00409   _AF0 (f0);
00410   _AD0 (d0);
00411   defCon (lexp, ::exp (limitexp));
00412   defCon (lcon, limitexp);
00413   defApp2 (ask, "<", f0->recreate(), lcon);
00414   defApp1 (exp, "exp", f0->recreate());
00415   defApp3 (ite, "?:", ask, exp, lexp);
00416   return times_reduce (d0, ite);
00417 }
00418 
00419 node * differentiate::sin (application * app, char * derivative) {
00420   _AF0 (f0);
00421   _AD0 (d0);
00422   return app_reduce ("cos", d0, f0->recreate());
00423 }
00424 
00425 node * differentiate::cos (application * app, char * derivative) {
00426   _AF0 (f0);
00427   _AD0 (d0);
00428   node * t1 = minus_reduce (d0);
00429   return app_reduce ("sin", t1, f0->recreate());
00430 }
00431 
00432 node * differentiate::tan (application * app, char * derivative) {
00433   _AF0 (f0);
00434   _AD0 (d0);
00435   defApp1 (sec, "sec", f0->recreate());
00436   defCon (two, 2);
00437   node * t1 = power_reduce (sec, two);
00438   return times_reduce (d0, t1);
00439 }
00440 
00441 node * differentiate::sec (application * app, char * derivative) {
00442   _AF0 (f0);
00443   _AD0 (d0);
00444   defApp1 (sin, "sin", f0->recreate());
00445   defApp1 (cos, "cos", f0->recreate());
00446   defCon (two, 2);
00447   node * t1 = power_reduce (cos, two);
00448   node * t2 = over_reduce (sin, t1);
00449   return times_reduce (d0, t2);
00450 }
00451 
00452 node * differentiate::cot (application * app, char * derivative) {
00453   _AF0 (f0);
00454   _AD0 (d0);
00455   defApp1 (cosec, "cosec", f0->recreate());
00456   defCon (two, 2);
00457   node * t1 = minus_reduce (d0);
00458   node * t2 = power_reduce (cosec, two);
00459   return times_reduce (t1, t2);
00460 }
00461 
00462 node * differentiate::cosec (application * app, char * derivative) {
00463   _AF0 (f0);
00464   _AD0 (d0);
00465   defApp1 (sin, "sin", f0->recreate());
00466   defApp1 (cos, "cos", f0->recreate());
00467   defCon (two, 2);
00468   node * t1 = minus_reduce (d0);
00469   node * t2 = power_reduce (sin, two);
00470   node * t3 = over_reduce (cos, t2);
00471   return times_reduce (t1, t3);
00472 }
00473 
00474 node * differentiate::abs (application * app, char * derivative) {
00475   _AF0 (f0);
00476   _AD0 (d0);
00477   return app_reduce ("sign", d0, f0->recreate());
00478 }
00479 
00480 node * differentiate::step (application *, char *) {
00481   retCon (0);
00482 }
00483 
00484 node * differentiate::sign (application *, char *) {
00485   retCon (0);
00486 }
00487 
00488 node * differentiate::arcsin (application * app, char * derivative) {
00489   _AF0 (f0);
00490   _AD0 (d0);
00491   node * sqr = sqr_reduce (f0->recreate());
00492   defCon (one, 1);
00493   node * t1 = minus_reduce (one, sqr);
00494   node * t2 = sqrt_reduce (t1);
00495   return over_reduce (d0, t2);
00496 }
00497 
00498 node * differentiate::square (application * app, char * derivative) {
00499   _AF0 (f0);
00500   _AD0 (d0);
00501   defCon (two, 2);
00502   node * t1 = times_reduce (two, d0);
00503   return times_reduce (t1, f0->recreate());
00504 }
00505 
00506 node * differentiate::sqr_reduce (node * f0) {
00507   if (isOne (f0)) {
00508     delete f0;
00509     retCon (1);
00510   } else if (isZero (f0)) {
00511     delete f0;
00512     retCon (0);
00513   } else if (isConst (f0)) {
00514     nr_double_t t = D(f0) * D(f0);
00515     delete f0;
00516     retCon (t);
00517   } else {
00518     retApp1 ("sqr", f0);
00519   }
00520 }
00521 
00522 node * differentiate::arccos (application * app, char * derivative) {
00523   _AF0 (f0);
00524   _AD0 (d0);
00525   node * sqr = sqr_reduce (f0->recreate());
00526   defCon (one, 1);
00527   node * t1 = minus_reduce (one, sqr);
00528   node * t2 = sqrt_reduce (t1);
00529   node * t3 = minus_reduce (d0);
00530   return over_reduce (t3, t2);
00531 }
00532 
00533 node * differentiate::arctan (application * app, char * derivative) {
00534   _AF0 (f0);
00535   _AD0 (d0);
00536   node * sqr = sqr_reduce (f0->recreate());
00537   defCon (one, 1);
00538   node * t1 = plus_reduce (one, sqr);
00539   return over_reduce (d0, t1);
00540 }
00541 
00542 node * differentiate::arccot (application * app, char * derivative) {
00543   _AF0 (f0);
00544   _AD0 (d0);
00545   node * sqr = sqr_reduce (f0->recreate());
00546   defCon (one, 1);
00547   node * t1 = plus_reduce (one, sqr);
00548   node * t2 = minus_reduce (d0);
00549   return over_reduce (t2, t1);
00550 }
00551 
00552 node * differentiate::arcsec (application * app, char * derivative) {
00553   _AF0 (f0);
00554   _AD0 (d0);
00555   node * sqr = sqr_reduce (f0->recreate());
00556   defCon (one, 1);
00557   node * t1 = minus_reduce (sqr, one);
00558   node * t2 = sqrt_reduce (t1);
00559   node * t3 = times_reduce (f0->recreate(), t2);
00560   return over_reduce (d0, t3);
00561 }
00562 
00563 node * differentiate::arccosec (application * app, char * derivative) {
00564   _AF0 (f0);
00565   _AD0 (d0);
00566   node * sqr = sqr_reduce (f0->recreate());
00567   defCon (one, 1);
00568   node * t1 = minus_reduce (sqr, one);
00569   node * t2 = sqrt_reduce (t1);
00570   node * t3 = times_reduce (f0->recreate(), t2);
00571   node * t4 = minus_reduce (d0);
00572   return over_reduce (t4, t3);
00573 }
00574 
00575 node * differentiate::sinh (application * app, char * derivative) {
00576   _AF0 (f0);
00577   _AD0 (d0);
00578   return app_reduce ("cosh", d0, f0->recreate());
00579 }
00580 
00581 node * differentiate::cosh (application * app, char * derivative) {
00582   _AF0 (f0);
00583   _AD0 (d0);
00584   return app_reduce ("sinh", d0, f0->recreate());
00585 }
00586 
00587 node * differentiate::tanh (application * app, char * derivative) {
00588   _AF0 (f0);
00589   _AD0 (d0);
00590   defApp1 (cosh, "cosh", f0->recreate());
00591   defCon (two, 2);
00592   node * t1 = power_reduce (cosh, two);
00593   return over_reduce (d0, t1);
00594 }
00595 
00596 node * differentiate::coth (application * app, char * derivative) {
00597   _AF0 (f0);
00598   _AD0 (d0);
00599   defApp1 (sinh, "sinh", f0->recreate());
00600   defCon (two, 2);
00601   node * t1 = power_reduce (sinh, two);
00602   node * t2 = minus_reduce (d0);
00603   return over_reduce (t2, t1);
00604 }
00605 
00606 node * differentiate::artanh (application * app, char * derivative) {
00607   _AF0 (f0);
00608   _AD0 (d0);
00609   node * sqr = sqr_reduce (f0->recreate());
00610   defCon (one, 1);
00611   node * t1 = minus_reduce (one, sqr);
00612   return over_reduce (d0, t1);
00613 }
00614 
00615 node * differentiate::arcoth (application * app, char * derivative) {
00616   _AF0 (f0);
00617   _AD0 (d0);
00618   node * sqr = sqr_reduce (f0->recreate());
00619   defCon (one, 1);
00620   node * t1 = minus_reduce (sqr, one);
00621   node * t2 = minus_reduce (d0);
00622   return over_reduce (t2, t1);
00623 }
00624 
00625 node * differentiate::arcosh (application * app, char * derivative) {
00626   _AF0 (f0);
00627   _AD0 (d0);
00628   node * sqr = sqr_reduce (f0->recreate());
00629   defCon (one, 1);
00630   node * t1 = minus_reduce (sqr, one);
00631   node * t2 = sqrt_reduce (t1);
00632   return over_reduce (d0, t2);
00633 }
00634 
00635 node * differentiate::arsinh (application * app, char * derivative) {
00636   _AF0 (f0);
00637   _AD0 (d0);
00638   node * sqr = sqr_reduce (f0->recreate());
00639   defCon (one, 1);
00640   node * t1 = plus_reduce (sqr, one);
00641   node * t2 = sqrt_reduce (t1);
00642   return over_reduce (d0, t2);
00643 }
00644 
00645 node * differentiate::arsech (application * app, char * derivative) {
00646   _AF0 (f0);
00647   _AD0 (d0);
00648   node * sqr = sqr_reduce (f0->recreate());
00649   defCon (one, 1);
00650   node * t1 = minus_reduce (one, sqr);
00651   node * t2 = sqrt_reduce (t1);
00652   node * t3 = times_reduce (f0->recreate(), t2);
00653   node * t4 = minus_reduce (d0);
00654   return over_reduce (t4, t3);
00655 }
00656 
00657 node * differentiate::arcosech (application * app, char * derivative) {
00658   _AF0 (f0);
00659   _AD0 (d0);
00660   node * sqr = sqr_reduce (f0->recreate());
00661   defCon (one, 1);
00662   node * t1 = plus_reduce (one, sqr);
00663   node * t2 = sqrt_reduce (t1);
00664   node * t3 = times_reduce (f0->recreate(), t2);
00665   node * t4 = minus_reduce (d0);
00666   return over_reduce (t4, t3);
00667 }
00668 
00669 node * differentiate::ifthenelse (application * app, char * derivative) {
00670   _AF0 (f0);
00671   _AD1 (d1);
00672   _AD2 (d2);
00673   if (isConst (d1) && isConst (d2)) {
00674     if (D(d1) == D(d2)) {
00675       nr_double_t t = D(d1);
00676       delete d1; delete d2;
00677       retCon (t);
00678     }
00679   }
00680   retApp3 ("?:", f0->recreate(), d1, d2);
00681 }
00682 
00683 node * differentiate::sinc (application * app, char * derivative) {
00684   _AF0 (f0);
00685   _AD0 (d0);
00686   defApp1 (sinc, "sinc", f0->recreate());
00687   defApp1 (cos, "cos", f0->recreate());
00688   node * t1 = minus_reduce (cos, sinc);
00689   node * t2 = over_reduce (t1, f0->recreate());
00690   return times_reduce (d0, t2);
00691 }
00692 
00693 node * differentiate::norm (application * app, char * derivative) {
00694   _AF0 (f0);
00695   _AD0 (d0);
00696   defCon (two, 2);
00697   node * t1 = times_reduce (d0, two);
00698   return times_reduce (t1, f0->recreate());
00699 }
00700 
00701 node * differentiate::xhypot (application * app, char * derivative) {
00702   _AF0 (f0);
00703   _AF1 (f1);
00704   _AD0 (d0);
00705   _AD1 (d1);
00706   node * t1 = hypot_reduce (f0->recreate(), f1->recreate());
00707   node * t2 = times_reduce (d0, f0->recreate());
00708   node * t3 = times_reduce (d1, f1->recreate());
00709   node * t4 = plus_reduce (t2, t3);
00710   return over_reduce (t4, t1);
00711 }
00712 
00713 node * differentiate::hypot_reduce (node * f0, node * f1) {
00714   if (isZero (f0) && isZero (f1)) {
00715     delete f0; delete f1;
00716     retCon (0);
00717   } else if (isZero (f0)) {
00718     delete f0;
00719     return sqrt_reduce (sqr_reduce (f1));
00720   } else if (isZero (f1)) {
00721     delete f1;
00722     return sqrt_reduce (sqr_reduce (f0));
00723   } else if (isConst (f0) && isConst (f1)) {
00724     nr_double_t t = ::xhypot (D(f0), D(f1));
00725     delete f0; delete f1;
00726     retCon (t);
00727   } else {
00728     retApp2 ("hypot", f0, f1);
00729   }
00730 }
00731 
00732 #include "constants.h"
00733 
00734 node * differentiate::vt (application * app, char * derivative) {
00735   _AD0 (d0);
00736   defCon (con, kBoverQ);
00737   return times_reduce (d0, con);
00738 }
00739 
00740 // List of differentiators.
00741 struct differentiation_t eqn::differentiations[] = {
00742   { "+", differentiate::plus_binary,  2 },
00743   { "+", differentiate::plus_unary,   1 },
00744   { "-", differentiate::minus_binary, 2 },
00745   { "-", differentiate::minus_unary,  1 },
00746   { "*", differentiate::times,        2 },
00747   { "/", differentiate::over,         2 },
00748   { "^", differentiate::power,        2 },
00749 
00750   { "?:", differentiate::ifthenelse, 3 },
00751 
00752   { "ln",       differentiate::ln,        1 },
00753   { "log10",    differentiate::log10,     1 },
00754   { "log2",     differentiate::log2,      1 },
00755   { "sqrt",     differentiate::sqrt,      1 },
00756   { "exp",      differentiate::exp,       1 },
00757   { "sinc",     differentiate::sinc,      1 },
00758   { "norm",     differentiate::norm,      1 },
00759   { "sin",      differentiate::sin,       1 },
00760   { "cos",      differentiate::cos,       1 },
00761   { "tan",      differentiate::tan,       1 },
00762   { "sec",      differentiate::sec,       1 },
00763   { "cot",      differentiate::cot,       1 },
00764   { "cosec",    differentiate::cosec,     1 },
00765   { "abs",      differentiate::abs,       1 },
00766   { "step",     differentiate::step,      1 },
00767   { "sign",     differentiate::sign,      1 },
00768   { "arcsin",   differentiate::arcsin,    1 },
00769   { "arccos",   differentiate::arccos,    1 },
00770   { "arctan",   differentiate::arctan,    1 },
00771   { "arccot",   differentiate::arccot,    1 },
00772   { "arcsec",   differentiate::arcsec,    1 },
00773   { "arccosec", differentiate::arccosec,  1 },
00774   { "sqr",      differentiate::square,    1 },
00775   { "sinh",     differentiate::sinh,      1 },
00776   { "cosh",     differentiate::cosh,      1 },
00777   { "tanh",     differentiate::tanh,      1 },
00778   { "coth",     differentiate::coth,      1 },
00779   { "arsinh",   differentiate::arsinh,    1 },
00780   { "arcosh",   differentiate::arcosh,    1 },
00781   { "artanh",   differentiate::artanh,    1 },
00782   { "arcoth",   differentiate::arcoth,    1 },
00783   { "arsech",   differentiate::arsech,    1 },
00784   { "arcosech", differentiate::arcosech,  1 },
00785   { "hypot",    differentiate::xhypot,    2 },
00786   { "limexp",   differentiate::limexp,    1 },
00787   { "vt",       differentiate::vt,        1 },
00788 
00789   { NULL, NULL, 0 /* end of list */ }
00790 };