function [r1,r2,r3,r4,r5] = exponential_integral(z)
tol = eps/2;
gamma = 0.57721566490153286060651209008240243104215933593992; % E-M
% series (Abramowitz 5.1.11)
res = -gamma-log(z);
temp = z;
n = 1;
idx = abs(temp)/n > tol * abs(res);
while (any(idx))
  res(idx) = res(idx)+temp(idx)/n;
  n = n+1;
  temp(idx) = -(z(idx).*temp(idx))/n;
  idx = abs(temp)/n > tol * abs(res);
end
r1 = res;
% Ramanujan (wikipedia)
res = -gamma-log(z);
temp = -z/2;
res = res-exp(-z/2)*2.*temp;
n = 2;
temp = temp.*(z/(2*n));
sum2 = sum(1./(1:2:n+mod(n,2)-1));
n = 3;
idx = abs(exp(-z/2)*2.*(temp*sum2)) > tol * abs(res);
while (any(idx))
  res(idx) = res(idx)-exp(-z(idx)/2)*2.*(temp(idx)*sum2);
  temp(idx) = temp(idx).*(z(idx)/(2*n));
  sum2 = sum(1./(1:2:n+mod(n,2)-1));
  n = n+1;
  idx = abs(exp(-z/2)*2.*(temp*sum2)) > tol * abs(res);
end
r2 = res;
% continued fractions http://mathworld.wolfram.com/ContinuedFraction.html
% Abramowitz 5.1.22
p_1 = exp(-z);
q_1 = zeros(size(z));
p0 = q_1;
q0 = ones(size(z));
p = z.*p0+1*p_1;
q = z.*q0+1*q_1;
scale = max(abs([p(:);p0(:);q(:);q0(:)]));
[~,e] = log2(scale);
j = 2;
p_1 = p0/2^e;
q_1 = q0/2^e;
p0 = p/2^e;
q0 = q/2^e;
idx = abs(p0./q0-p_1./q_1) > tol * abs((p0./q0));
while(any(idx) && j < 1000)
  p(idx) = 1*p0(idx)+j/2*p_1(idx);
  q(idx) = 1*q0(idx)+j/2*q_1(idx);
  j = j+1;
  p_1(idx) = p0(idx);
  q_1(idx) = q0(idx);
  p0(idx) = p(idx);
  q0(idx) = q(idx);
  p(idx) = z(idx).*p0(idx)+(j-1)/2*p_1(idx);
  q(idx) = z(idx).*q0(idx)+(j-1)/2*q_1(idx);
  j = j+1;
  scale = max(abs([p(idx);p0(idx);q(idx);q0(idx)]));
  [~,e] = log2(scale);
  p_1(idx) = p0(idx)/2^e;
  q_1(idx) = q0(idx)/2^e;
  p0(idx) = p(idx)/2^e;
  q0(idx) = q(idx)/2^e;
  idx = abs(p0./q0-p_1./q_1) > tol * abs((p0./q0));
end
r3 = p0./q0;%*exp(-z);
% continued fraction Numerical recipes (6.3.5)
% see 5.2 (Lentz's algorithm)
p_1 = exp(-z);
q_1 = zeros(size(z));
p0 = q_1;
q0 = ones(size(z));
p = (z+1).*p0+1*p_1;
q = (z+1).*q0+1*q_1;
scale = max(abs([p(:);p0(:);q(:);q0(:)]));
[~,e] = log2(scale);
j = 2;
p_1 = p0/2^e;
q_1 = q0/2^e;
p0 = p/2^e;
q0 = q/2^e;
idx = abs(p0./q0-p_1./q_1) > tol * abs((p0./q0));
while(any(idx) && j < 1000)
  p(idx) = (z(idx)+j*2-1).*p0(idx)-(j-1)^2*p_1(idx);
  q(idx) = (z(idx)+j*2-1).*q0(idx)-(j-1)^2*q_1(idx);
  j = j+1;
  scale = max(abs([p(idx);p0(idx);q(idx);q0(idx)]));
  [~,e] = log2(scale);
  p_1(idx) = p0(idx)/2^e;
  q_1(idx) = q0(idx)/2^e;
  p0(idx) = p(idx)/2^e;
  q0(idx) = q(idx)/2^e;
  idx = abs(p0./q0-p_1./q_1) > tol * abs((p0./q0));
end
r4 = p./q;

