function [E1s,E1b,E1cf] = different_expint(z)
sz = size(z);
z = z(:);
# First series expansion
  gm =  0.577215664901532860606512090082402431042159335;
  E1s = -gm - log(z);
  res = -z;
  ssum = res;
  k = 1;
  fflag = ones(size(res));
  rrange = find(fflag);
  while(any(fflag))
    rrange = rrange(find(fflag));
    res_tmp = res(rrange);
    res_tmp = res_tmp*k.*(-z(rrange))./((k+1)^2);
    ssum(rrange) = ssum(rrange) + res_tmp;
    k++;
    res(rrange) = res_tmp;
    fflag = abs(res_tmp) > eps*abs(ssum(rrange));
  endwhile
  E1s = E1s - ssum;
  E1s = reshape(E1s,sz);
# Asymptotic series expansion
  z_b = z;
  E1b = exp(-z_b)./z_b;
  res = ones(size(z_b));
  oldres = res;
  ssum = res;
  k = 0;
  fflag = ones(size(z_b));
  rrange = find(fflag);
  while (k<1e3 && any(fflag))
    rrange = rrange(find(fflag));
    res_tmp = res(rrange);
    oldres_tmp = oldres(rrange);
    oldres_tmp = res_tmp;
    res_tmp = res_tmp*(k+1)./(-z_b(rrange));
    ssum(rrange) = ssum(rrange) + res_tmp;
    k++;
    res(rrange) = res_tmp;
    oldres(rrange) = oldres_tmp;
    fflag = abs(oldres_tmp)>abs(res_tmp);
  endwhile
  E1b = E1b.*ssum;
  E1b = reshape(E1b,sz);
# continued fraction, Lent'z modified algorithm
f_new = 2^(-100)*ones(size(z));
C_new = f_new;
D_new = zeros(size(z));
k = 0;
fflag = ones(size(z));
rrange = find(fflag);
Delta = ones(size(z));
C_old = ones(size(z));
D_old = ones(size(z));
f_old = ones(size(z));
while (k<1e3 && any(fflag))
  rrange = rrange(find(fflag));
  f_new_tmp = f_new(rrange);
  C_new_tmp = C_new(rrange);
  D_new_tmp = D_new(rrange);
  f_old = f_new_tmp;
  C_old = C_new_tmp;
  D_old = D_new_tmp;
  b = z(rrange)*(mod(k,2)==0) + (mod(k,2)==1);
  a = exp(-z(rrange))*(k==0) + ceil((k)/2)*(k>=1);
  D_new_tmp = b + a.*D_old;
  D_new_tmp = D_new_tmp.*(D_new_tmp~=0) + 2^(-100)*(D_new_tmp==0);
  C_new_tmp = b + a./C_old;
  C_new_tmp = C_new_tmp.*(C_new_tmp~=0) + 2^(-100)*(C_new_tmp==0);
  D_new_tmp = 1./D_new_tmp;
  Delta(rrange) = C_new_tmp.*D_new_tmp;
  f_new(rrange) = f_old.*Delta(rrange);
  C_new(rrange) = C_new_tmp;
  D_new(rrange) = D_new_tmp;
  fflag = abs(Delta(rrange)-1) > 2^(-50);
  k++;
endwhile
E1cf = reshape(f_new,sz);
endfunction
