NB. Exponential Integral Function: n Ei x ---> E_n(x)
NB. If monadic, defaults to 1st exponential integral E_1(x).
NB. The order n must be a positive integer.  
NB. -------------------------------------------------------
Eit1=: _44178.5471728216546 57721.7247139444345 9938.31388962036811 1842.11088667999639 101.093806161905533 
Eit1=: Eit1, 5.03416184097568393
Eib1=: 76537.3323337613621 32597.1881290274727 6106.10794245759497 635.419418378381702 37.2298352833327428 1
Eit2=: 8.67745954838443744e_8 0.999995519301390301 11.8483105554945844 45.5930644253389823 69.9279451291003023 
Eit2=: Eit2, 42.5202034768840779 8.83671808803843939 0.40137766494066472
Eib2=: 1 12.848193537915665 56.4433569561803199 106.645183769913882 89.7311097125289802 31.4971849170440750 
Eib2=: Eib2, 3.79559003762122243 0.0908804569188869219
Eit3=: _0.999999999999973414 _34.4061995006684895 _427.532671201988539 _2396.0194324749054 _6168.852100554763512
Eit3=: Eit3, _6576.09698748021179 _2106.07737142633289 _14.89908499729483169
Eib3=: 1 36.4061995006459804 494.3450702099036451 3190.27237489543304 10337.0753085840977 16324.1453557783503
Eib3=: Eib3, 11149.775287109662 2378.13899102160221

Ei=: (1&$:) : (4 : 0)
if. y < 0 do. z=. 0 
elseif. y=0 do. z=. %_1+x
elseif. (y<4)+.(x=1)+.((y+x)<12) do.   
NB. Get E_1(x) by rational polynomials in x or 1/x. 
  if. y <: 1 do.   
t_on_b=. (Eit1&p. % Eib1&p.) y
z=. (-^. y) + t_on_b
  elseif. y < 4 do.
t_on_b=. (Eit2&p. % Eib2&p.) %y
z=. (^-y)*t_on_b
  elseif. y >: 4 do.
t_on_b=. (Eit3&p. % Eib3&p.) %y
z=. (^-y)*(%y)*(1+t_on_b%y) 
  end.
  if. x > 1 do. n=. 2   NB. Use recursion relation for n > 1.
      e=. ^-y
      whilst. x>:(n=.>:n) do.
         z=. (e - y*z)%(n-1)
      end.
  end.
elseif. 1 do. NB. For (x+n)>= 12, use 70-term continued fraction:
   fk=. x&+ (2 35)&$@, 1&+   NB. Fork to generate ...
   a=. 1,}:(,|: fk i. 35)    NB. series a: 1,n,1,n+1,2,n+2,3, ...
   b=. ,35# (1 2)$ y,1       NB. series b: x,1,x,1, ...
   r=. a % 2*/\ (1,b)        NB. form series a_i/[b_(i-1)*b_i]
   R=. (1+{.r),(%1+(n=.1){r)
      whilst. n<69 do.
      R=. R, %1+((n=.>:n){r)*({:R)  
      end.
   z=.(^-y)* +/(}:P),(-:{:P=. */\R-1) NB. Take just half the last term
end.
)