integral.cpp

   1 #include "stdafx.h"
   2 #include "defs.h"
   3
   4 #define F p3
   5 #define X p4
   6 #define N p5
   7
   8 void
   9 eval_integral(void)
  10 {
  11         int i, n;
  12
  13         // evaluate 1st arg to get function F
  14
  15         p1 = cdr(p1);
  16         push(car(p1));
  17         eval();
  18
  19         // evaluate 2nd arg and then...
  20
  21         // example              result of 2nd arg       what to do
  22         //
  23         // integral(f)          nil                     guess X, N = nil
  24         // integral(f,2)        2                       guess X, N = 2
  25         // integral(f,x)        x                       X = x, N = nil
  26         // integral(f,x,2)      x                       X = x, N = 2
  27         // integral(f,x,y)      x                       X = x, N = y
  28
  29         p1 = cdr(p1);
  30         push(car(p1));
  31         eval();
  32
  33         p2 = pop();
  34         if (p2 == symbol(NIL)) {
  35                 guess();
  36                 push(symbol(NIL));
  37         } else if (isnum(p2)) {
  38                 guess();
  39                 push(p2);
  40         } else {
  41                 push(p2);
  42                 p1 = cdr(p1);
  43                 push(car(p1));
  44                 eval();
  45         }
  46
  47         N = pop();
  48         X = pop();
  49         F = pop();
  50
  51         while (1) {
  52
  53                 // N might be a symbol instead of a number
  54
  55                 if (isnum(N)) {
  56                         push(N);
  57                         n = pop_integer();
  58                         if (n == (int) 0x80000000)
  59                                 stop("nth integral: check n");
  60                 } else
  61                         n = 1;
  62
  63                 push(F);
  64
  65                 if (n >= 0) {
  66                         for (i = 0; i < n; i++) {
  67                                 push(X);
  68                                 integral();
  69                         }
  70                 } else {
  71                         n = -n;
  72                         for (i = 0; i < n; i++) {
  73                                 push(X);
  74                                 derivative();
  75                         }
  76                 }
  77
  78                 F = pop();
  79
  80                 // if N is nil then arglist is exhausted
  81
  82                 if (N == symbol(NIL))
  83                         break;
  84
  85                 // otherwise...
  86
  87                 // N            arg1            what to do
  88                 //
  89                 // number       nil             break
  90                 // number       number          N = arg1, continue
  91                 // number       symbol          X = arg1, N = arg2, continue
  92                 //
  93                 // symbol       nil             X = N, N = nil, continue
  94                 // symbol       number          X = N, N = arg1, continue
  95                 // symbol       symbol          X = N, N = arg1, continue
  96
  97                 if (isnum(N)) {
  98                         p1 = cdr(p1);
  99                         push(car(p1));
 100                         eval();
 101                         N = pop();
 102                         if (N == symbol(NIL))
 103                                 break;          // arglist exhausted
 104                         if (isnum(N))
 105                                 ;               // N = arg1
 106                         else {
 107                                 X = N;          // X = arg1
 108                                 p1 = cdr(p1);
 109                                 push(car(p1));
 110                                 eval();
 111                                 N = pop();      // N = arg2
 112                         }
 113                 } else {
 114                         X = N;                  // X = N
 115                         p1 = cdr(p1);
 116                         push(car(p1));
 117                         eval();
 118                         N = pop();              // N = arg1
 119                 }
 120         }
 121
 122         push(F);        // final result
 123 }
 124
 125 void
 126 integral(void)
 127 {
 128         save();
 129         p2 = pop();
 130         p1 = pop();
 131         if (car(p1) == symbol(ADD))
 132                 integral_of_sum();
 133         else if (car(p1) == symbol(MULTIPLY))
 134                 integral_of_product();
 135         else
 136                 integral_of_form();
 137         p1 = pop();
 138         if (find(p1, symbol(INTEGRAL)))
 139                 stop("integral: sorry, could not find a solution");
 140         push(p1);
 141         simplify();     // polish the result
 142         eval();         // normalize the result
 143         restore();
 144 }
 145
 146 void
 147 integral_of_sum(void)
 148 {
 149         p1 = cdr(p1);
 150         push(car(p1));
 151         push(p2);
 152         integral();
 153         p1 = cdr(p1);
 154         while (iscons(p1)) {
 155                 push(car(p1));
 156                 push(p2);
 157                 integral();
 158                 add();
 159                 p1 = cdr(p1);
 160         }
 161 }
 162
 163 void
 164 integral_of_product(void)
 165 {
 166         push(p1);
 167         push(p2);
 168         partition();
 169         p1 = pop();                     // pop variable part
 170         integral_of_form();
 171         multiply();                     // multiply constant part
 172 }
 173
 174 extern char *itab[];
 175
 176 void
 177 integral_of_form(void)
 178 {
 179         push(p1);
 180         push(p2);
 181         transform(itab);
 182         p3 = pop();
 183         if (p3 == symbol(NIL)) {
 184                 push_symbol(INTEGRAL);
 185                 push(p1);
 186                 push(p2);
 187                 list(3);
 188         } else
 189                 push(p3);
 190 }