Work around world age error for Duals in exp_generic

Project.toml

@@ -4,12 +4,14 @@ authors = ["Chris Rackauckas <[email protected]>"]
 version = "1.7.0"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
+ForwardDiff = "0.10"
 Requires = "1.0"
 julia = "1"
 

src/ExponentialUtilities.jl

@@ -1,5 +1,5 @@
 module ExponentialUtilities
-using LinearAlgebra, SparseArrays, Printf, Requires
+using LinearAlgebra, SparseArrays, Printf, Requires, ForwardDiff
 
 """
     @diagview(A,d) -> view of the `d`th diagonal of `A`.

src/exp.jl

@@ -105,18 +105,22 @@ function _exp!(A::StridedMatrix{T}; caches=nothing) where T <: LinearAlgebra.Bla
 end
 
 """
-    exp(x, vk=Val{13}()) 
-Generic exponential function, working on any `x` for which the functions 
-`LinearAlgebra.opnorm`, `+`, `*`, `^`, and `/` (including addition with 
-UniformScaling objects) are defined. Use the argument `vk` to adjust the 
+    exp(x, vk=Val{13}())
+Generic exponential function, working on any `x` for which the functions
+`LinearAlgebra.opnorm`, `+`, `*`, `^`, and `/` (including addition with
+UniformScaling objects) are defined. Use the argument `vk` to adjust the
 number of terms used in the Pade approximants at compile time.
 
 See "The Scaling and Squaring Method for the Matrix Exponential Revisited"
 by Higham, Nicholas J. in 2005 for algorithm details.
 """
-function exp_generic(x, vk=Val{13}()) 
+function exp_generic(x, vk=Val{13}())
     nx = opnorm(x, 1)
-    s = ceil(Int, log2(nx))
+    nxl2 = log2(nx)
+    if iszero(nx)
+        return x + oneunit(x)
+    end
+    s = ceil(Int, nxl2)
     if s >= 1
         exp_generic(x/(2^s), vk)^(2^s)
     else

test/runtests.jl

@@ -14,21 +14,34 @@ using ExponentialUtilities: getH, getV, _exp!
     @test A2 ≈ expA2
 end
 
-@testset "exp_generic" begin 
+@testset "exp_generic" begin
     for n in [5, 10, 30, 50, 100, 500]
         M = rand(n, n)
         @test exp(M) ≈ exp_generic(M)
 
         M′ = M / 10opnorm(M, 1)
         @test exp(M′) ≈ exp_generic(M′)
-        
+
         N = randn(n, n)
         @test exp(N) ≈ exp_generic(N)
 
         exp(n) ≈ exp_generic(n)
     end
 end
 
+# The issue was only triggered if exp_generic had been called before loading
+# ForwardDiff so it's importnat that the loading happens after the "exp_generic"
+# testset
+using ForwardDiff
+@testset "Issue 41" begin
+    @test ForwardDiff.derivative(exp_generic, 0.1) == exp_generic(0.1)
+end
+
+@testset "Issue 42" begin
+    @test exp_generic(0.0) == 1
+    @test ForwardDiff.derivative(exp_generic, 0.0) == 1
+end
+
 @testset "Phi" begin
     # Scalar phi
     K = 4