hermite.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Hermite polynomial
   4 //
   5 //      Input:          tos-2           x       (can be a symbol or expr)
   6 //
   7 //                      tos-1           n
   8 //
   9 //      Output:         Result on stack
  10 //
  11 //-----------------------------------------------------------------------------
  12
  13 #include "stdafx.h"
  14 #include "defs.h"
  15
  16 void
  17 hermite(void)
  18 {
  19         save();
  20         yyhermite();
  21         restore();
  22 }
  23
  24 // uses the recurrence relation H(x,n+1)=2*x*H(x,n)-2*n*H(x,n-1)
  25
  26 #define X p1
  27 #define N p2
  28 #define Y p3
  29 #define Y1 p4
  30 #define Y0 p5
  31
  32 void
  33 yyhermite(void)
  34 {
  35         int n;
  36
  37         N = pop();
  38         X = pop();
  39
  40         push(N);
  41         n = pop_integer();
  42
  43         if (n < 0) {
  44                 push_symbol(HERMITE);
  45                 push(X);
  46                 push(N);
  47                 list(3);
  48                 return;
  49         }
  50
  51         if (issymbol(X))
  52                 yyhermite2(n);
  53         else {
  54                 Y = X;                  // do this when X is an expr
  55                 X = symbol(SECRETX);
  56                 yyhermite2(n);
  57                 X = Y;
  58                 push(symbol(SECRETX));
  59                 push(X);
  60                 subst();
  61                 eval();
  62         }
  63 }
  64
  65 void
  66 yyhermite2(int n)
  67 {
  68         int i;
  69
  70         push_integer(1);
  71         push_integer(0);
  72
  73         Y1 = pop();
  74
  75         for (i = 0; i < n; i++) {
  76
  77                 Y0 = Y1;
  78
  79                 Y1 = pop();
  80
  81                 push(X);
  82                 push(Y1);
  83                 multiply();
  84
  85                 push_integer(i);
  86                 push(Y0);
  87                 multiply();
  88
  89                 subtract();
  90
  91                 push_integer(2);
  92                 multiply();
  93         }
  94 }
  95
  96 #if SELFTEST
  97
  98 static char *s[] = {
  99
 100         "hermite(x,n)",
 101         "hermite(x,n)",
 102
 103         "hermite(x,0)-1",
 104         "0",
 105
 106         "hermite(x,1)-2*x",
 107         "0",
 108
 109         "hermite(x,2)-(4*x^2-2)",
 110         "0",
 111
 112         "hermite(x,3)-(8*x^3-12*x)",
 113         "0",
 114
 115         "hermite(x,4)-(16*x^4-48*x^2+12)",
 116         "0",
 117
 118         "hermite(x,5)-(32*x^5-160*x^3+120*x)",
 119         "0",
 120
 121         "hermite(x,6)-(64*x^6-480*x^4+720*x^2-120)",
 122         "0",
 123
 124         "hermite(x,7)-(128*x^7-1344*x^5+3360*x^3-1680*x)",
 125         "0",
 126
 127         "hermite(x,8)-(256*x^8-3584*x^6+13440*x^4-13440*x^2+1680)",
 128         "0",
 129
 130         "hermite(x,9)-(512*x^9-9216*x^7+48384*x^5-80640*x^3+30240*x)",
 131         "0",
 132
 133         "hermite(x,10)-(1024*x^10-23040*x^8+161280*x^6-403200*x^4+302400*x^2-30240)",
 134         "0",
 135
 136         "hermite(a-b,10)-eval(subst(a-b,x,hermite(x,10)))",
 137         "0",
 138 };
 139
 140 void
 141 test_hermite(void)
 142 {
 143         test(__FILE__, s, sizeof s / sizeof (char *));
 144 }
 145
 146 #endif