eigenmath/laguerre.cpp

184 lines
2.1 KiB
C++

/* Laguerre function
Example
laguerre(x,3)
Result
1 3 3 2
- --- x + --- x - 3 x + 1
6 2
The computation uses the following recurrence relation.
L(x,0,k) = 1
L(x,1,k) = -x + k + 1
n*L(x,n,k) = (2*(n-1)+1-x+k)*L(x,n-1,k) - (n-1+k)*L(x,n-2,k)
In the "for" loop i = n-1 so the recurrence relation becomes
(i+1)*L(x,n,k) = (2*i+1-x+k)*L(x,n-1,k) - (i+k)*L(x,n-2,k)
*/
#include "stdafx.h"
#include "defs.h"
void
eval_laguerre(void)
{
// 1st arg
push(cadr(p1));
eval();
// 2nd arg
push(caddr(p1));
eval();
// 3rd arg
push(cadddr(p1));
eval();
p2 = pop();
if (p2 == symbol(NIL))
push_integer(0);
else
push(p2);
laguerre();
}
#define X p1
#define N p2
#define K p3
#define Y p4
#define Y0 p5
#define Y1 p6
void
laguerre(void)
{
int n;
save();
K = pop();
N = pop();
X = pop();
push(N);
n = pop_integer();
if (n < 0) {
push_symbol(LAGUERRE);
push(X);
push(N);
push(K);
list(4);
restore();
return;
}
if (issymbol(X))
laguerre2(n);
else {
Y = X; // do this when X is an expr
X = symbol(SECRETX);
laguerre2(n);
X = Y;
push(symbol(SECRETX));
push(X);
subst();
eval();
}
restore();
}
void
laguerre2(int n)
{
int i;
push_integer(1);
push_integer(0);
Y1 = pop();
for (i = 0; i < n; i++) {
Y0 = Y1;
Y1 = pop();
push_integer(2 * i + 1);
push(X);
subtract();
push(K);
add();
push(Y1);
multiply();
push_integer(i);
push(K);
add();
push(Y0);
multiply();
subtract();
push_integer(i + 1);
divide();
}
}
#if SELFTEST
static char *s[] = {
"laguerre(x,n)",
"laguerre(x,n,0)",
"laguerre(x,n,k)",
"laguerre(x,n,k)",
"laguerre(x,0)-1",
"0",
"laguerre(x,1)-(-x+1)",
"0",
"laguerre(x,2)-1/2*(x^2-4*x+2)",
"0",
"laguerre(x,3)-1/6*(-x^3+9*x^2-18*x+6)",
"0",
"laguerre(x,0,k)-1",
"0",
"laguerre(x,1,k)-(-x+k+1)",
"0",
"laguerre(x,2,k)-1/2*(x^2-2*(k+2)*x+(k+1)*(k+2))",
"0",
"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))",
"0",
"laguerre(a-b,10)-eval(subst(a-b,x,laguerre(x,10)))",
"0",
};
void
test_laguerre(void)
{
test(__FILE__, s, sizeof s / sizeof (char *));
}
#endif