module specfaux; % Auxiliary functions for %Special functions package for REDUCE. % Author: Winfried Neun Feb 1993 create!-package ('(specfaux sfbinom sfpolys sfbdata), '(contrib specfn)); load!-package 'specfn; endmodule; module sfbinom; % Procedures for computing Binomial coefficients % % Author: Winfried Neun, Feb 1993 algebraic operator binomial; algebraic << let { binomial (~n,~k) => factorial n / factorial k / factorial (n-k) when fixp n and fixp k and n >= k and k >= 0, binomial (~n,~k) => gamma(n+1) / gamma (k+1) / gamma(n-k+1) } >>; % Some rules for quotients of binomials are still missing endmodule; module sfpolys; % Assorted Polynomials % will be a package of its own one day % % Author: Winfried Neun, Feb 1993 % Bernoulli Polynomials (see e.g. Abramowitz Stegun , chapter 23 algebraic operator bernoullip; algebraic << let { bernoullip (~n,~x) => (for k:=0:n sum (binomial(n,k) * bernoulli(k) * x^(n-k))) when fixp n and n >=0, bernoullip (~n,0) => bernoulli n 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,{1-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,bernsteinp,jacobip,legendrep, laguerrep,chebyshevt,chebyshevu,gegenbauerp; let {hermitep (~n,~x)=> sub(z = x,first reverse hermite_base (z,n)) when fixp n and n > 0}; let {bernsteinp (~n,~x) => sub(z = x,first bernstein_base (z,n)) when fixp n and n > 0}; let{ legendrep (~n,~x) => sub(z = x, first reverse legendre_base (z,n,0,0)) when fixp n and n > 0, legendrep (~n,~m,~x) => (-1)^m *(1-x^2)^(m/2)* df(legendrep (n,x),x,m) when fixp n and n > 0 and fixp m and m > 0, jacobip (~n,~a,~b,~x) => sub(z = x ,first reverse legendre_base (z,n,a,b)) when fixp n and n > 0}; let{ laguerrep(~n,~x) => sub(z = x,first reverse laguerre_base(z,n,0)) when fixp n and n > 0, laguerrep(~n,~a,~x) => sub(z = x, first reverse laguerre_base(z,n,a)) when fixp n and n > 0}; let {chebyshevt (~n,~x) => sub(z = x, first reverse chebyshev_base_t(z,n)) when fixp n and n > 0}; let {chebyshevu (~n,~x) => sub(z = x, first reverse chebyshev_base_u(z,n)) when fixp n and n > 0}; let {gegenbauerp (~n,~a,~x) => sub(z = x, first reverse gegenbauer_base(z,n,a)) when fixp n and n > 0}; >>; endmodule; module sfbdata; % Generate necessary data for Bernoulli computation. % Author: Winfried Neun. fluid '(compute!-bernoulli); global '(!*force); flag('(force),'switch); flag('(on),'eval); on force; symbolic macro procedure mk!-bernoulli u; <>; % When I read in save!-bernoulli the macro mk!-bernoulli() will get % expanded. This is because of the RLISP flag "*force". The effect % will be that the definition of save!-bernoulli() is in effect % just bernoulli!-alist := '((....)) symbolic procedure save!-bernoulli(); bernoulli!-alist := mk!-bernoulli(); % I want to execute save!-bernoulli() just once to initialize the % table. That way even if I am running interpreted the painfully % slow initial calculation of the table gets done only once when % I first process this chunk of code. save!-bernoulli()$ compute!-bernoulli := t; off force; remflag('(on),'eval); endmodule; end;