123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011 |
- /* derivative */
- #include "stdafx.h"
- #include "defs.h"
- #define F p3
- #define X p4
- #define N p5
- void
- eval_derivative(void)
- {
- int i, n;
- // evaluate 1st arg to get function F
- p1 = cdr(p1);
- push(car(p1));
- eval();
- // evaluate 2nd arg and then...
- // example result of 2nd arg what to do
- //
- // d(f) nil guess X, N = nil
- // d(f,2) 2 guess X, N = 2
- // d(f,x) x X = x, N = nil
- // d(f,x,2) x X = x, N = 2
- // d(f,x,y) x X = x, N = y
- p1 = cdr(p1);
- push(car(p1));
- eval();
- p2 = pop();
- if (p2 == symbol(NIL)) {
- guess();
- push(symbol(NIL));
- } else if (isnum(p2)) {
- guess();
- push(p2);
- } else {
- push(p2);
- p1 = cdr(p1);
- push(car(p1));
- eval();
- }
- N = pop();
- X = pop();
- F = pop();
- while (1) {
- // N might be a symbol instead of a number
- if (isnum(N)) {
- push(N);
- n = pop_integer();
- if (n == (int) 0x80000000)
- stop("nth derivative: check n");
- } else
- n = 1;
- push(F);
- if (n >= 0) {
- for (i = 0; i < n; i++) {
- push(X);
- derivative();
- }
- } else {
- n = -n;
- for (i = 0; i < n; i++) {
- push(X);
- integral();
- }
- }
- F = pop();
- // if N is nil then arglist is exhausted
- if (N == symbol(NIL))
- break;
- // otherwise...
- // N arg1 what to do
- //
- // number nil break
- // number number N = arg1, continue
- // number symbol X = arg1, N = arg2, continue
- //
- // symbol nil X = N, N = nil, continue
- // symbol number X = N, N = arg1, continue
- // symbol symbol X = N, N = arg1, continue
- if (isnum(N)) {
- p1 = cdr(p1);
- push(car(p1));
- eval();
- N = pop();
- if (N == symbol(NIL))
- break; // arglist exhausted
- if (isnum(N))
- ; // N = arg1
- else {
- X = N; // X = arg1
- p1 = cdr(p1);
- push(car(p1));
- eval();
- N = pop(); // N = arg2
- }
- } else {
- X = N; // X = N
- p1 = cdr(p1);
- push(car(p1));
- eval();
- N = pop(); // N = arg1
- }
- }
- push(F); // final result
- }
- void
- derivative(void)
- {
- save();
- p2 = pop();
- p1 = pop();
- if (isnum(p2))
- stop("undefined function");
- if (istensor(p1))
- if (istensor(p2))
- d_tensor_tensor();
- else
- d_tensor_scalar();
- else
- if (istensor(p2))
- d_scalar_tensor();
- else
- d_scalar_scalar();
- restore();
- }
- void
- d_scalar_scalar(void)
- {
- if (issymbol(p2))
- d_scalar_scalar_1();
- else {
- // Example: d(sin(cos(x)),cos(x))
- // Replace cos(x) <- X, find derivative, then do X <- cos(x)
- push(p1); // sin(cos(x))
- push(p2); // cos(x)
- push(symbol(SECRETX)); // X
- subst(); // sin(cos(x)) -> sin(X)
- push(symbol(SECRETX)); // X
- derivative();
- push(symbol(SECRETX)); // X
- push(p2); // cos(x)
- subst(); // cos(X) -> cos(cos(x))
- }
- }
- void
- d_scalar_scalar_1(void)
- {
- // d(x,x)?
- if (equal(p1, p2)) {
- push(one);
- return;
- }
- // d(a,x)?
- if (!iscons(p1)) {
- push(zero);
- return;
- }
- if (isadd(p1)) {
- dsum();
- return;
- }
- if (car(p1) == symbol(MULTIPLY)) {
- dproduct();
- return;
- }
- if (car(p1) == symbol(POWER)) {
- dpower();
- return;
- }
- if (car(p1) == symbol(DERIVATIVE)) {
- dd();
- return;
- }
- if (car(p1) == symbol(LOG)) {
- dlog();
- return;
- }
- if (car(p1) == symbol(SIN)) {
- dsin();
- return;
- }
- if (car(p1) == symbol(COS)) {
- dcos();
- return;
- }
- if (car(p1) == symbol(TAN)) {
- dtan();
- return;
- }
- if (car(p1) == symbol(ARCSIN)) {
- darcsin();
- return;
- }
- if (car(p1) == symbol(ARCCOS)) {
- darccos();
- return;
- }
- if (car(p1) == symbol(ARCTAN)) {
- darctan();
- return;
- }
- if (car(p1) == symbol(SINH)) {
- dsinh();
- return;
- }
- if (car(p1) == symbol(COSH)) {
- dcosh();
- return;
- }
- if (car(p1) == symbol(TANH)) {
- dtanh();
- return;
- }
- if (car(p1) == symbol(ARCSINH)) {
- darcsinh();
- return;
- }
- if (car(p1) == symbol(ARCCOSH)) {
- darccosh();
- return;
- }
- if (car(p1) == symbol(ARCTANH)) {
- darctanh();
- return;
- }
- if (car(p1) == symbol(ABS)) {
- dabs();
- return;
- }
- if (car(p1) == symbol(SGN)) {
- dsgn();
- return;
- }
- if (car(p1) == symbol(HERMITE)) {
- dhermite();
- return;
- }
- if (car(p1) == symbol(ERF)) {
- derf();
- return;
- }
- if (car(p1) == symbol(ERFC)) {
- derfc();
- return;
- }
- /*
- if (car(p1) == symbol(BESSELJ)) {
- if (iszero(caddr(p1)))
- dbesselj0();
- else
- dbesseljn();
- return;
- }
- if (car(p1) == symbol(BESSELY)) {
- if (iszero(caddr(p1)))
- dbessely0();
- else
- dbesselyn();
- return;
- }
- */
- if (car(p1) == symbol(INTEGRAL) && caddr(p1) == p2) {
- derivative_of_integral();
- return;
- }
- dfunction();
- }
- void
- dsum(void)
- {
- int h = tos;
- p1 = cdr(p1);
- while (iscons(p1)) {
- push(car(p1));
- push(p2);
- derivative();
- p1 = cdr(p1);
- }
- add_all(tos - h);
- }
- void
- dproduct(void)
- {
- int i, j, n;
- n = length(p1) - 1;
- for (i = 0; i < n; i++) {
- p3 = cdr(p1);
- for (j = 0; j < n; j++) {
- push(car(p3));
- if (i == j) {
- push(p2);
- derivative();
- }
- p3 = cdr(p3);
- }
- multiply_all(n);
- }
- add_all(n);
- }
- //-----------------------------------------------------------------------------
- //
- // v
- // y = u
- //
- // log y = v log u
- //
- // 1 dy v du dv
- // - -- = - -- + (log u) --
- // y dx u dx dx
- //
- // dy v v du dv
- // -- = u (- -- + (log u) --)
- // dx u dx dx
- //
- //-----------------------------------------------------------------------------
- void
- dpower(void)
- {
- push(caddr(p1)); // v/u
- push(cadr(p1));
- divide();
- push(cadr(p1)); // du/dx
- push(p2);
- derivative();
- multiply();
- push(cadr(p1)); // log u
- logarithm();
- push(caddr(p1)); // dv/dx
- push(p2);
- derivative();
- multiply();
- add();
- push(p1); // u^v
- multiply();
- }
- void
- dlog(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- divide();
- }
- // derivative of derivative
- //
- // example: d(d(f(x,y),y),x)
- //
- // p1 = d(f(x,y),y)
- //
- // p2 = x
- //
- // cadr(p1) = f(x,y)
- //
- // caddr(p1) = y
- void
- dd(void)
- {
- // d(f(x,y),x)
- push(cadr(p1));
- push(p2);
- derivative();
- p3 = pop();
- if (car(p3) == symbol(DERIVATIVE)) {
- // sort dx terms
- push_symbol(DERIVATIVE);
- push_symbol(DERIVATIVE);
- push(cadr(p3));
- if (lessp(caddr(p3), caddr(p1))) {
- push(caddr(p3));
- list(3);
- push(caddr(p1));
- } else {
- push(caddr(p1));
- list(3);
- push(caddr(p3));
- }
- list(3);
- } else {
- push(p3);
- push(caddr(p1));
- derivative();
- }
- }
- // derivative of a generic function
- void
- dfunction(void)
- {
- p3 = cdr(p1); // p3 is the argument list for the function
- if (p3 == symbol(NIL) || find(p3, p2)) {
- push_symbol(DERIVATIVE);
- push(p1);
- push(p2);
- list(3);
- } else
- push(zero);
- }
- void
- dsin(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- cosine();
- multiply();
- }
- void
- dcos(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- sine();
- multiply();
- negate();
- }
- void
- dtan(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- cosine();
- push_integer(-2);
- power();
- multiply();
- }
- void
- darcsin(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push_integer(1);
- push(cadr(p1));
- push_integer(2);
- power();
- subtract();
- push_rational(-1, 2);
- power();
- multiply();
- }
- void
- darccos(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push_integer(1);
- push(cadr(p1));
- push_integer(2);
- power();
- subtract();
- push_rational(-1, 2);
- power();
- multiply();
- negate();
- }
- // Without simplify With simplify
- //
- // d(arctan(y/x),x) -y/(x^2*(y^2/x^2+1)) -y/(x^2+y^2)
- //
- // d(arctan(y/x),y) 1/(x*(y^2/x^2+1)) x/(x^2+y^2)
- void
- darctan(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push_integer(1);
- push(cadr(p1));
- push_integer(2);
- power();
- add();
- inverse();
- multiply();
- simplify();
- }
- void
- dsinh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- ycosh();
- multiply();
- }
- void
- dcosh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- ysinh();
- multiply();
- }
- void
- dtanh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- ycosh();
- push_integer(-2);
- power();
- multiply();
- }
- void
- darcsinh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push_integer(2);
- power();
- push_integer(1);
- add();
- push_rational(-1, 2);
- power();
- multiply();
- }
- void
- darccosh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push_integer(2);
- power();
- push_integer(-1);
- add();
- push_rational(-1, 2);
- power();
- multiply();
- }
- void
- darctanh(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push_integer(1);
- push(cadr(p1));
- push_integer(2);
- power();
- subtract();
- inverse();
- multiply();
- }
- void
- dabs(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- sgn();
- multiply();
- }
- void
- dsgn(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- dirac();
- multiply();
- push_integer(2);
- multiply();
- }
- void
- dhermite(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push_integer(2);
- push(caddr(p1));
- multiply();
- multiply();
- push(cadr(p1));
- push(caddr(p1));
- push_integer(-1);
- add();
- hermite();
- multiply();
- }
- void
- derf(void)
- {
- push(cadr(p1));
- push_integer(2);
- power();
- push_integer(-1);
- multiply();
- exponential();
- push_symbol(PI);
- push_rational(-1,2);
- power();
- multiply();
- push_integer(2);
- multiply();
- push(cadr(p1));
- push(p2);
- derivative();
- multiply();
- }
- void
- derfc(void)
- {
- push(cadr(p1));
- push_integer(2);
- power();
- push_integer(-1);
- multiply();
- exponential();
- push_symbol(PI);
- push_rational(-1,2);
- power();
- multiply();
- push_integer(-2);
- multiply();
- push(cadr(p1));
- push(p2);
- derivative();
- multiply();
- }
- /*
- void
- dbesselj0(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push_integer(1);
- besselj();
- multiply();
- push_integer(-1);
- multiply();
- }
- void
- dbesseljn(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push(caddr(p1));
- push_integer(-1);
- add();
- besselj();
- push(caddr(p1));
- push_integer(-1);
- multiply();
- push(cadr(p1));
- divide();
- push(cadr(p1));
- push(caddr(p1));
- besselj();
- multiply();
- add();
- multiply();
- }
- void
- dbessely0(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push_integer(1);
- besselj();
- multiply();
- push_integer(-1);
- multiply();
- }
- void
- dbesselyn(void)
- {
- push(cadr(p1));
- push(p2);
- derivative();
- push(cadr(p1));
- push(caddr(p1));
- push_integer(-1);
- add();
- bessely();
- push(caddr(p1));
- push_integer(-1);
- multiply();
- push(cadr(p1));
- divide();
- push(cadr(p1));
- push(caddr(p1));
- bessely();
- multiply();
- add();
- multiply();
- }
- */
- void
- derivative_of_integral(void)
- {
- push(cadr(p1));
- }
- #if SELFTEST
- static char *s[] = {
- "x=quote(x)",
- "",
- "f=quote(f)",
- "",
- "g=quote(g)",
- "",
- "d(a,x)",
- "0",
- "d(x,x)",
- "1",
- "d(x^2,x)",
- "2*x",
- "d(log(x),x)",
- "1/x",
- "d(exp(x),x)",
- "exp(x)",
- "d(a^x,x)",
- "a^x*log(a)",
- "d(x^x,x)-(x^x+x^x*log(x))",
- "0",
- "d(log(x^2+5),x)-(2*x/(5+x^2))",
- "0",
- "d(d(f(x),x),y)",
- "0",
- "d(d(f(x),y),x)",
- "0",
- "d(d(f(y),x),y)",
- "0",
- "d(d(f(y),y),x)",
- "0",
- "d((x*y*z,y,x+z),(x,y,z))",
- "((y*z,x*z,x*y),(0,1,0),(1,0,1))",
- "d(x+z,(x,y,z))",
- "(1,0,1)",
- "d(cos(theta)^2,cos(theta))",
- "2*cos(theta)",
- "d(f())",
- "d(f(),x)",
- "d(x^2)",
- "2*x",
- "d(t^2)",
- "2*t",
- "d(t^2 x^2)",
- "2*t^2*x",
- // trig functions
- "d(sin(x),x)-cos(x)",
- "0",
- "d(cos(x),x)+sin(x)",
- "0",
- "d(tan(x),x)-cos(x)^(-2)",
- "0",
- "d(arcsin(x),x)-1/sqrt(1-x^2)",
- "0",
- "d(arccos(x),x)+1/sqrt(1-x^2)",
- "0",
- "d(arctan(x),x)-1/(1+x^2)",
- "0",
- "d(arctan(y/x),x)",
- "-y/(x^2+y^2)",
- "d(arctan(y/x),y)",
- "x/(x^2+y^2)",
- // hyp functions
- "d(sinh(x),x)-cosh(x)",
- "0",
- "d(cosh(x),x)-sinh(x)",
- "0",
- "d(tanh(x),x)-cosh(x)^(-2)",
- "0",
- "d(arcsinh(x),x)-1/sqrt(x^2+1)",
- "0",
- "d(arccosh(x),x)-1/sqrt(x^2-1)",
- "0",
- "d(arctanh(x),x)-1/(1-x^2)",
- "0",
- "d(sin(cos(x)),x)+cos(cos(x))*sin(x)",
- "0",
- "d(sin(x)^2,x)-2*sin(x)*cos(x)",
- "0",
- "d(sin(cos(x)),cos(x))-cos(cos(x))",
- "0",
- "d(abs(x),x)",
- "sgn(x)",
-
- "d(sgn(x),x)",
- "2*dirac(x)",
- // generic functions
- "d(f(),x)",
- "d(f(),x)",
- "d(f(x),x)",
- "d(f(x),x)",
- "d(f(y),x)",
- "0",
- "d(g(f(x)),f(x))",
- "d(g(f(x)),f(x))",
- "d(g(f(x)),x)",
- "d(g(f(x)),x)",
- // other functions
- "d(erf(x))-2*exp(-x^2)/sqrt(pi)",
- "0",
- // arg lists
- "f=x^5*y^7",
- "",
- "d(f)",
- "5*x^4*y^7",
- "d(f,x)",
- "5*x^4*y^7",
- "d(f,x,0)",
- "x^5*y^7",
- "d(f,x,1)",
- "5*x^4*y^7",
- "d(f,x,2)",
- "20*x^3*y^7",
- "d(f,2)",
- "20*x^3*y^7",
- "d(f,2,y)",
- "140*x^3*y^6",
- "d(f,x,x,y,y)",
- "840*x^3*y^5",
- "f=quote(f)",
- "",
- };
- void
- test_derivative(void)
- {
- test(__FILE__, s, sizeof s / sizeof (char *));
- }
- #endif
|