laguerre.cpp

   1 /* Laguerre function
   2
   3 Example
   4
   5         laguerre(x,3)
   6
   7 Result
   8
   9            1   3    3   2
  10         - --- x  + --- x  - 3 x + 1
  11            6        2
  12
  13 The computation uses the following recurrence relation.
  14
  15         L(x,0,k) = 1
  16
  17         L(x,1,k) = -x + k + 1
  18
  19         n*L(x,n,k) = (2*(n-1)+1-x+k)*L(x,n-1,k) - (n-1+k)*L(x,n-2,k)
  20
  21 In the "for" loop i = n-1 so the recurrence relation becomes
  22
  23         (i+1)*L(x,n,k) = (2*i+1-x+k)*L(x,n-1,k) - (i+k)*L(x,n-2,k)
  24 */
  25
  26 #include "stdafx.h"
  27 #include "defs.h"
  28
  29 void
  30 eval_laguerre(void)
  31 {
  32         // 1st arg
  33
  34         push(cadr(p1));
  35         eval();
  36
  37         // 2nd arg
  38
  39         push(caddr(p1));
  40         eval();
  41
  42         // 3rd arg
  43
  44         push(cadddr(p1));
  45         eval();
  46
  47         p2 = pop();
  48         if (p2 == symbol(NIL))
  49                 push_integer(0);
  50         else
  51                 push(p2);
  52
  53         laguerre();
  54 }
  55
  56 #define X p1
  57 #define N p2
  58 #define K p3
  59 #define Y p4
  60 #define Y0 p5
  61 #define Y1 p6
  62
  63 void
  64 laguerre(void)
  65 {
  66         int n;
  67         save();
  68
  69         K = pop();
  70         N = pop();
  71         X = pop();
  72
  73         push(N);
  74         n = pop_integer();
  75
  76         if (n < 0) {
  77                 push_symbol(LAGUERRE);
  78                 push(X);
  79                 push(N);
  80                 push(K);
  81                 list(4);
  82                 restore();
  83                 return;
  84         }
  85
  86         if (issymbol(X))
  87                 laguerre2(n);
  88         else {
  89                 Y = X;                  // do this when X is an expr
  90                 X = symbol(SECRETX);
  91                 laguerre2(n);
  92                 X = Y;
  93                 push(symbol(SECRETX));
  94                 push(X);
  95                 subst();
  96                 eval();
  97         }
  98
  99         restore();
 100 }
 101
 102 void
 103 laguerre2(int n)
 104 {
 105         int i;
 106
 107         push_integer(1);
 108         push_integer(0);
 109
 110         Y1 = pop();
 111
 112         for (i = 0; i < n; i++) {
 113
 114                 Y0 = Y1;
 115
 116                 Y1 = pop();
 117
 118                 push_integer(2 * i + 1);
 119                 push(X);
 120                 subtract();
 121                 push(K);
 122                 add();
 123                 push(Y1);
 124                 multiply();
 125
 126                 push_integer(i);
 127                 push(K);
 128                 add();
 129                 push(Y0);
 130                 multiply();
 131
 132                 subtract();
 133
 134                 push_integer(i + 1);
 135                 divide();
 136         }
 137 }
 138
 139 #if SELFTEST
 140
 141 static char *s[] = {
 142
 143         "laguerre(x,n)",
 144         "laguerre(x,n,0)",
 145
 146         "laguerre(x,n,k)",
 147         "laguerre(x,n,k)",
 148
 149         "laguerre(x,0)-1",
 150         "0",
 151
 152         "laguerre(x,1)-(-x+1)",
 153         "0",
 154
 155         "laguerre(x,2)-1/2*(x^2-4*x+2)",
 156         "0",
 157
 158         "laguerre(x,3)-1/6*(-x^3+9*x^2-18*x+6)",
 159         "0",
 160
 161         "laguerre(x,0,k)-1",
 162         "0",
 163
 164         "laguerre(x,1,k)-(-x+k+1)",
 165         "0",
 166
 167         "laguerre(x,2,k)-1/2*(x^2-2*(k+2)*x+(k+1)*(k+2))",
 168         "0",
 169
 170         "laguerre(x,3,k)-1/6*(-x^3+3*(k+3)*x^2-3*(k+2)*(k+3)*x+(k+1)*(k+2)*(k+3))",
 171         "0",
 172
 173         "laguerre(a-b,10)-eval(subst(a-b,x,laguerre(x,10)))",
 174         "0",
 175 };
 176
 177 void
 178 test_laguerre(void)
 179 {
 180         test(__FILE__, s, sizeof s / sizeof (char *));
 181 }
 182
 183 #endif