Add differentials for exponential integrals (#328)

Author: devmotion

src/chainrules.jl

@@ -9,6 +9,12 @@ implemented currently:
 https://github.com/JuliaMath/SpecialFunctions.jl/issues/317
 """
 
+const INCOMPLETE_EXPINT_INFO = """
+derivatives of the exponential integral with respect to parameter `ν` are not implemented
+currently:
+https://github.com/JuliaMath/SpecialFunctions.jl/issues/321
+"""
+
 ChainRulesCore.@scalar_rule(airyai(x), airyaiprime(x))
 ChainRulesCore.@scalar_rule(airyaiprime(x), x * airyai(x))
 ChainRulesCore.@scalar_rule(airybi(x), airybiprime(x))
@@ -136,3 +142,24 @@ ChainRulesCore.@scalar_rule(
         -exp(- (x + Ω)) * x^(a - 1),
     )
 )
+
+# exponential integrals
+ChainRulesCore.@scalar_rule(expint(z), - exp(-z) / z)
+ChainRulesCore.@scalar_rule(
+    expint(ν, z),
+    (
+        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
+        - expint(ν - 1, z),
+    )
+)
+ChainRulesCore.@scalar_rule(expintx(z), Ω - inv(z))
+ChainRulesCore.@scalar_rule(
+    expintx(ν, z),
+    (
+        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
+        Ω - expintx(ν - 1, z),
+    )
+)
+ChainRulesCore.@scalar_rule(expinti(x), exp(x) / x)
+ChainRulesCore.@scalar_rule(sinint(x), sinc(x / π))
+ChainRulesCore.@scalar_rule(cosint(x), cos(x) / x)

test/chainrules.jl

@@ -110,4 +110,27 @@
             end
         end
     end
+
+    @testset "exponential integrals" begin
+        for x in (1.5, 2.5, 10.5, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
+            test_scalar(expint, x)
+            test_scalar(expintx, x)
+
+            for nu in (-1.5, 2.2, 4.0)
+                # ensure all complex if any complex for FiniteDifferences
+                _x, _nu = promote(x, nu)
+
+                test_frule(expint, _nu, _x)
+                test_rrule(expint, _nu, _x)
+
+                test_frule(expintx, _nu, _x)
+                test_rrule(expintx, _nu, _x)
+            end
+
+            isreal(x) || continue
+            test_scalar(expinti, x)
+            test_scalar(sinint, x)
+            test_scalar(cosint, x)
+        end
+    end
 end