added Logarithmic integral function

docs/src/functions_list.md

@@ -78,6 +78,7 @@ ncF
 expint
 expinti
 expintx
+li
 sinint
 cosint
 ```

docs/src/functions_overview.md

@@ -39,6 +39,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`expintx(x)`](@ref SpecialFunctions.expintx)  | scaled [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``e^z \operatorname{E}_\nu(z)`` |
 | [`sinint(x)`](@ref SpecialFunctions.sinint)    | [sine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Sine_integral) ``\operatorname{Si}(x)``      |
 | [`cosint(x)`](@ref SpecialFunctions.cosint)    | [cosine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Cosine_integral) ``\operatorname{Ci}(x)``  |
+| [`li(x)`](@ref SpecialFunctions.li)            | [logarithmic integral function](https://en.wikipedia.org/wiki/Logarithmic_integral_function) ``\operatorname{li}(x)``  |
 
 
 ## Error Functions, Dawson’s and Fresnel Integrals

src/SpecialFunctions.jl

@@ -90,6 +90,7 @@ export
     expintx,
     sinint,
     cosint,
+    li,
     lbinomial
 
 include("bessel.jl")

src/expint.jl

@@ -610,3 +610,32 @@ function expinti(x::BigFloat)
     ccall((:mpfr_eint, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, Base.MPFR.ROUNDING_MODE[])
     return z
 end
+
+##############################################################################
+# Logarithmic integral function li
+
+@doc raw"""
+    expinti(x::Real)
+
+Computes the Logarithmic integral function
+```math
+\operatorname{li}(x) = \int_{0}^x \frac{1}{\ln{t}} \mathrm{d}t,
+```
+which is equivalent to ``\operatorname{Ei}(\ln{x})`` where
+``\operatorname{Ei}`` is the `expinti` function.
+
+External links: [Wikipedia](https://en.wikipedia.org/wiki/Logarithmic_integral_function)
+"""
+function li(x::Real)
+    if x < 0
+        throw(DomainError(x, "negative argument, convert to complex first"))
+    elseif x == 0
+        return 0
+    elseif x == 1
+        return -Inf
+    elseif x == Inf
+        return Inf
+    else
+        return expinti(log(x))
+    end
+end

test/expint.jl

@@ -245,3 +245,10 @@ expinti_real(x) = invoke(expinti, Tuple{Real}, x)
         end
     end
 end
+
+@testset "Logarithmic integral function" begin
+    @test li(0.0) ≈ 0
+    @test li(1.0) ≈ -Inf
+    @test li(Inf) === Inf
+    @test li(2) ≈ 1.0451637801174927
+end