ei Function

  pure real(wp) function ei(x)
    !! Exponential integral \(\mathrm{Ei}(x)\).
    !
    !! \(\lbrace x \in \mathbb{R} \mid x \neq 0 \rbrace\)

    real(wp), intent(in) :: x

    real(wp) :: r
    integer(i2) :: n

    if (x == 0.0_wp) then
      ei = ninf()
    else if (x < 0.0_wp) then
      ei = -e1x(-x)
    else if (x <= 40.0_wp) then
      ei = x
      r = x
      do n = 2, 101
        r = r * x * (n-1) / n**2
        ei = ei + r
        if (ismall(r, ei)) exit
      end do
      ei = ei + gm + log(x)
    else if (x <= 709.0_wp) then
      ei = 1.0_wp
      r = 1.0_wp
      do n = 1, 20
        r = r * n / x
        ei = ei + r
      end do
      ei = exp(x) * ei / x 
    else
      ei = pinf()
    end if
  end function ei