--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__expint_lentz__.cc	Thu Sep 07 16:24:15 2017 +0200
@@ -0,0 +1,109 @@
+// Copyright (C) 2017 Michele Ginesi
+//
+// This file is part of Octave.
+//
+// Octave is free software; you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 3 of the License, or
+// (at your option) any later version.
+//
+// Octave is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Octave; see the file COPYING.  If not, see
+// <http://www.gnu.org/licenses/>.
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include "defun.h"
+#include <complex>
+
+DEFUN (__expint_lentz__, args, , "Continued fraction for exponential integral function")
+{
+  int nargin = args.length ();
+
+  if (nargin != 1)
+    print_usage ();
+  else
+    {
+      octave_value x_arg = args(0);
+      if (x_arg.is_single_type ())
+        {
+          const std::complex<float> x = args(0).float_complex_value ();
+          static const std::complex<float> tiny = pow (2, -100);
+          static const float eps = std::numeric_limits<float>::epsilon();
+          std::complex<float> cone(1.0, 0.0);
+          std::complex<float> czero(0.0, 0.0);
+          std::complex<float> f = tiny;
+          std::complex<float> C = f;
+          std::complex<float> D = czero;
+          std::complex<float> alpha_m = cone;
+          std::complex<float> beta_m = x;
+          std::complex<float> Delta = czero;
+          int m = 1;
+          int maxit = 100;
+          while((std::abs(Delta - cone) > eps) & (m < maxit))
+            {
+               D = beta_m + alpha_m * D;
+               if (D == czero)
+                 D = tiny;
+               C = beta_m + alpha_m / C;
+               if (C == czero)
+                 C = tiny;
+               D = cone / D;
+               Delta = C * D;
+               f *= Delta;
+               alpha_m = floor ((m + 1) / 2);
+               if ((m % 2) == 0)
+                 beta_m = x;
+               else
+                 beta_m = cone;
+               m++;
+             }
+           if (! error_state)
+            return octave_value (f);
+        }
+      else
+        {
+          const std::complex<double> x = args(0).complex_value ();
+          static const std::complex<double> tiny = pow (2, -100);
+          static const double eps = std::numeric_limits<double>::epsilon();
+          std::complex<double> cone(1.0, 0.0);
+          std::complex<double> czero(0.0, 0.0);
+          std::complex<double> f = tiny;
+          std::complex<double> C = f;
+          std::complex<double> D = czero;
+          std::complex<double> alpha_m = cone;
+          std::complex<double> beta_m = x;
+          std::complex<double> Delta = czero;
+          int m = 1;
+          int maxit = 200;
+          while((std::abs(Delta - cone) > eps) & (m < maxit))
+            {
+               D = beta_m + alpha_m * D;
+               if (D == czero)
+                 D = tiny;
+               C = beta_m + alpha_m / C;
+               if (C == czero)
+                 C = tiny;
+               D = cone / D;
+               Delta = C * D;
+               f *= Delta;
+               alpha_m = floor ((m + 1) / 2);
+               if ((m % 2) == 0)
+                 beta_m = x;
+               else
+                 beta_m = cone;
+               m++;
+             }
+           if (! error_state)
+            return octave_value (f);
+        }
+      }
+  return octave_value_list ();
+}

--- a/libinterp/corefcn/module.mk	Thu Sep 07 16:13:16 2017 +0200
+++ b/libinterp/corefcn/module.mk	Thu Sep 07 16:24:15 2017 +0200
@@ -110,6 +110,7 @@
   %reldir%/__betainc_lentz__.cc \
   %reldir%/__contourc__.cc \
   %reldir%/__dsearchn__.cc \
+  %reldir%/__expint_lentz__.cc \
   %reldir%/__gammainc_lentz__.cc \
   %reldir%/__ichol__.cc \
   %reldir%/__ilu__.cc \

--- a/scripts/specfun/expint.m	Thu Sep 07 16:13:16 2017 +0200
+++ b/scripts/specfun/expint.m	Thu Sep 07 16:24:15 2017 +0200
@@ -75,6 +75,19 @@
 ## @ifnottex
 ## @w{@code{E_1 (-x) = -Ei (x) - i*pi}}.
 ## @end ifnottex
+##
+## References:
+##
+## @nospell{M. Abramowitz and I.A. Stegun},
+## @cite{Handbook of Mathematical Functions}
+## 1964.
+##
+## @nospell{N. Bleistein and R.A. Handelsman},
+## @cite{Asymptotic expansions of integrals}
+## 1986.
+##
+## @seealso{exp}
+##
 ## @end deftypefn
 
 function E1 = expint (x)
@@ -143,39 +156,22 @@
   ## Abramowitz, Stegun, "Handbook of Mathematical Functions",
   ## formula 5.1.22, p 229.
   ## modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165.
-  f_new = 2^-100 * ones (size (x_cf), class (x_cf));
-  C_new = f_new;
-  D_new = zeros (size (x_cf), class (x_cf));
-  k = 0;
-  fflag = true (size (x_cf));
-  Delta = C_old = D_old = f_old = ones (size (x_cf), class (x_cf));
-  while (k < 1e3 && any (fflag))
-    x_cf_tmp = x_cf(fflag);
-    C_new_tmp = C_new(fflag);
-    D_new_tmp = D_new(fflag);
-    f_old = f_new(fflag);
-    C_old = C_new_tmp;
-    D_old = D_new_tmp;
-    b = x_cf_tmp*(mod (k,2) == 0) + (mod (k,2) == 1);
-    a = exp (-x_cf_tmp)*(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(fflag) = C_new_tmp.*D_new_tmp;
-    x_cf(fflag) = x_cf_tmp;
-    f_new(fflag) = f_old.*Delta(fflag);
-    C_new(fflag) = C_new_tmp;
-    D_new(fflag) = D_new_tmp;
-    fflag = abs (Delta-1) > tol;
-    k += 1;
-  endwhile
+
+  e1_cf = exp(- x_cf);
+
+  for ii = 1:length(x_cf)
+    e1_cf (ii) *= __expint_lentz__ (x_cf(ii));
+  endfor
 
-  e1_cf = f_new;
-  ## Asymptotic series, from Fikioris, Tastsoglou, Bakas,
-  ## "Selected Asymptotic Methods with Application to Magnetics and Antennas"
-  ## formula A.10 p 161.
+  # Remove spurious imaginary part if needed
+
+  if (isreal (x_cf) && x_cf >= 0)
+    e1_cf = real (e1_cf);
+  endif
+
+  ## Asymptotic series, from N. Bleistein and R.A. Handelsman
+  ## "Asymptotic expansion of integrals"
+  ## pages 1 -- 4.
   e1_a = exp (-x_a) ./ x_a;
   oldres = ssum = res = ones (size (x_a), class (x_a));
   k = 0;
@@ -189,7 +185,7 @@
     k += 1;
     res(fflag) = res_tmp;
     oldres(fflag) = oldres_tmp;
-    fflag = abs (oldres) > abs (res);
+    fflag = abs (x_a) > k;
   endwhile
   e1_a .*= ssum;
 
@@ -249,7 +245,7 @@
 %!            9.40470496160143e-05 - 2.44265223110736e-04i, ...
 %!            6.64487526601190e-09 - 7.87242868014498e-09i, ...
 %!            3.10273337426175e-13 - 2.28030229776792e-13i];
-%! assert (expint (X), y_exp, -1e-12);
+%! assert (expint (X), y_exp, -1e-14);
 
 ## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078)
 %!test
@@ -270,6 +266,10 @@
 %!assert (class (expint (uint64 (1))), "double")
 %!assert (issparse (expint (sparse (1))))
 
+## Test on the correct Image set
+%!assert (isreal (expint (linspace (0, 100))))
+%!assert (!isreal (expint (-1)))
+
 ## Test input validation
 %!error expint ()
 %!error expint (1,2)