reduce-historical/r38/packages/specfn/sfother.red

module sfother;   % Rulesets for the Struve H and L functions, Lommel
		  % 1 and 2 functions and Whittaker M and W functions.

% Author: Chris Cannam, Nov 1992.

% The aim is to re-express in terms % of other (more `standard') special
% functions. No numerical approximation code.

% Neither imports nor exports functions.
% This module contains only rulesets.

algebraic (operator struveh, struvel);

algebraic (struve!*rules := {

df(struveh(~n,~z),z)  =>
   (2/pi) - struveh(1,z) when numberp n and n = 0,

df(struveh(~n,~x),x) => (x*StruveH(-1 + n,x)- n*StruveH(n,x))/x,

df((z**n)*struveh(~n,~z),z)  => (z**n)*struveh(n-1,z),

df((z**(-n))*struveh(~n,~z),z)  =>
   (1/(sqrt(pi)*(2**n)*gamma(n+(3/2)))) - (z**(-n))*struveh(n+1,z),


struveh(~n,~z)  =>
   ((-1)**n)*besselj(-n,z)
      when numberp n and impart n = 0
	 and n < 0 and (n*2)=floor(n*2) and not evenp floor(n*2),

struveh(~n,~z)  =>
   ((2/(pi*z))**(1/2))*(1-cos z) when numberp n and n=1/2,

struveh(~n,~z)  =>
   ((z/(pi*2))**(1/2)) * (1+(2/(z**2))) -
      ((2/(pi*z))**(1/2)) * (sin z + ((cos z)/z))
      when numberp n and n=3/2,

struveh(~n,~x) => (x*0.5)^(n+1)*struve_compute_term(n,x,h)
          when numberp x and numberp n and symbolic !*rounded,

struvel(~n,~x) => struve_compute_term(n,x,l)
          when numberp x and numberp n and symbolic !*rounded,

struvel(~n,~z)  =>
   besseli(-n,z)
      when numberp n and impart n = 0
	 and n < 0 and (n*2)=floor(n*2) and not evenp floor(n*2),

struvel(~n,~z)  =>
   -i*(e**((-i*n*pi)/2))*struveh(n,i*z) when symbolic !*complex,

df(struvel(~n,~x),x) => (x*StruveL(-1 + n,x)- n*StruveL(n,x))/x

})$

algebraic (let struve!*rules);



algebraic (operator lommel1, lommel2);

algebraic (lommel!*rules := {

lommel1(~a,~b,~z)  =>
   -(2**a)*besselj(a,z)*gamma(a+1)+z**a
      when numberp a and numberp b and a = b+1,

lommel1(~a,~b,~z)  =>
   lommel1(a,-b,z)
      when numberp b and b < 0 and a neq b and a neq (b+1),

lommel1(~a,~b,~z)  =>
   (sqrt(pi)*(2**a)*gamma((2*a + 1)/2)*struveh(a,z))/2 when a = b,

lommel2(~a,~b,~z)  => z**b when numberp a and numberp b and a = b+1,

lommel2(~a,~b,~z)  => lommel2(a,-b,z)
      when numberp b and b < 0 and a neq b and a neq (b+1),

lommel2(~a,~b,~z)  =>
   (sqrt(pi)*(2**a)*gamma((2*a + 1)/2)*(-bessely(a,z)+struveh(a,z)))/2
      when a = b

})$

algebraic (let lommel!*rules);



algebraic (operator whittakerm, whittakerw);

algebraic (whittaker!*rules := {

whittakerm(~k,~m,~z)  =>
   exp(-z/2)*(z**(1/2+m))*kummerm(1/2+m-k,1+2*m,z),

whittakerw(~k,~m,~z)  =>
   exp(-z/2)*(z**(1/2+m))*kummeru(1/2+m-k,1+2*m,z),

df(WhittakerM(~n,~m,~z),z)  => 1/(2*z)*
        ((1+2*m-2*n)*WhittakerM(n-1,m,z) + (2*n-z)*WhittakerM(n,m,z)),

df(WhittakerW(~n,~m,~z),z)  => 1/(4*z)*
        ((1-4*m^2-4*n+4*n^2)*WhittakerW(n-1,m,z)
		 + (4*n-2*z)*WhittakerW(n,m,z))
% AS (8.5.4)

})$

algebraic (let whittaker!*rules);

%Handbook of Mathematical Functions - page 496

algebraic procedure struve_compute_term(n,x,h_or_l);

begin scalar dmode!*!*;
  lisp(dmode!*!* :=  dmode!*);
 return
  begin scalar pre,term,k,precis,result,!*complex,!*rounded,
   dmode!*,expo,!*msg;
   lisp (dmode!* := dmode!*!*);
  if h_or_l = l
	then << on complex;
                off rounded;
                expo := e^(-i*n*pi/2);
                on rounded;
		return  (-i*expo*struveh(n,i*x))>>
  else <<
   pre := precision 0;
   precis := 10.0^(-pre-2);
   result := 0;

    << if n > -2 then <<k:=1, term := 2^(n+2)/(pi *
                    	(for i:= 1 :n+1 product(2i-1))) ;
                	result := term >>
	else for kk:=0:-(n+2) do << k:=kk+1;
			term := (-1)^kk*(1/2*x)^(2*kk)/
				(Gamma(kk+3/2) * Gamma(kk+n+3/2));
			result := result + term>>;
        while abs(term) > precis do
          <<  term:= term*(-0.25)*(x^2)/((k+0.5)*(k+n+0.5));
              result := result + term;
              k := k+1>>;
  >>;  >>;
  return result;
end; end;

symbolic operator struve_compute_term;

% Lambert's W  (Omega) function.
% see: "On Lambert's W function" by R. Corless, G. Gonnet et. al.
% only the principal branch is implemented

algebraic <<

    % Remove autoload properties.
 lisp null remprop('lambert_w,'simpfn); 
 lisp null remflag('(lambert_w),'full);

 operator lambert_w;

 let {	lambert_w(0) => 0,
	lambert_w(-1/e) => -1,
        sum((- ~n)^(n-1)/factorial n *~z^n,n,1,infinity) 
			=> lambert_w(z),
        df(lambert_w(~z),z) => 1/((1 + lambert_w(z))*e^lambert_w z),
	log(lambert_w(~z)) => log(z) - lambert_w z,
	e^(lambert_w ~z) => ~z/lambert_w z,
        int(lambert_w(~z),z) => z*(lambert_w z -1 +1/lambert_w z),
	lambert_w(~z) =>  num_lambert_w(z)
		when numberp z and lisp !*rounded};

procedure num_lambert_w(z);
  
     if z=0 then 0 else

       begin scalar wjnew,wj,accu,expwj,oldprec,!*complex,olddmode!*;

	 lisp setq(olddmode!* ,dmode!*);
         on complex;
         oldprec := precision 5;
         accu := 10^(- lisp !:prec!:);
         if (abs z) <= 1 then  % starting point for iteration
		 if z >= -1/e then wj := 0 else wj := log(z)
         else wj := log(z) - log(log(z));
         wjnew := 100;
         while abs(wjnew) > accu do <<
	     expwj := exp(wj);
             wjnew := - (wj*expwj -z)/
		(expwj*(wj+1)-(1/2(wj+2)*(wj*expwj -z))/(wj+1));
             wj := wj + wjnew >>;

	 precision oldprec;
         accu := 10^(- lisp !:prec!:);
         while abs(wjnew) > accu do <<
             expwj := exp(wj);
             wjnew := - (wj*expwj -z)/
                (expwj*(wj+1)-(1/2(wj+2)*(wj*expwj -z))/(wj+1));
             wj := wj + wjnew >>;
        lisp setq(dmode!*,olddmode!*);
         return wj;
    end;

>>;

endmodule;

end;