Unexpected result for exponentialintegralen function

Fred5962 · **Joined:** Fri Aug 02, 2024 10:53 am **Posts:** 5

The function exponentialintegralen(x) replies -1.0 when x>170.0, while the function should decrease asymptotically to zero.
There is a threshold effect at x=170.0 (replies 8,64881582108917E-77).
Thank you for any clarification.

Fred5962 · **Joined:** Fri Aug 02, 2024 10:53 am **Posts:** 5

Typical polynomial and rational approximations show the obvious asymptotical trend of the exponential integral towards zero for large values of x (references in the attached graphic file).
Why does the result of AlgLib switch to E1(x) = -1 when x>170?!

Fred5962 · **Joined:** Fri Aug 02, 2024 10:53 am **Posts:** 5

The answer is found in specialfunctions.cpp file, at line 6749:

Code:

if( ((n<0||ae_fp_less(x,(double)(0)))||[b]ae_fp_greater(x,(double)(170)))[/b]||(ae_fp_eq(x,(double)(0))&&n<2) )
    {
        result = (double)(-1);
        return result;
    }

When the function exponentialintegralen is used with invalid parameters (x<0 or n<0), it returns (-1) to show that input parameters are out of the validity domain. Correct.
But the middle expression:

Code:

ae_fp_greater(x,(double)(170))

in the test must not be here.

The correct code can be found in the CEPHES script (where MAXLOG would be 170):

Code:

/* Cephes Math Library Release 2.8:  June, 2000
   Copyright 1985, 2000 by Stephen L. Moshier */

if( x > MAXLOG )
   return( 0.0 );

Why has this original CEPHES code been changed?

Fred5962 · **Joined:** Fri Aug 02, 2024 10:53 am **Posts:** 5

When the original script of CEPHES is kept - i.e. without the test ae_fp_greater(x,(double)(170)) - the function returns zero as expected when x>738.527209850, with a double float precision (64 bits).
Convergence is reached after 5 iterations.
The test with MAXLOG = 738.527209850 would even save them...
NB : always prefer exp(-x) to 1/exp(x) for large values of x, to avoid any overflow error.

Fred5962 · **Joined:** Fri Aug 02, 2024 10:53 am **Posts:** 5

Strange experience with this AlgLib forum, since after one month, I had zero interaction with the community, even with the administrator.
The issue that I have detected and located is no more and no less than a bug. I have brought a reasonable solution with a technical reference to a robust library.
Will this be fixed in a future version?

forum.alglib.net

Forum rules

Unexpected result for exponentialintegralen function

Who is online