derivative.cpp

   1 /* derivative */
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 #define F p3
   7 #define X p4
   8 #define N p5
   9
  10 void
  11 eval_derivative(void)
  12 {
  13         int i, n;
  14
  15         // evaluate 1st arg to get function F
  16
  17         p1 = cdr(p1);
  18         push(car(p1));
  19         eval();
  20
  21         // evaluate 2nd arg and then...
  22
  23         // example      result of 2nd arg       what to do
  24         //
  25         // d(f)         nil                     guess X, N = nil
  26         // d(f,2)       2                       guess X, N = 2
  27         // d(f,x)       x                       X = x, N = nil
  28         // d(f,x,2)     x                       X = x, N = 2
  29         // d(f,x,y)     x                       X = x, N = y
  30
  31         p1 = cdr(p1);
  32         push(car(p1));
  33         eval();
  34
  35         p2 = pop();
  36         if (p2 == symbol(NIL)) {
  37                 guess();
  38                 push(symbol(NIL));
  39         } else if (isnum(p2)) {
  40                 guess();
  41                 push(p2);
  42         } else {
  43                 push(p2);
  44                 p1 = cdr(p1);
  45                 push(car(p1));
  46                 eval();
  47         }
  48
  49         N = pop();
  50         X = pop();
  51         F = pop();
  52
  53         while (1) {
  54
  55                 // N might be a symbol instead of a number
  56
  57                 if (isnum(N)) {
  58                         push(N);
  59                         n = pop_integer();
  60                         if (n == (int) 0x80000000)
  61                                 stop("nth derivative: check n");
  62                 } else
  63                         n = 1;
  64
  65                 push(F);
  66
  67                 if (n >= 0) {
  68                         for (i = 0; i < n; i++) {
  69                                 push(X);
  70                                 derivative();
  71                         }
  72                 } else {
  73                         n = -n;
  74                         for (i = 0; i < n; i++) {
  75                                 push(X);
  76                                 integral();
  77                         }
  78                 }
  79
  80                 F = pop();
  81
  82                 // if N is nil then arglist is exhausted
  83
  84                 if (N == symbol(NIL))
  85                         break;
  86
  87                 // otherwise...
  88
  89                 // N            arg1            what to do
  90                 //
  91                 // number       nil             break
  92                 // number       number          N = arg1, continue
  93                 // number       symbol          X = arg1, N = arg2, continue
  94                 //
  95                 // symbol       nil             X = N, N = nil, continue
  96                 // symbol       number          X = N, N = arg1, continue
  97                 // symbol       symbol          X = N, N = arg1, continue
  98
  99                 if (isnum(N)) {
 100                         p1 = cdr(p1);
 101                         push(car(p1));
 102                         eval();
 103                         N = pop();
 104                         if (N == symbol(NIL))
 105                                 break;          // arglist exhausted
 106                         if (isnum(N))
 107                                 ;               // N = arg1
 108                         else {
 109                                 X = N;          // X = arg1
 110                                 p1 = cdr(p1);
 111                                 push(car(p1));
 112                                 eval();
 113                                 N = pop();      // N = arg2
 114                         }
 115                 } else {
 116                         X = N;                  // X = N
 117                         p1 = cdr(p1);
 118                         push(car(p1));
 119                         eval();
 120                         N = pop();              // N = arg1
 121                 }
 122         }
 123
 124         push(F); // final result
 125 }
 126
 127 void
 128 derivative(void)
 129 {
 130         save();
 131         p2 = pop();
 132         p1 = pop();
 133         if (isnum(p2))
 134                 stop("undefined function");
 135         if (istensor(p1))
 136                 if (istensor(p2))
 137                         d_tensor_tensor();
 138                 else
 139                         d_tensor_scalar();
 140         else
 141                 if (istensor(p2))
 142                         d_scalar_tensor();
 143                 else
 144                         d_scalar_scalar();
 145         restore();
 146 }
 147
 148 void
 149 d_scalar_scalar(void)
 150 {
 151         if (issymbol(p2))
 152                 d_scalar_scalar_1();
 153         else {
 154                 // Example: d(sin(cos(x)),cos(x))
 155                 // Replace cos(x) <- X, find derivative, then do X <- cos(x)
 156                 push(p1);               // sin(cos(x))
 157                 push(p2);               // cos(x)
 158                 push(symbol(SECRETX));  // X
 159                 subst();                // sin(cos(x)) -> sin(X)
 160                 push(symbol(SECRETX));  // X
 161                 derivative();
 162                 push(symbol(SECRETX));  // X
 163                 push(p2);               // cos(x)
 164                 subst();                // cos(X) -> cos(cos(x))
 165         }
 166 }
 167
 168 void
 169 d_scalar_scalar_1(void)
 170 {
 171         // d(x,x)?
 172
 173         if (equal(p1, p2)) {
 174                 push(one);
 175                 return;
 176         }
 177
 178         // d(a,x)?
 179
 180         if (!iscons(p1)) {
 181                 push(zero);
 182                 return;
 183         }
 184
 185         if (isadd(p1)) {
 186                 dsum();
 187                 return;
 188         }
 189
 190         if (car(p1) == symbol(MULTIPLY)) {
 191                 dproduct();
 192                 return;
 193         }
 194
 195         if (car(p1) == symbol(POWER)) {
 196                 dpower();
 197                 return;
 198         }
 199
 200         if (car(p1) == symbol(DERIVATIVE)) {
 201                 dd();
 202                 return;
 203         }
 204
 205         if (car(p1) == symbol(LOG)) {
 206                 dlog();
 207                 return;
 208         }
 209
 210         if (car(p1) == symbol(SIN)) {
 211                 dsin();
 212                 return;
 213         }
 214
 215         if (car(p1) == symbol(COS)) {
 216                 dcos();
 217                 return;
 218         }
 219
 220         if (car(p1) == symbol(TAN)) {
 221                 dtan();
 222                 return;
 223         }
 224
 225         if (car(p1) == symbol(ARCSIN)) {
 226                 darcsin();
 227                 return;
 228         }
 229
 230         if (car(p1) == symbol(ARCCOS)) {
 231                 darccos();
 232                 return;
 233         }
 234
 235         if (car(p1) == symbol(ARCTAN)) {
 236                 darctan();
 237                 return;
 238         }
 239
 240         if (car(p1) == symbol(SINH)) {
 241                 dsinh();
 242                 return;
 243         }
 244
 245         if (car(p1) == symbol(COSH)) {
 246                 dcosh();
 247                 return;
 248         }
 249
 250         if (car(p1) == symbol(TANH)) {
 251                 dtanh();
 252                 return;
 253         }
 254
 255         if (car(p1) == symbol(ARCSINH)) {
 256                 darcsinh();
 257                 return;
 258         }
 259
 260         if (car(p1) == symbol(ARCCOSH)) {
 261                 darccosh();
 262                 return;
 263         }
 264
 265         if (car(p1) == symbol(ARCTANH)) {
 266                 darctanh();
 267                 return;
 268         }
 269
 270         if (car(p1) == symbol(ABS)) {
 271                 dabs();
 272                 return;
 273         }
 274
 275         if (car(p1) == symbol(SGN)) {
 276                 dsgn();
 277                 return;
 278         }
 279
 280         if (car(p1) == symbol(HERMITE)) {
 281                 dhermite();
 282                 return;
 283         }
 284
 285         if (car(p1) == symbol(ERF)) {
 286                 derf();
 287                 return;
 288         }
 289
 290         if (car(p1) == symbol(ERFC)) {
 291                 derfc();
 292                 return;
 293         }
 294 /*
 295         if (car(p1) == symbol(BESSELJ)) {
 296                 if (iszero(caddr(p1)))
 297                         dbesselj0();
 298                 else
 299                         dbesseljn();
 300                 return;
 301         }
 302
 303         if (car(p1) == symbol(BESSELY)) {
 304                 if (iszero(caddr(p1)))
 305                         dbessely0();
 306                 else
 307                         dbesselyn();
 308                 return;
 309         }
 310 */
 311         if (car(p1) == symbol(INTEGRAL) && caddr(p1) == p2) {
 312                 derivative_of_integral();
 313                 return;
 314         }
 315
 316         dfunction();
 317 }
 318
 319 void
 320 dsum(void)
 321 {
 322         int h = tos;
 323         p1 = cdr(p1);
 324         while (iscons(p1)) {
 325                 push(car(p1));
 326                 push(p2);
 327                 derivative();
 328                 p1 = cdr(p1);
 329         }
 330         add_all(tos - h);
 331 }
 332
 333 void
 334 dproduct(void)
 335 {
 336         int i, j, n;
 337         n = length(p1) - 1;
 338         for (i = 0; i < n; i++) {
 339                 p3 = cdr(p1);
 340                 for (j = 0; j < n; j++) {
 341                         push(car(p3));
 342                         if (i == j) {
 343                                 push(p2);
 344                                 derivative();
 345                         }
 346                         p3 = cdr(p3);
 347                 }
 348                 multiply_all(n);
 349         }
 350         add_all(n);
 351 }
 352
 353 //-----------------------------------------------------------------------------
 354 //
 355 //           v
 356 //      y = u
 357 //
 358 //      log y = v log u
 359 //
 360 //      1 dy   v du           dv
 361 //      - -- = - -- + (log u) --
 362 //      y dx   u dx           dx
 363 //
 364 //      dy    v  v du           dv
 365 //      -- = u  (- -- + (log u) --)
 366 //      dx       u dx           dx
 367 //
 368 //-----------------------------------------------------------------------------
 369
 370 void
 371 dpower(void)
 372 {
 373         push(caddr(p1));        // v/u
 374         push(cadr(p1));
 375         divide();
 376
 377         push(cadr(p1));         // du/dx
 378         push(p2);
 379         derivative();
 380
 381         multiply();
 382
 383         push(cadr(p1));         // log u
 384         logarithm();
 385
 386         push(caddr(p1));        // dv/dx
 387         push(p2);
 388         derivative();
 389
 390         multiply();
 391
 392         add();
 393
 394         push(p1);               // u^v
 395
 396         multiply();
 397 }
 398
 399 void
 400 dlog(void)
 401 {
 402         push(cadr(p1));
 403         push(p2);
 404         derivative();
 405         push(cadr(p1));
 406         divide();
 407 }
 408
 409 //      derivative of derivative
 410 //
 411 //      example: d(d(f(x,y),y),x)
 412 //
 413 //      p1 = d(f(x,y),y)
 414 //
 415 //      p2 = x
 416 //
 417 //      cadr(p1) = f(x,y)
 418 //
 419 //      caddr(p1) = y
 420
 421 void
 422 dd(void)
 423 {
 424         // d(f(x,y),x)
 425
 426         push(cadr(p1));
 427         push(p2);
 428         derivative();
 429
 430         p3 = pop();
 431
 432         if (car(p3) == symbol(DERIVATIVE)) {
 433
 434                 // sort dx terms
 435
 436                 push_symbol(DERIVATIVE);
 437                 push_symbol(DERIVATIVE);
 438                 push(cadr(p3));
 439
 440                 if (lessp(caddr(p3), caddr(p1))) {
 441                         push(caddr(p3));
 442                         list(3);
 443                         push(caddr(p1));
 444                 } else {
 445                         push(caddr(p1));
 446                         list(3);
 447                         push(caddr(p3));
 448                 }
 449
 450                 list(3);
 451
 452         } else {
 453                 push(p3);
 454                 push(caddr(p1));
 455                 derivative();
 456         }
 457 }
 458
 459 // derivative of a generic function
 460
 461 void
 462 dfunction(void)
 463 {
 464         p3 = cdr(p1);   // p3 is the argument list for the function
 465
 466         if (p3 == symbol(NIL) || find(p3, p2)) {
 467                 push_symbol(DERIVATIVE);
 468                 push(p1);
 469                 push(p2);
 470                 list(3);
 471         } else
 472                 push(zero);
 473 }
 474
 475 void
 476 dsin(void)
 477 {
 478         push(cadr(p1));
 479         push(p2);
 480         derivative();
 481         push(cadr(p1));
 482         cosine();
 483         multiply();
 484 }
 485
 486 void
 487 dcos(void)
 488 {
 489         push(cadr(p1));
 490         push(p2);
 491         derivative();
 492         push(cadr(p1));
 493         sine();
 494         multiply();
 495         negate();
 496 }
 497
 498 void
 499 dtan(void)
 500 {
 501         push(cadr(p1));
 502         push(p2);
 503         derivative();
 504         push(cadr(p1));
 505         cosine();
 506         push_integer(-2);
 507         power();
 508         multiply();
 509 }
 510
 511 void
 512 darcsin(void)
 513 {
 514         push(cadr(p1));
 515         push(p2);
 516         derivative();
 517         push_integer(1);
 518         push(cadr(p1));
 519         push_integer(2);
 520         power();
 521         subtract();
 522         push_rational(-1, 2);
 523         power();
 524         multiply();
 525 }
 526
 527 void
 528 darccos(void)
 529 {
 530         push(cadr(p1));
 531         push(p2);
 532         derivative();
 533         push_integer(1);
 534         push(cadr(p1));
 535         push_integer(2);
 536         power();
 537         subtract();
 538         push_rational(-1, 2);
 539         power();
 540         multiply();
 541         negate();
 542 }
 543
 544 //                              Without simplify        With simplify
 545 //
 546 //      d(arctan(y/x),x)        -y/(x^2*(y^2/x^2+1))    -y/(x^2+y^2)
 547 //
 548 //      d(arctan(y/x),y)        1/(x*(y^2/x^2+1))       x/(x^2+y^2)
 549
 550 void
 551 darctan(void)
 552 {
 553         push(cadr(p1));
 554         push(p2);
 555         derivative();
 556         push_integer(1);
 557         push(cadr(p1));
 558         push_integer(2);
 559         power();
 560         add();
 561         inverse();
 562         multiply();
 563         simplify();
 564 }
 565
 566 void
 567 dsinh(void)
 568 {
 569         push(cadr(p1));
 570         push(p2);
 571         derivative();
 572         push(cadr(p1));
 573         ycosh();
 574         multiply();
 575 }
 576
 577 void
 578 dcosh(void)
 579 {
 580         push(cadr(p1));
 581         push(p2);
 582         derivative();
 583         push(cadr(p1));
 584         ysinh();
 585         multiply();
 586 }
 587
 588 void
 589 dtanh(void)
 590 {
 591         push(cadr(p1));
 592         push(p2);
 593         derivative();
 594         push(cadr(p1));
 595         ycosh();
 596         push_integer(-2);
 597         power();
 598         multiply();
 599 }
 600
 601 void
 602 darcsinh(void)
 603 {
 604         push(cadr(p1));
 605         push(p2);
 606         derivative();
 607         push(cadr(p1));
 608         push_integer(2);
 609         power();
 610         push_integer(1);
 611         add();
 612         push_rational(-1, 2);
 613         power();
 614         multiply();
 615 }
 616
 617 void
 618 darccosh(void)
 619 {
 620         push(cadr(p1));
 621         push(p2);
 622         derivative();
 623         push(cadr(p1));
 624         push_integer(2);
 625         power();
 626         push_integer(-1);
 627         add();
 628         push_rational(-1, 2);
 629         power();
 630         multiply();
 631 }
 632
 633 void
 634 darctanh(void)
 635 {
 636         push(cadr(p1));
 637         push(p2);
 638         derivative();
 639         push_integer(1);
 640         push(cadr(p1));
 641         push_integer(2);
 642         power();
 643         subtract();
 644         inverse();
 645         multiply();
 646 }
 647
 648 void
 649 dabs(void)
 650 {
 651         push(cadr(p1));
 652         push(p2);
 653         derivative();
 654         push(cadr(p1));
 655         sgn();
 656         multiply();
 657 }
 658
 659 void
 660 dsgn(void)
 661 {
 662         push(cadr(p1));
 663         push(p2);
 664         derivative();
 665         push(cadr(p1));
 666         dirac();
 667         multiply();
 668         push_integer(2);
 669         multiply();
 670 }
 671
 672 void
 673 dhermite(void)
 674 {
 675         push(cadr(p1));
 676         push(p2);
 677         derivative();
 678         push_integer(2);
 679         push(caddr(p1));
 680         multiply();
 681         multiply();
 682         push(cadr(p1));
 683         push(caddr(p1));
 684         push_integer(-1);
 685         add();
 686         hermite();
 687         multiply();
 688 }
 689
 690 void
 691 derf(void)
 692 {
 693         push(cadr(p1));
 694         push_integer(2);
 695         power();
 696         push_integer(-1);
 697         multiply();
 698         exponential();
 699         push_symbol(PI);
 700         push_rational(-1,2);
 701         power();
 702         multiply();
 703         push_integer(2);
 704         multiply();
 705         push(cadr(p1));
 706         push(p2);
 707         derivative();
 708         multiply();
 709
 710 }
 711
 712 void
 713 derfc(void)
 714 {
 715         push(cadr(p1));
 716         push_integer(2);
 717         power();
 718         push_integer(-1);
 719         multiply();
 720         exponential();
 721         push_symbol(PI);
 722         push_rational(-1,2);
 723         power();
 724         multiply();
 725         push_integer(-2);
 726         multiply();
 727         push(cadr(p1));
 728         push(p2);
 729         derivative();
 730         multiply();
 731
 732 }
 733 /*
 734 void
 735 dbesselj0(void)
 736 {
 737         push(cadr(p1));
 738         push(p2);
 739         derivative();
 740         push(cadr(p1));
 741         push_integer(1);
 742         besselj();
 743         multiply();
 744         push_integer(-1);
 745         multiply();
 746 }
 747
 748 void
 749 dbesseljn(void)
 750 {
 751         push(cadr(p1));
 752         push(p2);
 753         derivative();
 754         push(cadr(p1));
 755         push(caddr(p1));
 756         push_integer(-1);
 757         add();
 758         besselj();
 759         push(caddr(p1));
 760         push_integer(-1);
 761         multiply();
 762         push(cadr(p1));
 763         divide();
 764         push(cadr(p1));
 765         push(caddr(p1));
 766         besselj();
 767         multiply();
 768         add();
 769         multiply();
 770 }
 771
 772
 773 void
 774 dbessely0(void)
 775 {
 776         push(cadr(p1));
 777         push(p2);
 778         derivative();
 779         push(cadr(p1));
 780         push_integer(1);
 781         besselj();
 782         multiply();
 783         push_integer(-1);
 784         multiply();
 785 }
 786
 787
 788 void
 789 dbesselyn(void)
 790 {
 791         push(cadr(p1));
 792         push(p2);
 793         derivative();
 794         push(cadr(p1));
 795         push(caddr(p1));
 796         push_integer(-1);
 797         add();
 798         bessely();
 799         push(caddr(p1));
 800         push_integer(-1);
 801         multiply();
 802         push(cadr(p1));
 803         divide();
 804         push(cadr(p1));
 805         push(caddr(p1));
 806         bessely();
 807         multiply();
 808         add();
 809         multiply();
 810 }
 811 */
 812 void
 813 derivative_of_integral(void)
 814 {
 815         push(cadr(p1));
 816 }
 817
 818 #if SELFTEST
 819
 820 static char *s[] = {
 821
 822         "x=quote(x)",
 823         "",
 824
 825         "f=quote(f)",
 826         "",
 827
 828         "g=quote(g)",
 829         "",
 830
 831         "d(a,x)",
 832         "0",
 833
 834         "d(x,x)",
 835         "1",
 836
 837         "d(x^2,x)",
 838         "2*x",
 839
 840         "d(log(x),x)",
 841         "1/x",
 842
 843         "d(exp(x),x)",
 844         "exp(x)",
 845
 846         "d(a^x,x)",
 847         "a^x*log(a)",
 848
 849         "d(x^x,x)-(x^x+x^x*log(x))",
 850         "0",
 851
 852         "d(log(x^2+5),x)-(2*x/(5+x^2))",
 853         "0",
 854
 855         "d(d(f(x),x),y)",
 856         "0",
 857
 858         "d(d(f(x),y),x)",
 859         "0",
 860
 861         "d(d(f(y),x),y)",
 862         "0",
 863
 864         "d(d(f(y),y),x)",
 865         "0",
 866
 867         "d((x*y*z,y,x+z),(x,y,z))",
 868         "((y*z,x*z,x*y),(0,1,0),(1,0,1))",
 869
 870         "d(x+z,(x,y,z))",
 871         "(1,0,1)",
 872
 873         "d(cos(theta)^2,cos(theta))",
 874         "2*cos(theta)",
 875
 876         "d(f())",
 877         "d(f(),x)",
 878
 879         "d(x^2)",
 880         "2*x",
 881
 882         "d(t^2)",
 883         "2*t",
 884
 885         "d(t^2 x^2)",
 886         "2*t^2*x",
 887
 888 // trig functions
 889
 890         "d(sin(x),x)-cos(x)",
 891         "0",
 892
 893         "d(cos(x),x)+sin(x)",
 894         "0",
 895
 896         "d(tan(x),x)-cos(x)^(-2)",
 897         "0",
 898
 899         "d(arcsin(x),x)-1/sqrt(1-x^2)",
 900         "0",
 901
 902         "d(arccos(x),x)+1/sqrt(1-x^2)",
 903         "0",
 904
 905         "d(arctan(x),x)-1/(1+x^2)",
 906         "0",
 907
 908         "d(arctan(y/x),x)",
 909         "-y/(x^2+y^2)",
 910
 911         "d(arctan(y/x),y)",
 912         "x/(x^2+y^2)",
 913
 914 // hyp functions
 915
 916         "d(sinh(x),x)-cosh(x)",
 917         "0",
 918
 919         "d(cosh(x),x)-sinh(x)",
 920         "0",
 921
 922         "d(tanh(x),x)-cosh(x)^(-2)",
 923         "0",
 924
 925         "d(arcsinh(x),x)-1/sqrt(x^2+1)",
 926         "0",
 927
 928         "d(arccosh(x),x)-1/sqrt(x^2-1)",
 929         "0",
 930
 931         "d(arctanh(x),x)-1/(1-x^2)",
 932         "0",
 933
 934         "d(sin(cos(x)),x)+cos(cos(x))*sin(x)",
 935         "0",
 936
 937         "d(sin(x)^2,x)-2*sin(x)*cos(x)",
 938         "0",
 939
 940         "d(sin(cos(x)),cos(x))-cos(cos(x))",
 941         "0",
 942
 943         "d(abs(x),x)",
 944         "sgn(x)",
 945
 946         "d(sgn(x),x)",
 947         "2*dirac(x)",
 948
 949 // generic functions
 950
 951         "d(f(),x)",
 952         "d(f(),x)",
 953
 954         "d(f(x),x)",
 955         "d(f(x),x)",
 956
 957         "d(f(y),x)",
 958         "0",
 959
 960         "d(g(f(x)),f(x))",
 961         "d(g(f(x)),f(x))",
 962
 963         "d(g(f(x)),x)",
 964         "d(g(f(x)),x)",
 965
 966 // other functions
 967
 968         "d(erf(x))-2*exp(-x^2)/sqrt(pi)",
 969         "0",
 970
 971 // arg lists
 972
 973         "f=x^5*y^7",
 974         "",
 975
 976         "d(f)",
 977         "5*x^4*y^7",
 978
 979         "d(f,x)",
 980         "5*x^4*y^7",
 981
 982         "d(f,x,0)",
 983         "x^5*y^7",
 984
 985         "d(f,x,1)",
 986         "5*x^4*y^7",
 987
 988         "d(f,x,2)",
 989         "20*x^3*y^7",
 990
 991         "d(f,2)",
 992         "20*x^3*y^7",
 993
 994         "d(f,2,y)",
 995         "140*x^3*y^6",
 996
 997         "d(f,x,x,y,y)",
 998         "840*x^3*y^5",
 999
1000         "f=quote(f)",
1001         "",
1002 };
1003
1004 void
1005 test_derivative(void)
1006 {
1007         test(__FILE__, s, sizeof s / sizeof (char *));
1008 }
1009
1010 #endif