Artifact c717bc240d81d8cc8a059cce0ba4c007ced381da05d212b3a563a3099a828920:
- Executable file
r37/packages/specfn/sfpolys.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15716) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfpolys.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15716) [annotate] [blame] [check-ins using]
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;