From 54d5d10c7b94b3e41807feb8d1f042a52ec262e6 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 26 Jun 2026 17:36:20 +0200
Subject: [PATCH 01/14] Exponential via MatrixAlgebraKit

---
 src/tensors/linalg.jl | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index 40a0e0024..c45c9f4f6 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -421,7 +421,16 @@ function exp!(t::TensorMap)
     domain(t) == codomain(t) ||
         error("Exponential of a tensor only exist when domain == codomain.")
     for (c, b) in blocks(t)
-        copy!(b, LinearAlgebra.exp!(b))
+        MatrixAlgebraKit.exponential!(b, b)
+    end
+    return t
+end
+
+function exp!((τ, t)::Tuple{Number, TensorMap})
+    domain(t) == codomain(t) ||
+        error("Exponential of a tensor only exist when domain == codomain.")
+    for (c, b) in blocks(t)
+        MatrixAlgebraKit.exponential!((τ, b), b)
     end
     return t
 end

From caab0e48fc066979c915035c3ac9991d4db4caf9 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 26 Jun 2026 17:58:15 +0200
Subject: [PATCH 02/14] include case where t is real and tau is complex

---
 src/tensors/linalg.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index c45c9f4f6..f07f1b174 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -429,6 +429,13 @@ end
 function exp!((τ, t)::Tuple{Number, TensorMap})
     domain(t) == codomain(t) ||
         error("Exponential of a tensor only exist when domain == codomain.")
+    if eltype(t) <: Real && eltype(τ) <: Complex
+        t_complex = complex(t)
+        for ((cr, br), (cc, bc)) in zip(blocks(t), blocks(t_complex))
+            MatrixAlgebraKit.exponential!((τ, br), bc)
+        end
+        return t_complex
+    end
     for (c, b) in blocks(t)
         MatrixAlgebraKit.exponential!((τ, b), b)
     end

From 863ca6d28b709f16caafa96ab866b652d53c8214 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 26 Jun 2026 18:43:27 +0200
Subject: [PATCH 03/14] Add test

---
 test/tensors/exponential.jl | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 test/tensors/exponential.jl

diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
new file mode 100644
index 000000000..cb5befb4e
--- /dev/null
+++ b/test/tensors/exponential.jl
@@ -0,0 +1,25 @@
+using Test, TestExtras
+using TensorKit
+using Random
+
+spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1//2 => 1)]
+scalartypes = [Float64, ComplexF64]
+
+@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space = spaces, st1 = scalartypes, st2 = scalartypes
+    A = randn(st1, space, space)
+    τ = rand(st2)
+
+    @test exp!(copy(A)) == exp!((1.0, copy(A)))
+
+    A2 = exp!((τ,A))
+    if st1 <: Real && st2 <: Complex
+        @test objectid(A2) != objectid(A)
+    else
+        @test objectid(A2) == objectid(A)
+    end
+
+    expτA = exp!((τ,copy(A)))
+    expminτA = exp!((-τ,copy(A)))
+    @test expτA * expminτA ≈ id(scalartype(expτA), space)
+    @test expτA ≈ inv(expminτA)
+end

From bf8c3af6a3c254d21dbd599456c3d730668b00fc Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 26 Jun 2026 18:44:06 +0200
Subject: [PATCH 04/14] Fix formatting

---
 test/tensors/exponential.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index cb5befb4e..51fe1072b 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -2,24 +2,24 @@ using Test, TestExtras
 using TensorKit
 using Random
 
-spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1//2 => 1)]
+spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
 scalartypes = [Float64, ComplexF64]
 
-@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space = spaces, st1 = scalartypes, st2 = scalartypes
+@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
     A = randn(st1, space, space)
     τ = rand(st2)
 
     @test exp!(copy(A)) == exp!((1.0, copy(A)))
 
-    A2 = exp!((τ,A))
+    A2 = exp!((τ, A))
     if st1 <: Real && st2 <: Complex
         @test objectid(A2) != objectid(A)
     else
         @test objectid(A2) == objectid(A)
     end
 
-    expτA = exp!((τ,copy(A)))
-    expminτA = exp!((-τ,copy(A)))
+    expτA = exp!((τ, copy(A)))
+    expminτA = exp!((-τ, copy(A)))
     @test expτA * expminτA ≈ id(scalartype(expτA), space)
     @test expτA ≈ inv(expminτA)
 end

From 9db5a963cb986e26a5642a7011ed934daf890b4b Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 3 Jul 2026 15:37:16 +0200
Subject: [PATCH 05/14] Update correspondence with MAK and extend tests

---
 src/factorizations/matrixalgebrakit.jl | 34 ++++++++++++++++
 test/tensors/exponential.jl            | 55 ++++++++++++++++++++++----
 2 files changed, 81 insertions(+), 8 deletions(-)

diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index 83d53fb0b..87d20c042 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -8,6 +8,7 @@ for f in
         :eig_full, :eig_vals, :eigh_full, :eigh_vals,
         :left_polar, :right_polar,
         :project_hermitian, :project_antihermitian, :project_isometric,
+        :exponential,
     ]
     f! = Symbol(f, :!)
     @eval function MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AbstractTensorMap}
@@ -49,6 +50,7 @@ for f! in (
         :qr_null!, :lq_null!,
         :svd_vals!, :eig_vals!, :eigh_vals!,
         :project_hermitian!, :project_antihermitian!, :project_isometric!,
+        :exponential!,
     )
     @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
         $(f! in (:eig_vals!, :eigh_vals!, :project_hermitian!, :project_antihermitian!) && :(LinearAlgebra.checksquare(t)))
@@ -62,6 +64,31 @@ for f! in (
     end
 end
 
+# Exponential with Tuple
+function MAK.exponential!((τ, t)::Tuple{E, T}, N, alg::AbstractAlgorithm) where {E <: Number, T <: AbstractTensorMap}
+    LinearAlgebra.checksquare(t)
+    foreachblock(t, N) do _, (tblock, Nblock)
+        Nblock′ = exponential!((τ, tblock), Nblock, alg)
+        # deal with the case where the output is not the same as the input
+        Nblock === Nblock′ || copy!(Nblock, Nblock′)
+        return nothing
+    end
+    return N
+end
+
+# Default algorithm for exponential with Tuple
+MAK.exponential!((τ, t)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} = MAK.exponential!((τ, t), DefaultAlgorithm())
+MAK.exponential!(t::Tuple{E, T}, ::DefaultAlgorithm) where {E <: Number, T <: Diagonal} = MAK.exponential!(t, DiagonalAlgorithm())
+MAK.exponential!(t::Tuple{E, T1}, out::T2, ::DefaultAlgorithm) where {E <: Number, T1 <: Diagonal, T2 <: Diagonal} = MAK.exponential!(t, out, DiagonalAlgorithm())
+
+function MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T}
+    return MAK.default_algorithm(exponential!, blocktype(T); kwargs...)
+end
+
+function MAK.copy_input(::typeof(exponential), (τ, t)::Tuple{E, T}) where {E <: Number, T <: AbstractTensorMap}
+    return (τ, copy_oftype(t, factorisation_scalartype(exponential, t)))
+end
+
 MAK.zero!(t::AbstractTensorMap) = zerovector!(t)
 
 # Default algorithm
@@ -76,6 +103,7 @@ for f in [
         :left_polar, :right_polar,
         :left_orth, :right_orth, :left_null, :right_null,
         :project_hermitian, :project_antihermitian, :project_isometric,
+        :exponential,
     ]
     f! = Symbol(f, :!)
     @eval MAK.$f!(t::AbstractTensorMap, alg::DefaultAlgorithm) =
@@ -221,3 +249,9 @@ MAK.initialize_output(::typeof(project_antihermitian!), tsrc::AbstractTensorMap,
     tsrc
 MAK.initialize_output(::typeof(project_isometric!), tsrc::AbstractTensorMap, ::AbstractAlgorithm) =
     similar(tsrc)
+
+# Exponential
+# ----------------
+MAK.initialize_output(::typeof(exponential!), t::AbstractTensorMap, ::AbstractAlgorithm) = t
+MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{Number, AbstractTensorMap}, ::AbstractAlgorithm) = t
+MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{T1, AbstractTensorMap{T2}}, ::AbstractAlgorithm) where {T1 <: Complex, T2 <: Real} = similar(complex(t))
diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index 51fe1072b..fbb2d73e9 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -1,25 +1,64 @@
 using Test, TestExtras
-using TensorKit
+using TensorKit, MatrixAlgebraKit
 using Random
 
-spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
-scalartypes = [Float64, ComplexF64]
+spaces = [ℂ^4] #, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
+scalartypes = [Float64, ComplexF32] #, ComplexF64]
+algorithms = [MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()), MatrixFunctionViaEigh(DefaultAlgorithm())]
 
-@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+@timedtestset "exponential for Hermitian matrices with $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
     A = randn(st1, space, space)
+    A = project_hermitian!(A)
     τ = rand(st2)
 
-    @test exp!(copy(A)) == exp!((1.0, copy(A)))
+    expA = @constinferred exponential(A)
+    expτA = @constinferred exponential((τ, A))
 
-    A2 = exp!((τ, A))
+    for alg in algorithms
+        expA2 = @constinferred exponential(A, alg)
+        expτA2 = @constinferred exponential((τ, A), alg)
+
+        @test expA ≈ expA2
+        @test expτA ≈ expτA2
+    end
+end
+
+@timedtestset "exponential! for general matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+    A = randn(st1, space, space)
+    τ = rand(st2)
+
+    @test exponential!(copy(A)) == exponential!((1.0, copy(A)))
+
+    A2 = exponential!((τ, A))
+    if st1 <: Real && st2 <: Complex
+        @test objectid(A2) != objectid(A)
+    else
+        @test objectid(A2) == objectid(A)
+    end
+
+    expτA = exponential!((τ, copy(A)))
+    expminτA = exponential!((-τ, copy(A)))
+    @test expτA * expminτA ≈ id(scalartype(expτA), space)
+    @test expτA ≈ inv(expminτA)
+end
+
+@timedtestset "exponential! for diagonal matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+    A = DiagonalTensorMap(rand(st1, reduceddim(space)), space)
+    τ = rand(st2)
+
+    exponential!(copy(A))
+    @test exponential!(copy(A)) ≈ exponential!((1.0, copy(A))) #, DiagonalAlgorithm())
+
+    A2 = @constinferred exponential!((τ, A))
+    @test A2 isa DiagonalTensorMap
     if st1 <: Real && st2 <: Complex
         @test objectid(A2) != objectid(A)
     else
         @test objectid(A2) == objectid(A)
     end
 
-    expτA = exp!((τ, copy(A)))
-    expminτA = exp!((-τ, copy(A)))
+    expτA = exponential!((τ, copy(A)))
+    expminτA = exponential!((-τ, copy(A)))
     @test expτA * expminτA ≈ id(scalartype(expτA), space)
     @test expτA ≈ inv(expminτA)
 end

From 6e5d094618676743d4031b687ae5b45535a7c8a7 Mon Sep 17 00:00:00 2001
From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com>
Date: Fri, 3 Jul 2026 17:35:40 +0200
Subject: [PATCH 06/14] Update src/factorizations/matrixalgebrakit.jl

Co-authored-by: Lukas Devos <ldevos98@gmail.com>
---
 src/factorizations/matrixalgebrakit.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index 87d20c042..7de5cb740 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -254,4 +254,4 @@ MAK.initialize_output(::typeof(project_isometric!), tsrc::AbstractTensorMap, ::A
 # ----------------
 MAK.initialize_output(::typeof(exponential!), t::AbstractTensorMap, ::AbstractAlgorithm) = t
 MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{Number, AbstractTensorMap}, ::AbstractAlgorithm) = t
-MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{T1, AbstractTensorMap{T2}}, ::AbstractAlgorithm) where {T1 <: Complex, T2 <: Real} = similar(complex(t))
+MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{T1, AbstractTensorMap{T2}}, ::AbstractAlgorithm) where {T1 <: Complex, T2 <: Real} = similar(t, complex(eltype(t)))

From 121de729da6deba590024571bad09665c7c87057 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Fri, 3 Jul 2026 17:37:04 +0200
Subject: [PATCH 07/14] deprecate `exp!`

---
 src/tensors/linalg.jl | 26 +-------------------------
 1 file changed, 1 insertion(+), 25 deletions(-)

diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index f07f1b174..0852d9c27 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -416,31 +416,7 @@ function Base.:(/)(t1::AbstractTensorMap, t2::AbstractTensorMap)
     return t
 end
 
-# TensorMap exponentation:
-function exp!(t::TensorMap)
-    domain(t) == codomain(t) ||
-        error("Exponential of a tensor only exist when domain == codomain.")
-    for (c, b) in blocks(t)
-        MatrixAlgebraKit.exponential!(b, b)
-    end
-    return t
-end
-
-function exp!((τ, t)::Tuple{Number, TensorMap})
-    domain(t) == codomain(t) ||
-        error("Exponential of a tensor only exist when domain == codomain.")
-    if eltype(t) <: Real && eltype(τ) <: Complex
-        t_complex = complex(t)
-        for ((cr, br), (cc, bc)) in zip(blocks(t), blocks(t_complex))
-            MatrixAlgebraKit.exponential!((τ, br), bc)
-        end
-        return t_complex
-    end
-    for (c, b) in blocks(t)
-        MatrixAlgebraKit.exponential!((τ, b), b)
-    end
-    return t
-end
+@deprecate exp!(t) exponential!(t)
 
 # Sylvester equation with TensorMap objects:
 function LinearAlgebra.sylvester(A::AbstractTensorMap, B::AbstractTensorMap, C::AbstractTensorMap)

From 624dada4be296e7ecb58e92687fb06bd3d211f74 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Wed, 8 Jul 2026 11:30:34 +0200
Subject: [PATCH 08/14] Remove extra definitions that are now in
 MatrixAlgebraKit

---
 src/factorizations/matrixalgebrakit.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index 7de5cb740..a1aa77c62 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -78,8 +78,6 @@ end
 
 # Default algorithm for exponential with Tuple
 MAK.exponential!((τ, t)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} = MAK.exponential!((τ, t), DefaultAlgorithm())
-MAK.exponential!(t::Tuple{E, T}, ::DefaultAlgorithm) where {E <: Number, T <: Diagonal} = MAK.exponential!(t, DiagonalAlgorithm())
-MAK.exponential!(t::Tuple{E, T1}, out::T2, ::DefaultAlgorithm) where {E <: Number, T1 <: Diagonal, T2 <: Diagonal} = MAK.exponential!(t, out, DiagonalAlgorithm())
 
 function MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T}
     return MAK.default_algorithm(exponential!, blocktype(T); kwargs...)

From 68aff913420a1053d3daae20386d729c1f1a1d8c Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 10 Jul 2026 10:20:08 -0400
Subject: [PATCH 09/14] bump MAK version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 10d69842f..505e078ee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -54,7 +54,7 @@ FiniteDifferences = "0.12"
 GPUArrays = "11.4.1"
 LRUCache = "1.6"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.6.8"
+MatrixAlgebraKit = "0.6.9"
 Mooncake = "0.5.27"
 OhMyThreads = "0.8.0"
 Printf = "1"

From cf07354a3212f6be818d4e13c1903712afbc424b Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 10 Jul 2026 14:19:20 -0400
Subject: [PATCH 10/14] update implementations

---
 docs/src/man/linearalgebra.md          |  2 +-
 src/TensorKit.jl                       |  2 +-
 src/factorizations/diagonal.jl         |  5 ++++-
 src/factorizations/matrixalgebrakit.jl | 23 ++++++++++++-----------
 src/tensors/linalg.jl                  |  2 +-
 test/tensors/exponential.jl            |  4 ++--
 6 files changed, 21 insertions(+), 17 deletions(-)

diff --git a/docs/src/man/linearalgebra.md b/docs/src/man/linearalgebra.md
index cd545d65a..92f895ce3 100644
--- a/docs/src/man/linearalgebra.md
+++ b/docs/src/man/linearalgebra.md
@@ -40,7 +40,7 @@ Note that, because the adjoint interchanges domain and codomain, we have `space(
 
 When tensor map instances are endomorphisms, i.e. they have the same domain and codomain, there is a multiplicative identity which can be obtained as `one(t)` or `one!(t)`, where the latter overwrites the contents of `t`.
 The multiplicative identity on a space `V` can also be obtained using `id(A, V)` as discussed [above](@ref ss_tensor_construction), such that for a general homomorphism `t′`, we have `t′ == id(codomain(t′)) * t′ == t′ * id(domain(t′))`.
-Returning to the case of endomorphisms `t`, we can compute the trace via `tr(t)` and exponentiate them using `exp(t)`, or if the contents of `t` can be destroyed in the process, `exp!(t)`.
+Returning to the case of endomorphisms `t`, we can compute the trace via `tr(t)` and exponentiate them using `exp(t)` (or the more general `exponential(t)`), or if the contents of `t` can be destroyed in the process, `exponential!(t)`.
 Furthermore, there are a number of tensor factorizations for both endomorphisms and general homomorphisms that we discuss on the [Tensor factorizations](@ref ss_tensor_factorization) page.
 
 Finally, there are a number of operations that also belong in this paragraph because of their analogy to common matrix operations.
diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index 4d197c901..ac52c7641 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -82,7 +82,7 @@ export left_orth, right_orth, left_null, right_null,
     qr_full, qr_compact, qr_null, lq_full, lq_compact, lq_null,
     qr_full!, qr_compact!, qr_null!, lq_full!, lq_compact!, lq_null!,
     svd_compact!, svd_full!, svd_trunc!, svd_compact, svd_full, svd_trunc, svd_vals, svd_vals!,
-    exp, exp!,
+    exp, exp!, exponential, exponential!,
     eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eigh_vals!, eigh_vals,
     eig_full!, eig_full, eig_trunc!, eig_trunc, eig_vals!, eig_vals,
     eigen, eigen!,
diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl
index dae550ea1..dfbf22407 100644
--- a/src/factorizations/diagonal.jl
+++ b/src/factorizations/diagonal.jl
@@ -8,11 +8,14 @@ MAK.diagview(t::DiagonalTensorMap) = SectorVector(t.data, TensorKit.diagonalbloc
 for f in (
         :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null,
         :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full,
-        :eigh_trunc, :eigh_vals, :left_polar, :right_polar,
+        :eigh_trunc, :eigh_vals, :left_polar, :right_polar, :exponential,
     )
     @eval MAK.copy_input(::typeof($f), d::DiagonalTensorMap) = copy(d)
 end
 
+MAK.copy_input(::typeof(exponential), (τ, d)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} =
+    (τ, typeof(τ) <: Complex ? complex(d) : d)
+
 for f! in (:qr_full!, :qr_compact!)
     @eval function MAK.initialize_output(
             ::typeof($f!), d::AbstractTensorMap, ::DiagonalAlgorithm
diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index a1aa77c62..2a5e5c9b3 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -19,6 +19,11 @@ for f in
     end
 end
 
+MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T <: AbstractTensorMap} =
+    MAK.default_algorithm(exponential!, blocktype(T); kwargs...)
+MAK.copy_input(::typeof(exponential), (τ, t)::Tuple{E, T}) where {E <: Number, T <: AbstractTensorMap} =
+    (τ, copy_oftype(t, factorisation_scalartype(exponential, t)))
+
 _select_truncation(f, ::AbstractTensorMap, trunc::TruncationStrategy) = trunc
 function _select_truncation(::typeof(left_null!), ::AbstractTensorMap, trunc::NamedTuple)
     return MAK.null_truncation_strategy(; trunc...)
@@ -53,7 +58,7 @@ for f! in (
         :exponential!,
     )
     @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
-        $(f! in (:eig_vals!, :eigh_vals!, :project_hermitian!, :project_antihermitian!) && :(LinearAlgebra.checksquare(t)))
+        $(f! in (:eig_vals!, :eigh_vals!, :project_hermitian!, :project_antihermitian!, :exponential!) && :(LinearAlgebra.checksquare(t)))
         foreachblock(t, N) do _, (tblock, Nblock)
             Nblock′ = $f!(tblock, Nblock, alg)
             # deal with the case where the output is not the same as the input
@@ -76,16 +81,6 @@ function MAK.exponential!((τ, t)::Tuple{E, T}, N, alg::AbstractAlgorithm) where
     return N
 end
 
-# Default algorithm for exponential with Tuple
-MAK.exponential!((τ, t)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} = MAK.exponential!((τ, t), DefaultAlgorithm())
-
-function MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T}
-    return MAK.default_algorithm(exponential!, blocktype(T); kwargs...)
-end
-
-function MAK.copy_input(::typeof(exponential), (τ, t)::Tuple{E, T}) where {E <: Number, T <: AbstractTensorMap}
-    return (τ, copy_oftype(t, factorisation_scalartype(exponential, t)))
-end
 
 MAK.zero!(t::AbstractTensorMap) = zerovector!(t)
 
@@ -121,6 +116,12 @@ for f in [
         MAK.$f!(t, out, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
 end
 
+# resolve `DefaultAlgorithm` for the tuple form at the tensor level, mirroring the loop below
+MAK.exponential!((τ, t)::Tuple{E, T}, alg::DefaultAlgorithm) where {E <: Number, T <: AbstractTensorMap} =
+    MAK.exponential!((τ, t), MAK.select_algorithm(exponential!, t, nothing; alg.kwargs...))
+MAK.exponential!((τ, t)::Tuple{E, T}, out, alg::DefaultAlgorithm) where {E <: Number, T <: AbstractTensorMap} =
+    MAK.exponential!((τ, t), out, MAK.select_algorithm(exponential!, t, nothing; alg.kwargs...))
+
 
 # Singular value decomposition
 # ----------------------------
diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl
index 0852d9c27..4e5e4cce2 100644
--- a/src/tensors/linalg.jl
+++ b/src/tensors/linalg.jl
@@ -41,7 +41,7 @@ function compose(A::AbstractTensorMap, B::AbstractTensorMap)
 end
 Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) = compose(t1, t2)
 
-Base.exp(t::AbstractTensorMap) = exp!(copy(t))
+Base.exp(t::AbstractTensorMap) = exponential(t)
 function Base.:^(t::AbstractTensorMap, p::Integer)
     return p < 0 ? Base.power_by_squaring(inv(t), -p) : Base.power_by_squaring(t, p)
 end
diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index fbb2d73e9..007ce374f 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -2,8 +2,8 @@ using Test, TestExtras
 using TensorKit, MatrixAlgebraKit
 using Random
 
-spaces = [ℂ^4] #, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
-scalartypes = [Float64, ComplexF32] #, ComplexF64]
+spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
+scalartypes = [Float64, ComplexF32, ComplexF64]
 algorithms = [MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()), MatrixFunctionViaEigh(DefaultAlgorithm())]
 
 @timedtestset "exponential for Hermitian matrices with $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes

From bcb99cf91d681e3b9df5c6d3df3bd9fabc4672e1 Mon Sep 17 00:00:00 2001
From: sanderdemeyer <sander.demeyer@ugent.be>
Date: Sun, 12 Jul 2026 17:39:12 +0200
Subject: [PATCH 11/14] remove double test

---
 test/tensors/exponential.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index 007ce374f..65aff8c51 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -39,7 +39,6 @@ end
     expτA = exponential!((τ, copy(A)))
     expminτA = exponential!((-τ, copy(A)))
     @test expτA * expminτA ≈ id(scalartype(expτA), space)
-    @test expτA ≈ inv(expminτA)
 end
 
 @timedtestset "exponential! for diagonal matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
@@ -60,5 +59,4 @@ end
     expτA = exponential!((τ, copy(A)))
     expminτA = exponential!((-τ, copy(A)))
     @test expτA * expminτA ≈ id(scalartype(expτA), space)
-    @test expτA ≈ inv(expminτA)
 end

From bdb3e3144a0ff1684e6b9826453cdd3ca7a6f73a Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Mon, 13 Jul 2026 11:32:17 -0400
Subject: [PATCH 12/14] refactor exponential tests in the style of other tests

---
 test/tensors/exponential.jl | 120 ++++++++++++++++++++++--------------
 1 file changed, 74 insertions(+), 46 deletions(-)

diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index 65aff8c51..9a034e29b 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -1,62 +1,90 @@
 using Test, TestExtras
-using TensorKit, MatrixAlgebraKit
+using TensorKit
+using MatrixAlgebraKit: DefaultAlgorithm, MatrixFunctionViaLA, MatrixFunctionViaEig,
+    MatrixFunctionViaEigh, MatrixFunctionViaTaylor
 using Random
 
-spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
-scalartypes = [Float64, ComplexF32, ComplexF64]
-algorithms = [MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()), MatrixFunctionViaEigh(DefaultAlgorithm())]
+spacelist = default_spacelist(fast_tests)
+scalartypes = (Float32, Float64, ComplexF64)
 
-@timedtestset "exponential for Hermitian matrices with $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
-    A = randn(st1, space, space)
-    A = project_hermitian!(A)
-    τ = rand(st2)
+# algorithms that agree on Hermitian input
+hermitian_algs = (
+    MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()),
+    MatrixFunctionViaEigh(DefaultAlgorithm()), MatrixFunctionViaTaylor(),
+)
+# algorithms valid for general (non-Hermitian) input
+general_algs = (
+    MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()),
+    MatrixFunctionViaTaylor(),
+)
 
-    expA = @constinferred exponential(A)
-    expτA = @constinferred exponential((τ, A))
+for V in spacelist
+    I = sectortype(first(V))
+    Istr = TensorKit.type_repr(I)
+    println("---------------------------------------")
+    println("Matrix exponentials with symmetry: $Istr")
+    println("---------------------------------------")
+    @timedtestset "Matrix exponentials with symmetry: $Istr" verbose = true begin
+        V1, V2, V3, V4, V5 = V
+        W = V1 ⊗ V2 ⊗ V3
+        Vd = fuse(V1 ⊗ V2)
 
-    for alg in algorithms
-        expA2 = @constinferred exponential(A, alg)
-        expτA2 = @constinferred exponential((τ, A), alg)
+        # explicit (non-diagonal) algorithms only apply to dense endomorphisms;
+        # `DiagonalTensorMap` inputs are exercised through the default algorithm below
+        @testset "exponential for Hermitian matrices" begin
+            for T1 in scalartypes, T2 in scalartypes,
+                    A in (randn(T1, V1, V1), randn(T1, W, W))
 
-        @test expA ≈ expA2
-        @test expτA ≈ expτA2
-    end
-end
+                A = project_hermitian!(A)
+                τ = rand(T2)
 
-@timedtestset "exponential! for general matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
-    A = randn(st1, space, space)
-    τ = rand(st2)
+                expA = @constinferred exponential(A)
+                expτA = @constinferred exponential((τ, A))
 
-    @test exponential!(copy(A)) == exponential!((1.0, copy(A)))
+                for alg in hermitian_algs
+                    @test expA ≈ @constinferred exponential(A, alg)
+                    @test expτA ≈ @constinferred exponential((τ, A), alg)
+                end
+            end
+        end
 
-    A2 = exponential!((τ, A))
-    if st1 <: Real && st2 <: Complex
-        @test objectid(A2) != objectid(A)
-    else
-        @test objectid(A2) == objectid(A)
-    end
+        @testset "exponential! for general matrices" begin
+            for T1 in scalartypes, T2 in scalartypes,
+                    A in (
+                        randn(T1, V1, V1), randn(T1, W, W),
+                        DiagonalTensorMap(randn(T1, reduceddim(Vd)), Vd),
+                    )
 
-    expτA = exponential!((τ, copy(A)))
-    expminτA = exponential!((-τ, copy(A)))
-    @test expτA * expminτA ≈ id(scalartype(expτA), space)
-end
+                τ = rand(T2)
 
-@timedtestset "exponential! for diagonal matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
-    A = DiagonalTensorMap(rand(st1, reduceddim(space)), space)
-    τ = rand(st2)
+                expA = @constinferred exponential(A)
+                if !(A isa DiagonalTensorMap) # diagonal only supports the default algorithm == DiagonalAlgorithm
+                    for alg in general_algs
+                        @test expA ≈ exponential(A, alg)
+                    end
+                end
 
-    exponential!(copy(A))
-    @test exponential!(copy(A)) ≈ exponential!((1.0, copy(A))) #, DiagonalAlgorithm())
+                @test exponential!(copy(A)) ≈ exponential!((1.0, copy(A)))
 
-    A2 = @constinferred exponential!((τ, A))
-    @test A2 isa DiagonalTensorMap
-    if st1 <: Real && st2 <: Complex
-        @test objectid(A2) != objectid(A)
-    else
-        @test objectid(A2) == objectid(A)
-    end
+                # exp(A)² == exp(2A)
+                @test expA * expA ≈ exponential((2, A))
+
+                # in-place semantics: aliases the input, unless a complex scalar forces
+                # a real tensor to widen to complex
+                Ain = copy(A)
+                A2 = @constinferred exponential!((τ, Ain))
+                A isa DiagonalTensorMap && @test A2 isa DiagonalTensorMap
+                if T1 <: Real && T2 <: Complex
+                    @test A2 !== Ain
+                else
+                    @test A2 === Ain
+                end
 
-    expτA = exponential!((τ, copy(A)))
-    expminτA = exponential!((-τ, copy(A)))
-    @test expτA * expminτA ≈ id(scalartype(expτA), space)
+                # exp(τA) and exp(-τA) are inverse
+                expτA = exponential!((τ, copy(A)))
+                expmτA = exponential!((-τ, copy(A)))
+                @test expτA * expmτA ≈ id(scalartype(expτA), domain(A))
+            end
+        end
+    end
 end

From 647aaf07e7b1164b6e0499f442e2d6386e2fd108 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Mon, 13 Jul 2026 13:56:00 -0400
Subject: [PATCH 13/14] fix `copy_input` shouldn't alias

---
 src/factorizations/diagonal.jl         | 5 +----
 src/factorizations/matrixalgebrakit.jl | 2 +-
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl
index dfbf22407..dae550ea1 100644
--- a/src/factorizations/diagonal.jl
+++ b/src/factorizations/diagonal.jl
@@ -8,14 +8,11 @@ MAK.diagview(t::DiagonalTensorMap) = SectorVector(t.data, TensorKit.diagonalbloc
 for f in (
         :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null,
         :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full,
-        :eigh_trunc, :eigh_vals, :left_polar, :right_polar, :exponential,
+        :eigh_trunc, :eigh_vals, :left_polar, :right_polar,
     )
     @eval MAK.copy_input(::typeof($f), d::DiagonalTensorMap) = copy(d)
 end
 
-MAK.copy_input(::typeof(exponential), (τ, d)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} =
-    (τ, typeof(τ) <: Complex ? complex(d) : d)
-
 for f! in (:qr_full!, :qr_compact!)
     @eval function MAK.initialize_output(
             ::typeof($f!), d::AbstractTensorMap, ::DiagonalAlgorithm
diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index 2a5e5c9b3..9bb8531b7 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -22,7 +22,7 @@ end
 MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T <: AbstractTensorMap} =
     MAK.default_algorithm(exponential!, blocktype(T); kwargs...)
 MAK.copy_input(::typeof(exponential), (τ, t)::Tuple{E, T}) where {E <: Number, T <: AbstractTensorMap} =
-    (τ, copy_oftype(t, factorisation_scalartype(exponential, t)))
+    (τ, copy_oftype(t, (E <: Complex ? complex : identity)(factorisation_scalartype(exponential, t))))
 
 _select_truncation(f, ::AbstractTensorMap, trunc::TruncationStrategy) = trunc
 function _select_truncation(::typeof(left_null!), ::AbstractTensorMap, trunc::NamedTuple)

From 8712b744925dbc321f7e9ba08897b9576a4b143a Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Mon, 13 Jul 2026 17:00:30 -0400
Subject: [PATCH 14/14] some more small fixes

---
 src/factorizations/utility.jl | 2 +-
 test/tensors/exponential.jl   | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index 4d5e1bb00..0653cd7c6 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -5,7 +5,7 @@ end
 factorisation_scalartype(f, t) = factorisation_scalartype(t)
 
 function copy_oftype(t::AbstractTensorMap, T::Type{<:Number})
-    return copy!(similar(t, T, space(t)), t)
+    return copy!(similar(t, T), t)
 end
 
 function _reverse!(t::AbstractTensorMap; dims = :)
diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl
index 9a034e29b..16a6f7cca 100644
--- a/test/tensors/exponential.jl
+++ b/test/tensors/exponential.jl
@@ -81,9 +81,10 @@ for V in spacelist
                 end
 
                 # exp(τA) and exp(-τA) are inverse
+                # the inverse roundtrip is ill-conditioned; only assert at full (Float64) precision
                 expτA = exponential!((τ, copy(A)))
                 expmτA = exponential!((-τ, copy(A)))
-                @test expτA * expmτA ≈ id(scalartype(expτA), domain(A))
+                real(scalartype(expτA)) == Float64 && @test expτA * expmτA ≈ id(scalartype(expτA), domain(A))
             end
         end
     end