expintegral_e: series failed

 Double integral of oscilatoring function
expintegral_e: series failed A.Domarkas

 Maxima sin integral funkcion expintegral_si(x)=integrate(sin(t)/t,t,0,x)
fail for large x. See
http://permalink.gmane.org/gmane.comp.mathematics.maxima.general/36455

 This can fix, if we use asymptotic expansion

Si(x)~pi/2-sin(x)/x*(1!/x-3!/x^3+5!/x^5-..)-cos(x)/x*(1!-2!/x^2+4!/x^4-..),
for x>>1
see
http://maji.utsi.edu/courses/06_690_perturbations_2/handout_p6_SpFuncExps.pdf

(%i1) si(x):=if x<0 then -si(-x)
elseif x>=0 and x<=20 then float(expintegral_si(x))
elseif x>20 and x<10^20 then
float(%pi/2-sin(x)/x*(1!/x-3!/x^3+5!/x^5-(7!/x^7))
-cos(x)/x*(1!-2!/x^2+4!/x^4-(6!/x^6)))
elseif x>=10^20 then float(%pi/2)
else expintegral_si(x)$

 Compare with function si() from Mpmath  Python library
http://code.google.com/p/mpmath/

 Examples:
(%i2) si(0);
(%o2) 0.0
(%i3) si(1);
(%o3) 0.946083070367183
(%i4) si(-1);
(%o4) -0.946083070367183
(%i5) si(7);
(%o5) 1.45459661424809
(%i6) si(10^4);
(%o6) 1.570891545385962
(%i7) si(10^22);
(%o7) 1.570796326794897
(%i8) si(-10^22);
(%o8) -1.570796326794897
(%i9) si(2+3*%i);
(%o9) expintegral_si(3*%i+2)
(%i10) float(%);
(%o10) 1.399196580646055*%i+4.547513889562289
(%i11) si(inf);
(%o11) 1.570796326794897
(%i12) si(minf);
(%o12) -1.570796326794897
(%i13) log(2);
(%o13) log(2)
(%i14) float(%), numer;
(%o14) 0.693147180559945

 Theorem.
integrate(integrate(cos(x*y*log(2)),y,1,n),x,1,n)=(si(log(2)*n^2)-2*si(log(2)*n)+si(log(2)))/log(2),
where si(x)=integrate(sin(t)/t,t,0,x).

(%i15) r:'integrate(integrate(cos(x*y*log(2)),y,1,n),x,1,n);
(%o15) integrate(sin(log(2)*n*x)/(log(2)*x)-sin(log(2)*x)/(log(2)*x),x,1,n)
(%i16) r1:changevar(r, y=x*log(2), y, x);
(%o16) integrate((sin(n*y)-sin(y))/y,y,log(2),log(2)*n)/log(2)
(%i17) declare(integrate,linear)$
(%i18) r2:expand(r1);
(%o18)
integrate(sin(n*y)/y,y,log(2),log(2)*n)/log(2)-integrate(sin(y)/y,y,log(2),log(2)*n)/log(2)
(%i19) S:changevar(part(r2,1), x=n*y, x, y)+changevar(part(r2,2), x=y, x,
y);
(%o19)
integrate(sin(x)/x,x,log(2)*n,log(2)*n^2)/log(2)-integrate(sin(x)/x,x,log(2),log(2)*n)/log(2)
 Note, that si(x) is primitive of sin(x)/x,  i.e. diff(si(x),x)=sin(x)/x .
Then by  Newton-Leibniz formula
(%i20) ('si(n^2*log(2))-'si(n*log(2))-('si(n*log(2))-'si(log(2))))/log(2);
define(F(n),(%))$
(%o20) (si(log(2)*n^2)-2*si(log(2)*n)+si(log(2)))/log(2)

 Examples:
(%i22) F(10); ev(%, nouns),numer;
(%o22) (si(100*log(2))-2*si(10*log(2))+si(log(2)))/log(2)
(%o23) -0.959759281504239
(%i24) F(10^2); ev(%, nouns),numer;
(%o24) (si(10000*log(2))-2*si(100*log(2))+si(log(2)))/log(2)
(%o25) -1.251679029234925
(%i26) F(10^3); ev(%, nouns),numer;
(%o26) (si(1000000*log(2))-2*si(1000*log(2))+si(log(2)))/log(2)
(%o27) -1.294205651997974
(%i28) F(10^4); ev(%, nouns),numer;
(%o28) (si(100000000*log(2))-2*si(10000*log(2))+si(log(2)))/log(2)
(%o29) -1.292308315915049
(%i30) F(10^40)$ ev(%, nouns),numer;
(%o31) -1.292490307176227

 Note that, the standart quadrature functions from Maxima package quadpack
,
or "quad"  from Mpmath  Python library  in this case return wrong results.

Aleksas D.