Fortran implementation of exponential integral

Help
1

Looking for a Fortran version of scipy.special.expi — SciPy v1.9.0 Manual

Thanks!

2 Likes
3

Have you tried the code in TOMS 683:
https://jblevins.org/mirror/amiller/toms683.f90
It’s a double complex implementation.

You can find more in the Miller library:
https://jblevins.org/mirror/amiller/
It has served me well over the years.

4 Likes
4

Thanks for the suggestions! I found what I needed in the Miller library

1 Like
5

Though already solved, some comparison of Ei(x) for x = 1.0… (just for fun :mag_right:)

  $ gfortran-10 toms715.f90 test.f90 && ./a.out
  
  program main
    implicit none
    double precision, external :: EI
    print "(f25.17)", EI( 1.0d0 )    !!  1.89511781635593657
  end
  • C++
  $ g++ -std=c++17 test.cpp

  #include <iostream>
  #include <cmath>
  int main() {
    printf( "%25.17f\n", std::expint( 1.0 ) ); 
    // 1.89511781635593701
  }
  • Python
  > import scipy.special
  > scipy.special.expi( 1.0 )   # 1.8951178163559368
  • Julia
  > using SpecialFunctions
  > expinti( 1.0 )   # 1.8951178163559366
2 Likes
6

Oh, yeah. There’s always going to occur some round-off error.
You can verify the accuracy of the code using, e.g., MPMath (https://mpmath.org):
~>python
Python 3.10.6 (main, Aug 2 2022, 00:00:00) [GCC 12.1.1 20220507 (Red Hat 12.1.1-1)] on linux
Type “help”, “copyright”, “credits” or “license” for more information.

from mpmath import *
nprint(ei(1), n= 30)
1.89511781635593679062878891273

MPMath is a great tool to verify accuracy. An even better tool is the Arb library (https://arblib.org), which, besides giving arbitrary-precision solutions, also provides a maximum value for the error.

1 Like