expintegral_ei() bug on convergence



Hi Raymond san,

> Yes, this was basically what I had in mind, except the condition would
have been abs(z) > 1 and,
> maybe, abs(phase(z)) < 3.1, so that we stay away from the negative axis.

The condition (and (> (abs z) 1) (< (abs (phase z)) (* pi 0.99))) did not
pass the rtest_expintegral. The condition that passes the rtest_expintegral
is:
(and (> (abs z) 2) (< (abs (phase z)) (* pi 0.9)))

If you think this is OK, then I will modify bfloat and frac versions of
expintegral-e() in a same manner and create a patch for everyone's review.

For the z on the negative real axis, I don't know the correct behavior of
the expintegral_e(). Just I experimented with couple of z values such as
z=-100, z=-200, z=-300 with expintegral_e(1, z) in Maxima and e1(z) in
mpmath in Python. They coincide with each other.

(%i52) expintegral_e1(-100.0);
(%o52) -3.141592653589793*%i-2.715552744853881e+41
(%i53) expintegral_e1(-200.0);
(%o53) -3.141592653589793*%i-3.631235233159359e+84
(%i54) expintegral_e1(-300.0);
(%o54) -3.141592653589793*%i-6.49648250808866e+127

>>> e1(-100.0)
(-2.71555274485388e+41 - 3.14159265358979j)
>>> e1(-200.0)
(-3.63123523315936e+84 - 3.14159265358979j)
>>> e1(-300.0)
(-6.49648250808867e+127 - 3.14159265358979j)

Thanks and best regards,
Yasuaki Honda, Chiba, Japan

2012?4?6?2:10 Raymond Toy <toy.raymond at gmail.com>:

>
>
> 2012/4/5 ???? <yasuaki.honda at gmail.com>
>
>> Hi Raymond san, all,
>>
>> I am looking at the expintegral.lisp line 534 where the continued
>> fraction or power series
>> is chosen.
>>
>> According to DLMF Section 6.18 (http://dlmf.nist.gov/6.18 which is an
>> updated publication of A&S),
>>
>> "For small or moderate values of x and z, the expansion in power series
>> (?6.6) or in series of spherical Bessel functions (?6.10(ii)) can be used.
>> ....
>> Lastly, the continued fraction (6.9.1) can be used if z is
>> bounded away from the origin. Convergence becomes slow when z is near the
>> negative real axis, however."
>>
>> So, I would like to propose:
>> if z is close to origin or z is close to the negative real axis, we
>> expand in power series,
>> otherwise, we expand in continued fraction. The actual patch comes below.
>>
>
> Yes, this was basically what I had in mind, except the condition would
> have been abs(z) > 1 and, maybe, abs(phase(z)) < 3.1, so that we stay away
> from the negative axis.
>
> This fixes your immediate problem, I think, but we are still left with the
> issue of what do we do if we want expintegral_e on the negative axis.  I
> haven't figured that out yet....
>
> Ray
>
>