taylor.cpp

   1 /* Taylor expansion of a function
   2
   3         push(F)
   4         push(X)
   5         push(N)
   6         push(A)
   7         taylor()
   8 */
   9
  10 #include "stdafx.h"
  11 #include "defs.h"
  12
  13 void
  14 eval_taylor(void)
  15 {
  16         // 1st arg
  17
  18         p1 = cdr(p1);
  19         push(car(p1));
  20         eval();
  21
  22         // 2nd arg
  23
  24         p1 = cdr(p1);
  25         push(car(p1));
  26         eval();
  27         p2 = pop();
  28         if (p2 == symbol(NIL))
  29                 guess();
  30         else
  31                 push(p2);
  32
  33         // 3rd arg
  34
  35         p1 = cdr(p1);
  36         push(car(p1));
  37         eval();
  38         p2 = pop();
  39         if (p2 == symbol(NIL))
  40                 push_integer(24); // default number of terms
  41         else
  42                 push(p2);
  43
  44         // 4th arg
  45
  46         p1 = cdr(p1);
  47         push(car(p1));
  48         eval();
  49         p2 = pop();
  50         if (p2 == symbol(NIL))
  51                 push_integer(0); // default expansion point
  52         else
  53                 push(p2);
  54
  55         taylor();
  56 }
  57
  58 #define F p1
  59 #define X p2
  60 #define N p3
  61 #define A p4
  62 #define C p5
  63
  64 void
  65 taylor(void)
  66 {
  67         int i, k;
  68
  69         save();
  70
  71         A = pop();
  72         N = pop();
  73         X = pop();
  74         F = pop();
  75
  76         push(N);
  77         k = pop_integer();
  78         if (k == (int) 0x80000000) {
  79                 push_symbol(TAYLOR);
  80                 push(F);
  81                 push(X);
  82                 push(N);
  83                 push(A);
  84                 list(5);
  85                 restore();
  86                 return;
  87         }
  88
  89         push(F);        // f(a)
  90         push(X);
  91         push(A);
  92         subst();
  93         eval();
  94
  95         push_integer(1);
  96         C = pop();
  97
  98         for (i = 1; i <= k; i++) {
  99
 100                 push(F);        // f = f'
 101                 push(X);
 102                 derivative();
 103                 F = pop();
 104
 105                 if (iszero(F))
 106                         break;
 107
 108                 push(C);        // c = c * (x - a)
 109                 push(X);
 110                 push(A);
 111                 subtract();
 112                 multiply();
 113                 C = pop();
 114
 115                 push(F);        // f(a)
 116                 push(X);
 117                 push(A);
 118                 subst();
 119                 eval();
 120
 121                 push(C);
 122                 multiply();
 123                 push_integer(i);
 124                 factorial();
 125                 divide();
 126
 127                 add();
 128         }
 129
 130         restore();
 131 }
 132
 133 #if SELFTEST
 134
 135 static char *s[] = {
 136
 137         "taylor(1/(5+4*cos(x)),x,6,0)-(1/9+2/81*x^2+5/1458*x^4+49/131220*x^6)",
 138         "0",
 139
 140         "taylor(1/(5+4*cos(x)),x,6)-(1/9+2/81*x^2+5/1458*x^4+49/131220*x^6)",
 141         "0",
 142 };
 143
 144 void
 145 test_taylor(void)
 146 {
 147         test(__FILE__, s, sizeof s / sizeof (char *));
 148 }
 149
 150 #endif