reduce-historical: Artifact [4be2756ef6]

Artifact 4be2756ef6946e74509f09eaf335814b6190f7fcaa0347d3f8dab98200018bfd:

Executable file r37/packages/specfn/sfint.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: 6389) [annotate] [blame] [check-ins using] [more...]
Executable file r38/packages/specfn/sfint.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: 6389) [annotate] [blame] [check-ins using]

module sfint;     % Assorted Integral Functions, Ei, Si, Ci, Li etc.
		  % Includes rules and limited rounded mode evaluation
		  % for Ei, Si, si (called s_i here!), Ci, Chi, Shi, 
		  % Fresnel_S, Fresnel_C and Erf;
		  % the numerical part to be improved!
%
% Author: Winfried Neun, Jun 1993
% email : neun@ZIB.de

% added Li , Winfried Neun, Oct 1998

% For Math References, see e.g. Abramowitz/Stegun, chapter 5 and 7

% Exponential Integral etc.

algebraic operator Fresnel_C, Fresnel_S, erfc,erfi;

% algebraic operator Ci, Si, s_i, shi, chi;
% FJW: Ci, Si also defined in int (driver.red), so ...
symbolic((algebraic operator Ci, Si) where !*msg=nil);
algebraic operator s_i, shi, chi, li;

algebraic let limit(Si(~tt),~tt,infinity) => Pi/2;

algebraic
let { int(sin(~tt)/~tt,~tt,0,~z) => Si (z),
      int(cos(~a*~x)*sin(~b*~x)/x,x) => 1/2*Si(a*x+b*x)-1/2*Si(a*x-b*x),
      int(cos(~a*~x)*sin(~x)/x,x) => 1/2*Si(a*x+x)-1/2*Si(a*x-x),
      int(cos(~x)*sin(~b*~x)/x,x) => 1/2*Si(x+b*x)-1/2*Si(x-b*x),
      int(cos(~x)*sin(~x)/x,x) => 1/2*Si(2*x),
      Si(0) => 0,
      Si(-~x) => (- Si(x)),
      df(Si(~x),~x) => sin(x)/x  ,
      Si(~x) => compute_int_functions(x,Si)
		 when numberp x and lisp !*rounded
    };

algebraic
let { int(sin(~tt)/~tt,~tt,~z,infinity) => - s_i (z),
      limit(s_i(~tt),x,infinity) => 0,
      s_i(~x) => Si(x) - pi/2,
      df(s_i(~x),~x) => sin(x)/x
    };

algebraic
let { int(exp(~tt)/~tt,~tt,-infinity,~z) =>  Ei (z),
      df(Ei(~x),~x) => exp(x)/x,
      Ei(~x) => compute_int_functions(x,Ei)
	 when numberp x and abs(x) <= 20 and lisp !*rounded
    };

algebraic
let { int(1/ln(~tt),~tt,0,~z) =>  li (z),
      li (~z) => Ei(log z)
    };

algebraic
let { int(cos(~tt)/~tt,~tt,~z,infinity) => - Ci (z),
      int((cos(~tt) -1)/~tt,~tt,0,~z) => Ci (z) + psi(1) -log(z),
% psi(1) may be replaced by euler_gamma one day ...
      Ci(-~x) => - Ci(x) -i*pi,
      df(Ci(~x),~x) => cos(x)/x,
      Ci(~x) => compute_int_functions(x,Ci)
	 when numberp x and abs(x) <= 20 and lisp !*rounded
    };

algebraic
let { int(sinh(~tt)/~tt,~tt,0,~z) => shi (z),
      df(Shi(~x),~x) => sinh(x)/x  ,
      shi(~x) => compute_int_functions(x,shi)
	 when numberp x and abs(x) <= 20 and lisp !*rounded
    };

algebraic
let { int((cosh(~tt) -1)/~tt,~tt,0,~z) => Chi (z) + psi(1) -log(z),
% psi(1) may be replaced by euler_gamma one day ...
      df(Chi(~x),~x) => cosh(x)/x  ,
      Chi(~x) => compute_int_functions(x,Chi)
	 when numberp x and abs(x) <= 20 and lisp !*rounded
    };

algebraic
let { int(sin(Pi/2*~tt^2),~tt,0,~z) => Fresnel_S (z),
      Fresnel_S(-~x) => (- Fresnel_S (x)),
      Fresnel_S(i* ~x) => (-i*Fresnel_S (x)),
      limit(Fresnel_S(~tt),~tt,infinity) => 1/2,
      df(Fresnel_S(~x),~x) => sin(Pi/2*x^2) ,
      Fresnel_S (~x) => compute_int_functions(x,Fresnel_S) 
	      when numberp x and abs(x) <= 10 and lisp !*rounded };

algebraic
let { int(cos(Pi/2*~tt^2),~tt,0,~z) => Fresnel_C (z),
      Fresnel_C(-~x) => (- Fresnel_C (x)),
      Fresnel_C(i* ~x) => (i*Fresnel_C (x)),
      limit(Fresnel_C(~tt),~tt,infinity) => 1/2,
      df(Fresnel_C(~x),~x) => cos(Pi/2*x^2) ,
      Fresnel_C (~x) => compute_int_functions(x,Fresnel_C)
	       when numberp x and abs(x) <= 10 and lisp !*rounded };

algebraic
let { limit (erf(~x),~x,infinity) => 1,
      limit (erfc(~x),~x,infinity) => 0,
      erfc (~x) => 1-erf(x),
      erfi(~z)  => -i * erf(i*z),
      int(1/e^(~tt^2),~tt,0,~z) => erf(z)/2*sqrt(pi),
      int(1/e^(~tt^2),~tt,~z,infinity) => erfc(z)/2*sqrt(pi),
      erf (~x) => compute_int_functions(x,erf) 
		when numberp x and abs(x)<5 and lisp !*rounded };

algebraic procedure compute_int_functions(x,f);
   begin scalar pre,!*uncached,scale,term,n,precis,result,interm;
   pre := precision 0; precision pre;
   precis := 10^(-2 * pre);
   lisp (!*uncached := t);
   if f = Si then
	   if x < 0 then result := 
			- compute_int_functions(-x,f) else
	    << n:=1; term := x; result := x;
	       while abs(term) > precis do
		 << term := -1 * (term * x*x)/(2n * (2n+1));
		    result := result + term/(2n+1);
		    n := n + 1>>; >>
    else if f = Ci then
	   if x < 0 then result :=
			 - compute_int_functions(-x,f) -i*pi else
	      << n:=1; term := 1; result := euler!*constant + log(x);
	       while abs(term) > precis do
		 << term := -1 * (term * x*x)/((2n-1) * 2n);
		    result := result + term/(2n);
		    n := n + 1>>; >>
    else  if f = Ei then
	  << n:=1; term := 1; result := euler!*constant + log(x);
	       while abs(term) > precis do
		 << term := (term * x)/n;
		    result := result + term/n;
		    n := n + 1>>; >>
    else  if f = Shi then
	    << n:=1; term := x; result := x;
	       while abs(term) > precis do
		 << term := (term * x*x)/(2n * (2n+1));
		    result := result + term/(2n+1);
		    n := n + 1>>; >>
    else  if f = Chi then
	    << n:=1; term := 1; result := euler!*constant + log(x);
	       while abs(term) > precis do
		 << term := (term * x*x)/((2n-1) * 2n);
		    result := result + term/(2n);
		    n := n + 1>>; >>
    else if f = erf then
	   if x < 0 then result := - compute_int_functions(-x,f) else
	    << n:=1; term := x; result := x;
               if floor(x*7) > pre then precision  floor(x*7);
               interm := -1 *  x*x;
	       while abs(term) > precis do
		 << term := (term * interm)/n;
		    result := result + term/(2n+1);
		    n := n + 1>>;  
               precision pre;
	       result := result*2/sqrt(pi) ;>>
    else  if f = Fresnel_S then
	    << if x > 4.0 then precision max(pre,40);
	       if x > 6.0 then precision max(pre,80);
	       n:=1; term := x^3*pi/2; result := term/3;
			 interm := x^4*(pi/2)^2;
	       while abs(term) > precis do
		 << term := -1 * (term * interm)/(2n * (2n+1));
		    result := result + term/(4n+3);
		    n := n + 1>>; >>
    else  if f = Fresnel_C then
	    << if x > 4.0 then precision max(pre,40);
	       if x > 6.0 then precision max(pre,80);
	       n:=1; term := x; result := x; interm := x^4*(pi/2)^2;
	       while abs(term) > precis do
		 << term := -1 * (term * interm)/(2n * (2n-1));
		    result := result + term/(4n+1);
		    n := n + 1>>;  >>;
    precision pre;
    return result
  end;

endmodule;

end;