bugGNU Octave - Bugs: bug #47738, expint is inaccurate for certain...


bug #47738: expint is inaccurate for certain imaginary inputs

Submitted by:  Colin Macdonald <cbm>
Submitted on:  Tue Apr 19 21:18:35 2016  
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: FixedAssigned to: None
Originator Name: Open/Closed: Closed
Release: devOperating System: Any

Add a New Comment (Rich MarkupRich Markup):

You are not logged in

Please log in, so followups can be emailed to you.


(Jump to the original submission Jump to the original submission)

Sun Jan 22 16:37:18 2017, comment #18:

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.

Rik <rik5>
Project Administrator
Fri Jan 13 20:49:30 2017, comment #17:

@Michele: Thanks for contributing.
@Marco: Thanks for reviewing.

I'll look it over this weekend and push it.

Rik <rik5>
Project Administrator
Fri Jan 13 09:00:24 2017, comment #16:

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.

Marco Caliari <caliari>
Project Member
Tue Jan 10 18:45:10 2017, comment #15:

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.

Michele Ginesi <m_ginesi>
Tue Jan 10 18:40:22 2017, comment #14:

Dear all,

(file #39428)

Michele Ginesi <m_ginesi>
Tue Dec 27 19:42:53 2016, comment #13:

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

Rik <rik5>
Project Administrator
Fri Dec 16 08:31:56 2016, comment #12:

@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 m-files 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 Caliari <caliari>
Project Member
Wed Dec 14 17:17:14 2016, comment #11:

@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.

Rik <rik5>
Project Administrator
Wed Dec 14 14:26:04 2016, comment #10:

Dear Marco,
I've made some comparison and these are the results: the worst (relative) error for myexpint is 4.7825e-14, while for matlab is 1.0334e-13.
In the parts of the complex plane where Matlab make a better approximation, the worst error for Matlab is 2.8209e-15, for myexpint is 4.7825e-14.
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.

Michele Ginesi <m_ginesi>
Mon Dec 5 08:25:13 2016, comment #9:

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.

Marco Caliari <caliari>
Project Member
Wed Nov 23 08:32:40 2016, comment #8:

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.

Michele Ginesi <m_ginesi>
Mon May 2 07:38:31 2016, comment #7:

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)

Marco Caliari <caliari>
Project Member
Wed Apr 27 17:35:58 2016, comment #6:

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 BSD-style 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 2e-8, which looks poor.

Colin Macdonald <cbm>
Wed Apr 20 19:16:08 2016, comment #5:

What about using (5.1.56) from this page:

It has a better range and more digits.

Doug Stewart <dastew>
Wed Apr 20 18:22:57 2016, comment #4:

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.

Mike Miller <mtmiller>
Project Administrator
Wed Apr 20 16:20:53 2016, comment #3:

This is clearly about the series expansion in expint.m being inaccurate. The relevant code from expint.m is

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.

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 1e-6. 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.

Rik <rik5>
Project Administrator
Wed Apr 20 16:07:57 2016, comment #2:

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 Objective-C files by default. Of course I just overrode that and made a local association.

Rik <rik5>
Project Administrator
Tue Apr 19 21:20:52 2016, comment #1:

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

Colin Macdonald <cbm>
Tue Apr 19 21:18:35 2016, original submission:

For some imaginary inputs, the error can be as large as 10^-6.


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.

Colin Macdonald <cbm>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):

Attached Files
file #39428:  patch_expint.txt added by m_ginesi (17kB - text/plain)
file #39051:  different_expint.m added by m_ginesi (2kB - text/x-objcsrc)
file #39052:  testdifferent.png added by m_ginesi (70kB - image/png)
file #39053:  myexpint.m added by m_ginesi (5kB - text/x-objcsrc)
file #39054:  octave_vs_matlab.png added by m_ginesi (55kB - image/png)
file #37040:  expint.ref added by caliari (65kB - application/octet-stream)
file #37041:  exponential_integral.m added by caliari (2kB - text/x-objcsrc)
file #37042:  testexpint.m added by caliari (871B - text/x-objcsrc)
file #36966:  expint_test.m added by cbm (254B - application/vnd.wolfram.mathematica.package)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by m_ginesi (Updated the item)
  • -unavailable- added by caliari (Updated the item)
  • -unavailable- added by dastew (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by cbm (Submitted the item)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 15 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sun Jan 22 16:37:18 2017rik5StatusPatch Submitted=>Fixed
    Tue Jan 10 18:40:22 2017m_ginesiAttached File-=>Added patch_expint.txt, #39428
    Tue Dec 27 19:42:53 2016rik5StatusConfirmed=>Patch Submitted
    Wed Nov 23 08:24:16 2016m_ginesiAttached File-=>Added different_expint.m, #39051
      Attached File-=>Added testdifferent.png, #39052
      Attached File-=>Added myexpint.m, #39053
      Attached File-=>Added octave_vs_matlab.png, #39054
    Mon May 2 07:38:31 2016caliariAttached File-=>Added expint.ref, #37040
      Attached File-=>Added exponential_integral.m, #37041
      Attached File-=>Added testexpint.m, #37042
    Wed Apr 20 18:22:57 2016mtmillerStatusNone=>Confirmed
    Wed Apr 20 16:21:14 2016rik5Summaryexpint is inaccurate for imaginary inputs=>expint is inaccurate for certain imaginary inputs
    Tue Apr 19 21:19:29 2016cbmAttached File-=>Added expint_err_on_imag_axis.png, #36965
      Attached File-=>Added expint_test.m, #36966

    Back to the top

    Powered by Savane 3.1-cleanup1