I pushed the cset here (http://hg.savannah.gnu.org/hgweb/octave/rev/c56d30ea6cd4). It took a while to review because there were actually a lot differences from the Octave coding guidelines. Some examples for the next time to look : single quotes instead of double quotes, '~' instead of '!' for logical not operator, space between function name and opening parenthesis to differentiate between indexing and function calls, parentheses around the conditional in an if clause, and long lines greater than 80 characters.

@Michele: Thanks for contributing.
@Marco: Thanks for reviewing.
I'll look it over this weekend and push it.

The patch follows coding conventions now, as far as I know. Moreover the class of the input is preserved now. I think the patch is ready to be pushed.

Dear all,
in file #39428 there is my patch to expint. I think that all (or at least most of) the points to correct are OK now. Let me know if there is still something to correct.

Dear all,
(file #39428)

@Michele: Can you update your patch to use Octave coding conventions? Guidelines are at http://wiki.octave.org/Contribution_guidelines. The guidelines specifically for mfiles are at http://wiki.octave.org/Octave_style_guide.

@Rik: from the mathematical and algorithmic point of view, I think Michele made a good job. He should add a little bit of documentation and better follow the style of mfiles in Octave. After that, maybe with the help of some active developer in chat, he could try to make a proper patch for the development branch of Octave.

@Marco: Can you comment on myexpint by Michele? The accuracy and improvement in running time make it seem like a good idea to move the function onto the development branch of Octave.

Dear Marco,
I've made some comparison and these are the results: the worst (relative) error for myexpint is 4.7825e14, while for matlab is 1.0334e13.
In the parts of the complex plane where Matlab make a better approximation, the worst error for Matlab is 2.8209e15, for myexpint is 4.7825e14.
Talking about the execution time, Octave's expint needs 5.69309 sec, Matlab's expint needs 0.515487 sec and myexpint needs 0.420296 sec.

Dear Michele,
in file #39054, what do you mean by saying Matlab is better than Octave? How much better? Moreover, can you test the efficiency of your code by comparing it (tic toc) with the current code and the code in Matlab? The discussion in bug #47800 may help you in improving, if needed, your code.

Dear all,
I also made some experiment with this function. In file #39051 I used three strategies to evaluate expint: series expansion, asymptotic series expansion and continued fraction. The file #39052 show which strategy make the less error (exact values were computed with the octave symbolic package) in the complex square [100,100]x[100,100].
Then I wrote a function (file #39053) which divides the input and uses a different strategy depending on where the input is in the complex plane. In file #39054 there is a comparison between myexpint.m and the Matlab expint.

Dear all,
some time ago I made some experiments with expint. I implemented four formulas that I found on Abramowitz&Stegun and Numerical Recipes in Fortran77. The result was that the classical series expansion (Abramowitz 5.1.11) was the most accurate in the region around the segment [40,0] in the complex plane. The continued fraction expansion (Abramowitz 5.1.22) or its modification (Numerical Recipes) was the most accurate elsewhere. In the computational algorithm proposed by Numerical Recipes, they use the continued fraction expansion and, if it fails, they use the series.
I do not remember the details now, but I attach my experiments. exponential_integral.m computes the values in the four ways, expint.ref contains reference values (Octave's binary format) and testexpint.m compare the results.
(file #37040, file #37041, file #37042)

SciPy uses the specfun Fortran library (unrelated to our specfun pkg) corresponding to the book "Computation of Special Functions" by Zhang and Jin.
Largest relative error I've seen in SciPy is 30*eps, which is still large but at least better. I'll look at the licensing (SciPy is BSDstyle so should be fine). Or I can just reimplement the algorithm from their book. But it has a 500 loop so might be better to compile it.
I need to look further into this anyway, because sinint and cosint are also poor (cosint uses same implementation as expint but sinint has its own much worse implementation which blows up for x>300: https://github.com/cbm755/octsympy/issues/416)
I haven't tested Boost (@Mike do you know what accuracy they claim? Does it work for complex inputs?)
@Doug: it looks like the error could be as large as 2e8, which looks poor.

What about using (5.1.56) from this page:
https://books.google.ca/books?id=KiPCAgAAQBAJ&pg=PA231&lpg=PA231&dq=4.03640;+1.15198;+5.03637;+4.19160;&source=bl&ots=kfpQohRmno&sig=Y_P8kcXMiLXcs7I8Ifxc8s4Hcf0&hl=en&sa=X&ved=0ahUKEwi1vaTQ9Z3MAhXptoMKHcl2BToQ6AEIJTAC#v=onepage&q=4.03640%3B%201.15198%3B%205.03637%3B%204.19160%3B&f=false
It has a better range and more digits.

I found a reference to this expansion in Abramowitz & Stegun with exactly those terms to six significant figures (5.1.55). There are other papers out there on approximations of the exponential integral. Boost uses http://www.jstor.org/stable/2004541?seq=1#page_scan_tab_contents.

This is clearly about the series expansion in expint.m being inaccurate. The relevant code from expint.m is
elseif (abs (xt) < 10)
## Series Expansion for real (range [0,2]) or complex inputs (r < 10)
k = 1;
do
term = xt^k / (k*factorial (k));
y(t) += term;
until (abs (term) < eps (abs (y(t))) / 2  k++ >= 100)
y(t) = 0.57721566490153286 + log (xt) + y(t);
else
## FIXME: This expansion is accurate to only 1e13 at the beginning
## near 10+i, although it becomes more accurate as the magnitude
## of xt grows.
if (imag (xt) <= 0)
persistent a1 = 4.03640;
persistent a2 = 1.15198;
persistent b1 = 5.03637;
persistent b2 = 4.19160;
y(t) = (xt^2  a1*xt + a2) ...
/ ((xt^2  b1*xt + b2) * (xt) * exp (xt)) ...
 i*pi;
Your value of 10i is exactly the worst input for this function according to the FIXME note. You can get a sense of this by using a value that is close, but goes through a different code path. See below.
exp = 0.0454564330044554 + 0.0875512674239774i;
format long
x = 10i  eps(10)*i
x = 0.000000000000000 + 9.999999999999998i
obs = expint (x)
obs = 0.0454564330044520 + 0.0875512674239927i
err = exp  obs
err = 3.38618022510673e15  1.52655665885959e14i
abs (err)
ans = 1.56366153558805e14
Clearly that works.
If you look at the code for the expansion the constants a1,a2,b1,b2 are only accurate to 6 significant digits. It is therefore no surprise that the answer is also only acurate to approximately 1e6. I don't know where this expansion came from, but finding more accurate versions of the coefficients might be enough to make this work. Maybe some internet searching will yield the answer quickly.

The answer to the second question is that the '.m' file extension is not unique. On my Kubuntu 15.10 install it thinks the .m files are ObjectiveC files by default. Of course I just overrode that and made a local association.

p.s.: I wonder why does savannah thinks my .m file is a 'application/vnd.wolfram.mathematica.package'

For some imaginary inputs, the error can be as large as 10^6.
>> A = expint(10i)
A = 0.0454579311577793 + 0.0875455463635251i
>> S = double(expint(sym(10i)))
S = 0.0454564330044554 + 0.0875512674239774i
>> err = A  S
err = 1.49815332391168e06  5.72106045235166e06i
Matlab:
>> M = expint(10i)
M =
0.0454564330044556 + 0.087551267423977i
>> M  S
ans =
1.80411241501588e16  3.46944695195361e16i
I'll attach a figure of the error along line in the complex plane.
Looking at the code shows it switches between two implementations at abs(x) == 10. There is a FIXME there, but I think the inaccuracy is more serious than it suggests.
