I currently compare analytical vs. numerical evaluation of the complex exponential integral and find mismatches: The imaginary part differs by $\pm \pi$ and the real part has a large error when integration does not include the singularity.
While Wikipedia and other sources like http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf state that $E_i(z)$ might be ambiguous, $E_1(z)$ should not be.
I use $E_1(z)$ as implemented as expint in MATLAB and ExpIntegralE in Mathematica - both results match. The definition is as follows:
$$ E_1(z) = \int_z^{\infty} \frac{e^{-v}}{v} \, dv $$
I then compare the two functions when $z \in \mathbb{R}$:
Ts = 1e-4;
t=-3:Ts:3;t=t(:);
figure, plot(t, real([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));legend('cumtrapz','E_1');
figure, plot(t, imag([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));legend('cumtrapz','E_1');
While the real part matches, the imaginary part has the phase jump of $-\pi$:

Similarly, I can do this when $z$ is imaginary:
tau = [ -10:Ts:-Ts Ts:Ts:7 ];
w = 1i*pi*tau;w = w(:);
figure; plot(imag(w), real([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
figure; plot(imag(w), imag([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
The phase shift is $+\pi$ (I only show the imaginary part):

Furthermore, I observe a slight divergence in the real part when I do not include the singularity (i.e., when I start integration at or above zero). This is probably because the contributions cancel when I include the the singularity but it bothers me that the difference is so huge! And what could be done about it if I just want to start the integration from zero?
tau = [ -n:Ts:0-Ts Ts:Ts:N-n-Ts ];
w = 1i*pi*gamma*tau;w = w(:);
figure; plot(imag(w), real([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
figure; plot(imag(w), imag([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
I only show the real part:
