bug in exponential integral #4807
Description
This was reported by M. Yurko in an email to sage-support:
here are some examples of what I did:
#First, an example of the bug
Ei(20)
Output: 25615646.4145 + 6.28318530718*I
it should instead be just 25615646.4145
Ei(19)
Output: 9950907.25105
the error doesn't occur here
-exponential_integral_1(-20).n(digits=50)
Output: 25615652.664056588
here the bug doesn't occur, although the
code for exponential_integral_1 loses a lot
of accuracy by converting the number returned
from PARI into a float, so it has much less
accuracy
-pari(-20).eint1().n(digits = 50)
Output: 2.5615652664056588773746625520288944244384765625000e7 }}}
here the full accuracy from PARI is preserved
the following shows the difference in the speed
of the two methods
the current implementation
{{{
%time
for i in srange (1,10^6):
num = Ei(10)
Output: CPU time: 51.64 s, Wall time: 51.81 s
}}}
the time of PARI's implementation
{{{
%time
for i in srange (1,10^6):
num = pari(-i).eint1()
Output: CPU time: 20.12 s, Wall time: 20.32 s
}}}
PARI's implementation seems to be more than twice as fast
William Stein added this:
A quick remark: The Pari Ei only works with *real* input, whereas the scipy one works with complex input.
The Sage function will have to fixed in the meantime or
at least as a bare minimum have a big comment at the top explaining that there is a bug.
Component: calculus
Author: Michael Yurko
Reviewer: Mike Hansen
Merged: sage-4.2.alpha0
Issue created by migration from https://trac.sagemath.org/ticket/4807