module sfpolys; % Assorted Polynomials
% will be a package of its own one day
%
% Author: Winfried Neun, Feb 1993
%
% Revision 6. April 1995, using explicit formulae for the orthogonal
% polynomials (Abramowitz/Stegun 22.3)
% Revision December 1995 by Wolfram Koepf
% June 1996 by Wolfram Koepf : improved numeric codes
% added Fibonacci numbers and Polys, W.N, September 1998
fluid '(powlis1!*);
% Bernoulli Polynomials (see e.g. Abramowitz Stegun , chapter 23
algebraic operator bernoullip;
algebraic <<
Let { BernoulliP (~n,0) => Bernoulli n when fixp n and n >=0,
BernoulliP (~n,~x) => (for k:=0:n sum (binomial(n,k) *
Bernoulli(k) * x^(n-k)))
when fixp n and n >=0} >>;
% Euler Polynomials (see e.g. Abramowitz Stegun , chapter 23
algebraic operator EulerP ;
algebraic <<
Let { EulerP (~n,1/2) => Euler(n)/2^n when fixp n and n >=0,
EulerP (~n,~x) => (for k:=0:n sum (binomial(n,k) *
Euler(k)/2^k * (x -1/2)^(n-k)))
when fixp n and n >=0}
>>;
% Univariate orthogonal bases (for approximation etc).
% Author: H. Melenk, ZIB, Berlin
% Copyright (c): ZIB Berlin 1993, all rights resrved
algebraic procedure monomial_base(x,n);
for i:=0:n collect x**i;
algebraic procedure trigonometric_base(x,n);
1 . for i:=1:n join list(sin(i*x),cos(i*x));
algebraic procedure bernstein_base(x,n);
for i:=0:n collect
binomial(n,i)*(1-x)**(n-i)*x**i;
algebraic procedure legendre_base(x,n,a,b);
legendre_base1(x,n,{a/2-b/2 + (1+a/2+b/2)*x,1},1,a,b);
algebraic procedure legendre_base1(x,n,base,r,a,b);
if r>=n then reverse base else
legendre_base1(x,n,
(((2*r+a+b+1)*(a**2-b**2)+(2*r+a+b)*(2*r+1+a+b)*(2*r+2+a+b)*x)/
(2*(r+1)*(r+1+a+b)*(2*r+a+b))*first base -
2*(r+a)*(r+b)*(2r+2+a+b)/(2*(r+1)*(r+1+a+b)*
(2*r+a+b))*second base)
. base, r+1,a,b);
algebraic procedure laguerre_base(x,n,a);
laguerre_base1(x,n,{1-x+a,1},1,a);
algebraic procedure laguerre_base1(x,n,base,r,a);
if r>=n then reverse base else
laguerre_base1(x,n,
((1+2r-x+a)/(r+1)*first base - (r+a)/(r+1)*second base )
. base, r+1,a);
algebraic procedure hermite_base(x,n);
hermite_base1(x,n,{2*x,1},1);
algebraic procedure hermite_base1(x,n,base,r);
if r>=n then reverse base else
hermite_base1(x,n,
(2x*first base - 2r*second base)
. base, r+1);
algebraic procedure chebyshev_base_T(x,n);
chebyshev_base_T1(x,n,{x,1},1);
algebraic procedure chebyshev_base_T1(x,n,base,r);
if r>=n then reverse base else
chebyshev_base_T1(x,n,
(2x*first base - second base )
. base, r+1);
algebraic procedure chebyshev_base_U(x,n);
chebyshev_base_T1(x,n,{2x,1},1);
algebraic procedure gegenbauer_base1(x,n,base,r,a);
if r>=n then reverse base else
gegenbauer_base1(x,n,
(2*(r+a)/(r+1)*x*first base - (r+2*a-1)/(r + 1)*second base )
. base, r+1,a);
algebraic procedure gegenbauer_base(x,n,a);
gegenbauer_base1(x,n,{2*a*x,1},1,a);
algebraic <<
operator HERMITEP,JACOBIP,LEGENDREP,LEGENDREQ, !~f,
LAGUERREP,CHEBYSHEVT,CHEBYSHEVU,gegenbauerP;
let limit(~f(~n,~x),~x,~lim) => f(n,lim) when freeof (lim,infinity)
and member (f,{LEGENDREP,CHEBYSHEVT,CHEBYSHEVU,Hermitep,
laguerreP,BernoulliP,EulerP,LaguerreP});
let limit(~f(~n,~m,~x),~x,~lim) => f(n,m,lim) when freeof (lim,infinity)
and member (f,{LEGENDREP,LegendreQ,gegenbauerP,laguerreP});
let limit(~f(~n,~m,~mm,~x),~x,~lim) => f(n,m,mm,lim)
when freeof (lim,infinity) and member (f,{JacobiP});
let { % AS (22.4)
LegendreP(~n,0,0) => cos(n*Pi/2)*factorial(n)/(2^n*(factorial(n/2))^2),
% AS (8.6.1)
LegendreP(~n,~m,0) =>
2^m/sqrt(Pi)*cos((n+m)*Pi/2)*GAMMA((n+m+1)/2)/GAMMA((n-m+2)/2),
% AS (8.6.2)
LegendreQ(~n,~m,0) =>
2^(m-1)/sqrt(Pi)*sin((n+m)*Pi/2)*GAMMA((n+m+1)/2)
/GAMMA((n-m+2)/2),
% AS (8.6.1)
LegendreP(~n,0) =>
1/sqrt(Pi)*cos((n)*Pi/2)*GAMMA((n+1)/2)/GAMMA((n+2)/2),
LegendreP(~n,1) => 1,
LegendreP(~n,-1) => (-1)^n,
% AS (22.4)
GegenbauerP(~n,0,0) => 2*cos(n*Pi/2)/n,
GegenbauerP(~n,~a,0)=> cos(n*Pi/2)*GAMMA(a+n/2)
/(GAMMA(a)*factorial(n/2)),
ChebyshevT(~n,0) => cos(n*Pi/2),
ChebyshevU(~n,0) => cos(n*Pi/2),
ChebyshevT(~n,1) => 1,
ChebyshevU(~n,1) => n + 1 ,
ChebyshevT(~n,-1) => (-1)^n,
ChebyshevU(~n,-1) => (n+1)* (-1)^n,
LaguerreP(~n,~a,0) => binomial(n+a,n),
LaguerreP(~n,0) => 1,
LaguerreP(0,~x) => 1,
HermiteP(~n,0) => cos(n*Pi/2)*factorial(n)/factorial(n/2) }$
let { hermitep (~n,~x)=> (begin scalar b1,b2,bex,r;
r := 1; b1 := 2x; b2 := 1;
for i:= 1:(n-1) do <<
bex := 2x*b1 - 2*r*b2;
r := r+1; b2 := b1; b1 := bex;
>>;
return b1; end)
when fixp n and n > 0 and numberp x ,
% hermitep (~n,~x)=> sub(!=z = x,first reverse hermite_base (!=z,n))
% (factorial n * for ii:=0:floor(n/2) sum ((-1)^ii/(factorial ii *
% factorial(n -2ii)) * (2*x)^(n-2ii)))
% when fixp n and n > 0 and lisp !*rounded,
hermitep (~n,~x)=> (begin scalar k,tmp,result,Ratio,oldslash,
powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
% tmp:=subs(k=0,term);
tmp:=(2*x)**n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/4/(k+1)*(n-2*k)*(n-2*k-1)/x**2;
for k:=0:n/2 do
<<
% tmp:=tmp*Ratio;
tmp:=-tmp*1/4/(k+1)*(n-2*k)*(n-2*k-1)/x**2;
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end)
when fixp n and n > 0 ,
hermitep (0,~x)=> 1};
let{ legendreP (~n,~x) =>
% (1/2^n * for ii:=0:floor(n/2) sum (binomial(n,ii) *
% binomial(2n-2ii,n)*(-1)^ii *x^(n-2ii)))
(begin
scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
tmp:=2**(-n)*factorial(2*n)/factorial(n)**2*x**n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/2/x**2*(n-2*k-1)*(n-2*k)/(k+1)/(2*n-2*k-1);
for k:=0:n/2 do
<<
% tmp:=tmp*eval(Ratio);
tmp:=-tmp/2/x**2*(n-2*k-1)*(n-2*k)/(k+1)/(2*n-2*k-1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end) when fixp n and n > 0,
legendreP (~n,~m,~x) => (-1)^m *(1-x^2)^(m/2)*
sub(!=z = x,df(legendreP (n,!=z),!=z,m))
when fixp n and n > 0 and fixp m and m > 0,
jacobiP (~n,~a,~b,~x) =>
(1/2^n * for ii:=0:n sum (binomial(n+a,ii) *
binomial(n+b,n-ii)*(x-1)^(n-ii)*(x+1)^ii))
when fixp n and n > 0 and numberp a and a > -1
and numberp b and b > -1,
jacobiP (~n,~a,~b,~x) => sub(!=z = x
,first reverse legendre_base (!=z,n,a,b))
when fixp n and n > 0,
legendreP (0,~x) => 1,
legendreP (0,0,~x) => 1,
jacobiP (0,~a,~b,~x) => 1};
let{ laguerreP(~n,~x) => laguerreP(~n,0,~x) when fixp n and n > 0,
% (for ii:=0:n sum (binomial(n,n-ii) *
% (-1)^ii/factorial ii *x^(ii)))
% when fixp n and n > 0,
laguerreP(~n,~alpha,~x) => (begin scalar b1,b2,bex,r;
r := 1; b1 := 1-x+alpha; b2 := 1;
for i:= 1:(n-1) do <<
bex := (1+2r-x+alpha)/(r+1)*b1 -
(r+alpha)/(r+1)*b2;
r := r+1; b2 := b1; b1 := bex;
>>;
return b1; end)
when fixp n and n > 0 and numberp alpha and numberp x ,
laguerreP(~n,~alpha,~x) =>
% (for ii:=0:n sum (binomial(n+alpha,n-ii) *
% (-1)^ii/factorial ii *x^(ii)))
% when fixp n and n > 0,
(begin scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
% tmp:=subs(k=0,term);
if n=0 then return(1);
tmp:=(for j:=1:n product (j+alpha))/factorial(n);
% tmp:=prod(j+alpha,j,1,n)/factorial(n);
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/(alpha+k+1)*(n-k)*x/(k+1);
for k:=0:n do
<<
% tmp:=tmp*Ratio;
tmp:=-tmp/(alpha+k+1)*(n-k)*x/(k+1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end) when fixp n and n > 0,
laguerreP(0,~a,~x) => 1};
let {chebyshevT (~n,~x) =>
%(n/2*for ii:=0:floor(n/2) sum ((-1)^ii*factorial (n-ii-1) /
% (factorial(ii) *factorial(n -2ii))* (2*x)^(n-2ii)))
(begin
scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
if n=0 then return(1);
if n=1 then return(x);
% tmp:=subs(k=0,term);
tmp:=2**(n-1)*x**n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/4*(n-2*k)*(n-2*k-1)/x**2/(n-k-1)/(k+1);
for k:=0:n/2-1 do
<<
% tmp:=tmp*eval(Ratio);
tmp:=-tmp/4*(n-2*k)*(n-2*k-1)/x**2/(n-k-1)/(k+1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end) when fixp n and n > 0 and not numberp x,
chebyshevT (~n,~x) =>
(begin
if n=0 then return(1) else if n=1 then return(x) else
if (floor(n/2)=n/2) then return(2*ChebyshevT(n/2,x)^2-1)
else return(2*ChebyshevT((n-1)/2,x)*ChebyshevT((n+1)/2,x)-x)
end) when fixp n and n > 0 and numberp x,
chebyshevT (0,~x) => 1};
let {chebyshevU (~n,~x) =>
%(for ii:=0:floor(n/2) sum ((-1)^ii*factorial (n-ii) /
% (factorial(ii) *factorial(n -2ii))* (2*x)^(n-2ii)))
(begin
scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
if n=0 then return(1);
% tmp:=subs(k=0,term);
tmp:=2**n*x**n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/4/(n-k)*(n-2*k)*(n-2*k-1)/x**2/(k+1);
for k:=0:n/2 do
<<
% tmp:=tmp*eval(Ratio);
tmp:=-tmp/4/(n-k)*(n-2*k)*(n-2*k-1)/x**2/(k+1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end) when fixp n and n > 0 and not numberp x,
chebyshevU (~n,~x) =>
( begin
if n=0 then return(1) else if n=1 then return(2*x) else
if evenp n
then return(2*ChebyshevT(n/2,x)*ChebyshevU(n/2,x)-1)
else return(2*ChebyshevU((n-1)/2,x)*ChebyshevT((n+1)/2,x))
end) when fixp n and n > 0 and numberp x,
chebyshevU (0,~x) => 1};
let { gegenbauerP (~n,~a,~x) => (begin scalar b1,b2,bex,r;
r := 1; b1 := 2*a*x; b2 := 1;
for i:= 1:(n-1) do <<
bex := 2*(r+a)/(r+1)*x*b1 -
(r+2*a-1)/(r + 1)*b2;
r := r+1; b2 := b1; b1 := bex;
>>;
return b1; end)
when fixp n and n > 0 and numberp a and numberp x ,
gegenbauerP (~n,~a,~x) =>
% (1/Gamma(a)*for ii:=0:floor(n/2) sum
%((-1)^ii* gamma(a+n-ii)/(factorial ii *factorial(n-2ii))*
% (2*x)^(n-2ii)))
(begin scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
% tmp:=subs(k=0,term);
tmp:=(for j:=1:n product (a+j-1))/factorial(n)*2**n*x**n;
% tmp:=prod(a+j-1,j,1,n)/factorial(n)*2**n*x**n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/4/(a+n-k-1)*(n-2*k)*(n-2*k-1)/x**2/(k+1);
for k:=0:n/2 do
<<
% tmp:=tmp*eval(Ratio);
tmp:=-tmp/4/(a+n-k-1)*(n-2*k)*(n-2*k-1)/x**2/(k+1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash); % restore
return(result);
end)
when fixp n and n > 0 and not(a=0),
gegenbauerP (~n,0,~x) =>
%(for ii:=0:floor(n/2) sum
%((-1)^ii* factorial(n-ii-1)/(factorial ii *factorial(n-2ii))*
% (2*x)^(n-2ii)))
(begin
scalar k,tmp,result,Ratio,oldslash,powlis1!*;
lisp setq(oldslash,remprop('slash,'opmtch));
% tmp:=subs(k=0,term);
tmp:=2**n*x**n/n;
result:=tmp;
% Ratio:=ratio(term,k);
Ratio:=-1/4*(n-2*k)*(n-2*k-1)/x**2/(n-k-1)/(k+1);
for k:=0:n/2 do
<<
% tmp:=tmp*eval(Ratio);
tmp:=-tmp/4*(n-2*k)*(n-2*k-1)/x**2/(n-k-1)/(k+1);
result:=result+tmp;
>>;
lisp put('slash,'opmtch,oldslash);
return(result);
end) when fixp n and n > 0 ,
% gegenbauerP (~n,~a,~x) => sub(!=z = x,
% first reverse gegenbauer_base(!=z,n,a))
% when fixp n and n > 0,
gegenbauerP (0,~a,~x) => 1};
% rules for differentiation
let {% AS (8.5.4)
df(LegendreP(~a,~b,~z),z) => 1/(1-z^2)*
((a+b)*LegendreP(a-1,b,z) - a*z*LegendreP(a,b,z)),
df(LegendreP(~n,~z),z) => n/(1-z^2)*(LegendreP(n-1,z)-z*LegendreP(n,z)),
df(LegendreQ(~a,~b,~z),z) => 1/(1-z^2)*
((a+b)*LegendreQ(a-1,b,z) - a*z*LegendreQ(a,b,z)),
% AS (22.8)
df(JacobiP(~n,~a,~b,~z),z) => 1/((1-z^2)*(2*n+a+b))*
(2*(n+a)*(n+b)*JacobiP(n-1,a,b,z)+n*(a-b-(2*n+a+b)*z)
*JacobiP(n,a,b,z)),
df(GegenbauerP(~n,~a,~z),z) => 1/(1-z^2)*
((n+2*a-1)*GegenbauerP(n-1,a,z)-n*z*GegenbauerP(n,a,z)),
df(ChebyshevT(~n,~z),z) =>
1/(1-z^2)*(n*ChebyshevT(n-1,z)-n*z*ChebyshevT(n,z)),
df(ChebyshevU(~n,~z),z) => 1/(1-z^2)*
((n+1)*ChebyshevU(n-1,z)-n*z*ChebyshevU(n,z)),
df(LaguerreP(~n,~a,~z),z) =>
1/z*(-(n+a)*LaguerreP(n-1,a,z)+n*LaguerreP(n,a,z)),
df(LaguerreP(~n,~z),z) => 1/z*(-(n)*LaguerreP(n-1,z)+n*LaguerreP(n,z)),
df(HermiteP(~n,~z),z) => 2*n*HermiteP(n-1,z),
% AS (23.1.5)
df(BernoulliP(~n,~z),z) => n*BernoulliP(n-1,z),
% AS (23.1.5)
df(EulerP(~n,~z),z) => n*EulerP(n-1,z) };
>>;
% following ideas from John Abbott and Wolfram Koepf
% FIBONACCI NUMBERS (and Polynoms)
flag('(fibonacci fibonaccip),'opfn);
flag('(fibonacci),'integer);
put('fibonacci,'number!-of!-args,1);
symbolic procedure fibonacci(n);
if not fixp n then mk!*sq mksqnew list ('fibonacci , n)
else
begin integer i3,m1;
if n = 0 then return 0 else if abs(n)=1 then return 1;
if n < 0 then << m1 := -1; n := abs n >>;
i3 := fib_aux (n);
return if (m1 = -1) then << if evenp n then (-i3) else i3; >>
else i3;
end;
global '(fibonacci_alist);
symbolic << fibonacci_alist := '(( 0 . 0)
(1 . 1) (2 . 1) (3 . 2) (4 . 3) (5 . 5)
(6 . 8) (7 . 13) (8 . 21) (9 . 34)) >>;
symbolic procedure fib_aux (n);
begin scalar fi;
fi := atsoc (n,fibonacci_alist);
if fi then return cdr fi;
fi := fib_aux_aux n;
fibonacci_alist := ( n . fi) . fibonacci_alist;
return fi;
end;
symbolic procedure fib_aux_aux (n); % from Wolfram Koepf, Sep 1998
% d'apres Knuth & Ptachnik: Concrete Mathematics
if evenp n then (f*(f+2*fib_aux(n/2-1))) where f=fib_aux(n/2)
else (fib_aux ((n+1)/2)^2 + fib_aux((n-1)/2)^2);
symbolic procedure fibonaccip(n,x);
if or(not fixp n, not idp x)
then mk!*sq mksqnew ('fibonaccip . list(n,x))
else
begin integer i3,i2,i1,m1;
if n= 0 then return 0 else if n=1 then return 1;
m1 := 1;
if n < 0 then << m1 := -1; n := abs n >>;
i2 := 1; i1 :=0;
for i:=2:n do << i3 := reval list('plus,list('times, x ,i2),i1);
i1 := i2; i2 :=i3>>;
return reval (list('times,list('expt,m1,list('plus,n,1)),i3));
end;
endmodule;
end;