1module sfint;     % Assorted Integral Functions, Ei, Si, Ci, Li etc.
2                  % Includes rules and limited rounded mode evaluation
3                  % for Ei, Si, si (called s_i here!), Ci, Chi, Shi,
4                  % Fresnel_S, Fresnel_C and Erf;
5                  % the numerical part to be improved!
6%
7% Author: Winfried Neun, Jun 1993
8% email : neun@ZIB.de
9
10% Redistribution and use in source and binary forms, with or without
11% modification, are permitted provided that the following conditions are met:
12%
13%    * Redistributions of source code must retain the relevant copyright
14%      notice, this list of conditions and the following disclaimer.
15%    * Redistributions in binary form must reproduce the above copyright
16%      notice, this list of conditions and the following disclaimer in the
17%      documentation and/or other materials provided with the distribution.
18%
19% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
21% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
23% CONTRIBUTORS
24% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30% POSSIBILITY OF SUCH DAMAGE.
31%
32put('Ei, 'plain!-functionsymbol, '!Ei);
33put('Ci, 'plain!-functionsymbol, '!Ci);
34put('Si, 'plain!-functionsymbol, '!Si);
35put('Chi, 'plain!-functionsymbol, '!Chi);
36put('Shi, 'plain!-functionsymbol, '!Shi);
37
38put('Ei, 'prifn, 'plain!-symbol);
39put('Si, 'prifn, 'plain!-symbol);
40put('Ci, 'prifn, 'plain!-symbol);
41put('Shi, 'prifn, 'plain!-symbol);
42put('Chi, 'prifn, 'plain!-symbol);
43
44% added Li , Winfried Neun, Oct 1998
45
46% The next 2 declarations enable better checking of number of arguments
47% by simpiden
48
49flag('(Ei Ci Si s_i Chi Shi Li erf erfc Fresnel_C Fresnel_S), 'specfn);
50
51deflist('((Ei 1) (Si 1) (s_i 1) (Ci 1) (Chi 1) (Shi 1) (Li 1)
52          (erf 1) (erfc 1) (Fresnel_S 1) (Fresnel_C 1)
53         ), 'number!-of!-args);
54
55% For Math References, see e.g. Abramowitz/Stegun, chapter 5 and 7
56
57% Exponential Integral etc.
58
59%algebraic operator Fresnel_C, Fresnel_S;%, erfc,erfi;
60
61%algebraic operator li;%, s_i, shi, chi
62
63algebraic let limit(si(~tt),~tt,infinity) => pi/2;
64
65algebraic
66let { int(sin(~~a*~tt)/~tt,~tt,0,~z) => si (a*z) when a freeof tt
67%      int(cos(~a*~x)*sin(~b*~x)/x,x) => 1/2*si(a*x+b*x)-1/2*si(a*x-b*x),
68%      int(cos(~a*~x)*sin(~x)/x,x) => 1/2*si(a*x+x)-1/2*si(a*x-x),
69%      int(cos(~x)*sin(~b*~x)/x,x) => 1/2*si(x+b*x)-1/2*si(x-b*x),
70%     int(cos(~x)*sin(~x)/x,x) => 1/2*si(2*x)
71%      Si(0) => 0,
72%      Si(-~x) => (- Si(x)),
73%      df(Si(~x),~x) => sin(x)/x  ,
74%      Si(~x) => compute!:int!:functions(x,Si)
75%                 when numberp x and lisp !*rounded
76    };
77
78algebraic
79let { int(sin(~tt)/~tt,~tt,~z,infinity) => - s_i (z),
80      limit(s_i(~tt),x,infinity) => 0
81%      s_i(~x) => si(x) - pi/2,
82%      df(s_i(~x),~x) => sin(x)/x
83    };
84
85algebraic
86let { int(exp(~tt)/~tt,~tt,-infinity,~z) =>  Ei (z)
87%      df(Ei(~x),~x) => exp(x)/x,
88%      Ei(~x) => compute!:int!:functions(x,Ei)
89%         when numberp x and abs(x) <= 20 and lisp !*rounded
90    };
91
92algebraic
93let { int(1/ln(~tt),~tt,0,~z) =>  li (z)
94%      li (~z) => Ei(log z)
95    };
96
97algebraic
98let { int(cos(~tt)/~tt,~tt,~z,infinity) => - ci (z),
99      int((cos(~tt) -1)/~tt,~tt,0,~z) => ci (z) + psi(1) -log(z)
100% psi(1) may be replaced by euler_gamma one day ...
101
102%      Ci(-~x) => - Ci(x) -i*pi,
103%      df(Ci(~x),~x) => cos(x)/x,
104%      Ci(~x) => compute!:int!:functions(x,Ci)
105%         when numberp x and abs(x) <= 20 and lisp !*rounded
106    };
107
108algebraic
109let { int(sinh(~tt)/~tt,~tt,0,~z) => shi (z)
110%      df(shi(~x),~x) => sinh(x)/x  ,
111%      shi(~x) => compute!:int!:functions(x,shi)
112%         when numberp x and abs(x) <= 20 and lisp !*rounded
113    };
114
115algebraic
116let { int((cosh(~tt) -1)/~tt,~tt,0,~z) => chi (z) + psi(1) -log(z)
117% psi(1) may be replaced by euler_gamma one day ...
118%      df(chi(~x),~x) => cosh(x)/x  ,
119%      chi(~x) => compute!:int!:functions(x,chi)
120%         when numberp x and abs(x) <= 20 and lisp !*rounded
121    };
122
123algebraic
124let { int(sin(pi/2*~tt^2),~tt,0,~z) => Fresnel_S (z),
125%      Fresnel_S(-~x) => (- Fresnel_S (x)),
126%      Fresnel_S(i* ~x) => (-i*Fresnel_S (x)),
127      limit(Fresnel_S(~tt),~tt,infinity) => 1/2
128%      df(Fresnel_S(~x),~x) => sin(pi/2*x^2) ,
129%      Fresnel_S (~x) => compute!:int!:functions(x,Fresnel_S)
130%              when numberp x and abs(x) <= 10 and lisp !*rounded
131};
132
133algebraic
134let { int(cos(pi/2*~tt^2),~tt,0,~z) => Fresnel_C (z),
135%      Fresnel_C(-~x) => (- Fresnel_C (x)),
136%      Fresnel_C(i* ~x) => (i*Fresnel_C (x)),
137      limit(Fresnel_C(~tt),~tt,infinity) => 1/2
138%      df(Fresnel_C(~x),~x) => cos(pi/2*x^2) ,
139%      Fresnel_C (~x) => compute!:int!:functions(x,Fresnel_C)
140%               when numberp x and abs(x) <= 10 and lisp !*rounded
141};
142
143algebraic
144let { limit (erf(~x),~x,infinity) => 1,
145      limit (erfc(~x),~x,infinity) => 0,
146%% Moved to alg/elem.red
147%%      erfc (~x) => 1-erf(x),
148%%      erfi(~z)  => -i * erf(i*z),
149      int(1/e^(~tt^2),~tt,0,~z) => erf(z)/2*sqrt(pi) when z freeof infinity,
150      int(1/e^(~tt^2),~tt,~z,infinity) => erfc(z)/2*sqrt(pi) when z freeof infinity
151% RmS 2013-04-09: rule moved to alg/elem.red, compute!:int!:functions autoloaded
152%      erf (~x) => compute!:int!:functions(x,erf)
153%                when numberp x and abs(x)<5 and lisp !*rounded
154    };
155
156algebraic procedure compute!:int!:functions(x,f);
157   begin scalar pre,!*uncached,scale,term,n,precis,result,interm;
158   pre := precision 0; precision pre;
159   precis := 10^(-2 * pre);
160   lisp (!*uncached := t);
161   if f = Si then
162           if x < 0 then result :=
163                        - compute!:int!:functions(-x,f) else
164            << n:=1; term := x; result := x;
165               while abs(term) > precis do
166                 << term := -1 * (term * x*x)/(2n * (2n+1));
167                    result := result + term/(2n+1);
168                    n := n + 1>>; >>
169   else if f = Ci then
170           if x=0 then rerror('specfn,0,"Ci(0) is undefined")
171            else if realvaluedp x and x < 0
172             then result := compute!:int!:functions(-x,f) + i*pi
173	    else
174              << n:=1; term := 1; result := Euler_gamma + log(x);
175               while abs(term) > precis do
176                 << term := -1 * (term * x*x)/((2n-1) * 2n);
177                    result := result + term/(2n);
178                    n := n + 1>>; >>
179    else  if f = Ei then
180          << if x=0 then rerror('specfn,0,"Ei(0) is undefined")
181              else if not realvaluedp x
182               then symbolic rerror('specfn,0,"Ei undefined for complex argument");
183             n:=1; term := 1; result := Euler_gamma + log(abs(x));
184               while abs(term) > precis do
185                 << term := (term * x)/n;
186                    result := result + term/n;
187                    n := n + 1>>; >>
188    else  if f = Shi then
189            << n:=1; term := x; result := x;
190               while abs(term) > precis do
191                 << term := (term * x*x)/(2n * (2n+1));
192                    result := result + term/(2n+1);
193                    n := n + 1>>; >>
194    else  if f = Chi then
195          if x=0 then rerror('specfn,0,"Chi(0) is undefined")
196           else if realvaluedp x and x<0
197            then result := compute!:int!:functions(-x,f) + i*pi
198	   else << n:=1; term := 1; result := Euler_gamma + log(x);
199               while abs(term) > precis do
200                 << term := (term * x*x)/((2n-1) * 2n);
201                    result := result + term/(2n);
202                    n := n + 1>>; >>
203    else if f = erf then
204           if x < 0 then result := - compute!:int!:functions(-x,f) else
205            << n:=1; term := x; result := x;
206               if floor(x*7) > pre then precision  floor(x*7);
207               interm := -1 *  x*x;
208               while abs(term) > precis do
209                 << term := (term * interm)/n;
210                    result := result + term/(2n+1);
211                    n := n + 1>>;
212               precision pre;
213               result := result*2/sqrt(pi) ;>>
214    else  if f = Fresnel_S then
215            << if x > 4.0 then precision max(pre,40);
216               if x > 6.0 then precision max(pre,80);
217               n:=1; term := x^3*pi/2; result := term/3;
218                         interm := x^4*(pi/2)^2;
219               while abs(term) > precis do
220                 << term := -1 * (term * interm)/(2n * (2n+1));
221                    result := result + term/(4n+3);
222                    n := n + 1>>; >>
223    else  if f = Fresnel_C then
224            << if x > 4.0 then precision max(pre,40);
225               if x > 6.0 then precision max(pre,80);
226               n:=1; term := x; result := x; interm := x^4*(pi/2)^2;
227               while abs(term) > precis do
228                 << term := -1 * (term * interm)/(2n * (2n-1));
229                    result := result + term/(4n+1);
230                    n := n + 1>>;  >>;
231    precision pre;
232    return result
233  end;
234
235endmodule;
236
237end;
238