NB. Exponential integral function E_n(y): n EI y --> result NB. If monadic, defaults to 1st exponential integral E_1(y) NB. y may be any array, n should be a single value. EIn=: (1&$:) : (dyad define) i0=. I. 1.4>:yy=. ,y NB. indices where y<=1.4 i1=. I. 1.41.4 c0=. x&EnS i0{ yy NB. do small arguments c1=. x&EnCF i1{ yy NB. do large arguments NB. put results in original positions and restore shape: ($y)$ c1 i1} c0 i0} (#yy)#0 ) NB. For y<=1.4, we use the series expansion. NB. n EnS y (should be accurate to machine limit) NB. order "n" must be single value, y may be a vector gamma=: 0.577215664901532860606512 NB. Euler's constant EnS=: dyad define psi=. - gamma- +/%>: i. nm1=._1+ n=.x b=. (psi- ^.y)* (!nm1)%~ nm1^~ -y a=. (!k-1)* n-k=. >:i. nm1 a=. +/ a%~ (i. nm1)^~/ -y d=. k* !kk=. n+ <: k=. >: i. 16 r=. a+ b- +/ d%~ kk^~/ -y NB. remove E_1(0)=infinity for 1st integral: if. n=1 do. r=. 0 (I. y<:0)}r end. ) NB. For y>1.4 we use a continued fraction expansion. NB. Usage n EnCF y (error of a few times 10^-16). NB. order "n" must be single value, y may be a vector NB. Follows code on p 217 of Numerical Recipes EnCF=: dyad define nm1=. _1+ n=. x h=. d=. %b=. n+y c=. 1e300 i=. 0 while. (i=.>:i)<100 do. a=. -i*(i+nm1) d=. %(a*d)+ b=. b+2 c=. b+ a%c h=. h*c*d end. h* ^-y )