File r34.3/src/specfaux.red artifact 25cd4b675b part of check-in d9e362f11e

module specfaux;  % Auxiliary functions for
%Special functions package for REDUCE.

% Author:  Winfried Neun  Feb 1993

create!-package ('(specfaux
sfbinom
sfpolys
sfbdata),
'(contrib 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;
<<for i := 1:300 do retrieve!*bern i;
list('quote, bernoulli!-alist) >>;

% 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;

REDUCE Historical
REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]