diff --git a/docs/src/man/precompilation.md b/docs/src/man/precompilation.md
index a9ab9dbd..a11fcf20 100644
--- a/docs/src/man/precompilation.md
+++ b/docs/src/man/precompilation.md
@@ -17,18 +17,19 @@ faster precompile times for fast TTFX for a wider range of inputs.
 
 ## Defaults
 
-By default, precompilation is disabled, but can be enabled for "tensors" of type `Array{T,N}`, where `T` and `N` range over the following values:
+By default, precompilation is enabled for "tensors" of type `Array{T,N}`, where `T` and `N` range over the following values:
 
 * `T` is either `Float64` or `ComplexF64`
 * `tensoradd!` is precompiled up to `N = 5`
 * `tensortrace!` is precompiled up to `4` free output indices and `2` pairs of traced indices
 * `tensorcontract!` is precompiled up to `3` free output indices on both inputs, and `2` contracted indices
 
-To enable precompilation with these default settings, you can *locally* change the `"precompile_workload"` key in the preferences.
+To disable precompilation altogether, for example during development or when you prefer to have small binaries,
+you can *locally* change the `"precompile_workload"` key in the preferences.
 
 ```julia
 using TensorOperations, Preferences
-set_preferences!(TensorOperations, "precompile_workload" => true; force=true)
+set_preferences!(TensorOperations, "precompile_workload" => false; force=true)
 ```
 
 ## Custom settings
@@ -43,12 +44,31 @@ set_preferences!(TensorOperations, "setting" => value; force=true)
 
 Here **setting** and **value** can take on the following:
 
-* `"precomple_eltypes"`: a `Vector{String}` that evaluate to the desired values of `T<:Number`
+* `"precompile_eltypes"`: a `Vector{String}` that evaluate to the desired values of `T<:Number`
 * `"precompile_add_ndims"`: an `Int` to specify the maximum `N` for `tensoradd!`
 * `"precompile_trace_ndims"`: a `Vector{Int}` of length 2 to specify the maximal number of free and traced indices for `tensortrace!`.
 * `"precompile_contract_ndims"`: a `Vector{Int}` of length 2 to specify the maximal number of free and contracted indices for `tensorcontract!`.
 
-!!! note "Backends"
+## Reuse in downstream packages
 
-    Currently, there is no support for precompiling methods that do not use the default backend. If this is a
-    feature you would find useful, feel free to contact us or open an issue.
+The default workload is factored into three reusable functions, one per operation family, which
+a downstream package can call from its own `PrecompileTools.@compile_workload` to precompile for a
+different backend, allocator, or array type:
+
+* `TensorOperations.precompile_tensoradd(T, N, backend, allocator)`
+* `TensorOperations.precompile_tensortrace(T, (N1, N2), backend, allocator)`
+* `TensorOperations.precompile_tensorcontract(T, (N1, N2, N3), backend, allocator)`
+
+Each takes a single scalar type `T`, a single index-rank specification, and (optionally) a
+`backend` and `allocator` selecting the implementation to precompile for. The tensors are built
+by `TensorOperations.precompile_maketensor(T, N)`, which returns a `Array{T,N}`. For example, to
+precompile the `tensoradd!` path for a custom backend and allocator:
+
+```julia
+using PrecompileTools, TensorOperations
+@compile_workload begin
+    for T in (Float64, ComplexF64), N in 0:3
+        TensorOperations.precompile_tensoradd(T, N, MyBackend(), MyAllocator())
+    end
+end
+```
diff --git a/src/implementation/abstractarray.jl b/src/implementation/abstractarray.jl
index 017b16f3..c9afb84e 100644
--- a/src/implementation/abstractarray.jl
+++ b/src/implementation/abstractarray.jl
@@ -63,9 +63,11 @@ _stridedordiag(A::Diagonal) = A
 
 Check that `C` has `numind(pC)` indices and that `pC` constitutes a valid permutation.
 """
-function argcheck_index2tuple(C::AbstractArray, pC::Index2Tuple)
-    return ndims(C) == numind(pC) && isperm(linearize(pC)) ||
+argcheck_index2tuple(C::AbstractArray, pC::Index2Tuple) = argcheck_indextuple(C, linearize(pC))
+function argcheck_indextuple(C::AbstractArray, pC::IndexTuple)
+    ndims(C) == numind(pC) && isperm(pC) ||
         throw(IndexError(lazy"invalid permutation of length $(ndims(C)): $pC"))
+    return nothing
 end
 
 """
@@ -73,9 +75,11 @@ end
 
 Check that `C` and `A` have `numind(pA)` indices and that `pA` constitutes a valid permutation.
 """
-function argcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::Index2Tuple)
+argcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::Index2Tuple) =
+    argcheck_tensoradd(C, A, linearize(pA))
+function argcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::IndexTuple)
     ndims(C) == ndims(A) || throw(IndexError("non-matching number of dimensions"))
-    argcheck_index2tuple(A, pA)
+    argcheck_indextuple(A, pA)
     return nothing
 end
 
@@ -85,14 +89,16 @@ end
 Check that the partial trace of `A` over indices `q` and with permutation of the remaining
 indices `p` is compatible with output `C`.
 """
+argcheck_tensortrace(C::AbstractArray, A::AbstractArray, p::Index2Tuple, q::Index2Tuple) =
+    argcheck_tensortrace(C, A, linearize(p), q)
 function argcheck_tensortrace(
-        C::AbstractArray, A::AbstractArray, p::Index2Tuple, q::Index2Tuple
+        C::AbstractArray, A::AbstractArray, p::IndexTuple, q::Index2Tuple
     )
     ndims(C) == numind(p) ||
         throw(IndexError(lazy"invalid selection of length $(ndims(C)): $p"))
     2 * numin(q) == 2 * numout(q) == ndims(A) - ndims(C) ||
         throw(IndexError("invalid number of trace dimensions"))
-    argcheck_index2tuple(A, ((p[1]..., q[1]...), (p[2]..., q[2]...)))
+    argcheck_indextuple(A, (p..., q[1]..., q[2]...))
     return nothing
 end
 
@@ -108,7 +114,15 @@ function argcheck_tensorcontract(
         B::AbstractArray, pB::Index2Tuple,
         pAB::Index2Tuple
     )
-    argcheck_index2tuple(C, pAB)
+    return argcheck_tensorcontract(C, A, pA, B, pB, linearize(pAB))
+end
+function argcheck_tensorcontract(
+        C::AbstractArray,
+        A::AbstractArray, pA::Index2Tuple,
+        B::AbstractArray, pB::Index2Tuple,
+        pAB::IndexTuple
+    )
+    argcheck_indextuple(C, pAB)
     argcheck_index2tuple(A, pA)
     argcheck_index2tuple(B, pB)
     numout(pA) + numin(pB) == ndims(C) ||
@@ -123,27 +137,28 @@ end
 
 Check that `C` and `A` have compatible sizes for the addition specified by `pA`.
 """
-function dimcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::Index2Tuple)
+dimcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::Index2Tuple) = dimcheck_tensoradd(C, A, linearize(pA))
+function dimcheck_tensoradd(C::AbstractArray, A::AbstractArray, pA::IndexTuple)
     szA, szC = size(A), size(C)
-    TupleTools.getindices(szA, linearize(pA)) == szC ||
+    TupleTools.getindices(szA, pA) == szC ||
         throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))
     return nothing
 end
 
 """
-    dimcheck_tensorcontract(C::AbstractArray, A::AbstractArray,
-                            p::Index2Tuple, q::Index2Tuple)
+    dimcheck_tensorcontract(C::AbstractArray, A::AbstractArray, p::Index2Tuple, q::Index2Tuple)
 
-Check that `C` and `A` have compatible sizes for the trace and addition specified by `p` and
-`q`.
+Check that `C` and `A` have compatible sizes for the trace and addition specified by `p` and `q`.
 """
+dimcheck_tensortrace(C::AbstractArray, A::AbstractArray, p::Index2Tuple, q::Index2Tuple) =
+    dimcheck_tensortrace(C, A, linearize(p), q)
 function dimcheck_tensortrace(
-        C::AbstractArray, A::AbstractArray, p::Index2Tuple, q::Index2Tuple
+        C::AbstractArray, A::AbstractArray, p::IndexTuple, q::Index2Tuple
     )
     szA, szC = size(A), size(C)
     TupleTools.getindices(szA, q[1]) == TupleTools.getindices(szA, q[2]) ||
         throw(DimensionMismatch("non-matching sizes in traced dimensions"))
-    TupleTools.getindices(szA, linearize(p)) == szC ||
+    TupleTools.getindices(szA, p) == szC ||
         throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))
     return nothing
 end
@@ -163,11 +178,19 @@ function dimcheck_tensorcontract(
         B::AbstractArray, pB::Index2Tuple,
         pAB::Index2Tuple
     )
+    return dimcheck_tensorcontract(C, A, pA, B, pB, linearize(pAB))
+end
+function dimcheck_tensorcontract(
+        C::AbstractArray,
+        A::AbstractArray, pA::Index2Tuple,
+        B::AbstractArray, pB::Index2Tuple,
+        pAB::IndexTuple
+    )
     szA, szB, szC = size(A), size(B), size(C)
     TupleTools.getindices(szA, pA[2]) == TupleTools.getindices(szB, pB[1]) ||
         throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
     szAB = (TupleTools.getindices(szA, pA[1])..., TupleTools.getindices(szB, pB[2])...)
-    TupleTools.getindices(szAB, linearize(pAB)) == szC ||
+    TupleTools.getindices(szAB, pAB) == szC ||
         throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))
     return nothing
 end
diff --git a/src/implementation/blascontract.jl b/src/implementation/blascontract.jl
index bc606260..72bb9278 100644
--- a/src/implementation/blascontract.jl
+++ b/src/implementation/blascontract.jl
@@ -6,14 +6,13 @@
 function blas_contract!(C, A, pA, B, pB, pAB, α, β, backend, allocator)
     rpA = reverse(pA)
     rpB = reverse(pB)
-    indCinoBA = let N₁ = numout(pA), N₂ = numin(pB)
-        map(n -> ifelse(n > N₁, n - N₁, n + N₂), linearize(pAB))
+    # note: `pAB` is only ever consumed through `linearize`/`invperm` from here on, so it is
+    # canonicalized to an `IndexTuple` and the reversed permutation does not need to be
+    # repartitioned the way `pAB` is
+    pAB = linearize(pAB)
+    rpAB = let N₁ = numout(pA), N₂ = numin(pB)
+        map(n -> ifelse(n > N₁, n - N₁, n + N₂), pAB)
     end
-    tpAB = trivialpermutation(pAB)
-    rpAB = (
-        TupleTools.getindices(indCinoBA, tpAB[1]),
-        TupleTools.getindices(indCinoBA, tpAB[2]),
-    )
     cp = allocator_checkpoint!(allocator)
     if contract_memcost(C, A, pA, B, pB, pAB) <= contract_memcost(C, B, rpB, A, rpA, rpAB)
         C = _blas_contract!(C, A, pA, B, pB, pAB, α, β, backend, allocator)
@@ -28,15 +27,15 @@ function blas_contract!(
         C::StridedView{T, 2},
         A::StridedView{T, 2}, pA::Index2Tuple{1, 1},
         B::StridedView{T, 2}, pB::Index2Tuple{1, 1},
-        pAB::Index2Tuple{1, 1},
+        pAB::IndexTuple{2},
         α::Number, β::Number,
         backend, allocator
     ) where {T}
     A′ = pA == ((1,), (2,)) ? A : transpose(A)
     B′ = pB == ((1,), (2,)) ? B : transpose(B)
-    if pAB == ((1,), (2,))
+    if pAB == (1, 2)
         mul!(C, A′, B′, α, β)
-    elseif pAB == ((2,), (1,))
+    elseif pAB == (2, 1)
         mul!(C, transpose(B′), transpose(A′), α, β)
     end
     return C
@@ -49,18 +48,23 @@ function _blas_contract!(C, A, pA, B, pB, pAB, α, β, backend, allocator)
     A_, pA, flagA = makeblascontractable(A, pA, TC, backend, allocator)
     B_, pB, flagB = makeblascontractable(B, pB, TC, backend, allocator)
 
+    # the partition of `ipAB` is required by `isblasdestination` and `tensoralloc_add`, but
+    # the actual contraction only needs the linearized permutation
     ipAB = oindABinC(pAB, pA, pB)
+    ipAB′ = linearize(ipAB)
     flagC = isblasdestination(C, ipAB)
     if flagC
         C_ = C
-        _unsafe_blas_contract!(C_, A_, pA, B_, pB, ipAB, α, β)
+        _unsafe_blas_contract!(C_, A_, pA, B_, pB, ipAB′, α, β)
     else
         C_ = SV(tensoralloc_add(TC, C, ipAB, false, Val(true), allocator))
         _unsafe_blas_contract!(
-            C_, A_, pA, B_, pB, trivialpermutation(ipAB),
+            C_, A_, pA, B_, pB, trivialpermutation(ipAB′),
             one(TC), zero(TC)
         )
-        tensoradd!(C, C_, pAB, false, α, β, backend, allocator)
+        # `C` already exists, so only the permutation of `pAB` matters here and it can be
+        # handed over in the canonical `{N,0}` partition
+        tensoradd!(C, C_, (pAB, ()), false, α, β, backend, allocator)
         tensorfree!(C_.parent, allocator)
     end
     flagA || tensorfree!(A_.parent, allocator)
@@ -74,7 +78,7 @@ function _unsafe_blas_contract!(
         C::StridedView{T},
         A::StridedView{T}, pA,
         B::StridedView{T}, pB,
-        pAB, α, β
+        pAB::IndexTuple, α, β
     ) where {T <: BlasFloat}
     sizeA = size(A)
     sizeB = size(B)
@@ -84,7 +88,7 @@ function _unsafe_blas_contract!(
     osizeB = TupleTools.getindices(sizeB, pB[2])
 
     mul!(
-        sreshape(permutedims(C, linearize(pAB)), (prod(osizeA), prod(osizeB))),
+        sreshape(permutedims(C, pAB), (prod(osizeA), prod(osizeB))),
         sreshape(permutedims(A, linearize(pA)), (prod(osizeA), prod(csizeA))),
         sreshape(permutedims(B, linearize(pB)), (prod(csizeB), prod(osizeB))),
         α, β
diff --git a/src/implementation/diagonal.jl b/src/implementation/diagonal.jl
index 5f085c4c..06a9fd25 100644
--- a/src/implementation/diagonal.jl
+++ b/src/implementation/diagonal.jl
@@ -9,17 +9,24 @@ function tensorcontract!(
         α::Number, β::Number,
         ::StridedNative, allocator = DefaultAllocator()
     )
-    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
-    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
+    @nospecialize allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
+
+    argcheck_tensorcontract(C, A, pA, B, pB, pAB′)
+    dimcheck_tensorcontract(C, A, pA, B, pB, pAB′)
 
     if conjA && conjB
-        _diagtensorcontract!(SV(C), conj(SV(A)), pA, conj(SV(B.diag)), pB, pAB, α, β)
+        _diagtensorcontract!(SV(C), conj(SV(A)), pA, conj(SV(B.diag)), pB, pAB′, α′, β′)
     elseif conjA
-        _diagtensorcontract!(SV(C), conj(SV(A)), pA, SV(B.diag), pB, pAB, α, β)
+        _diagtensorcontract!(SV(C), conj(SV(A)), pA, SV(B.diag), pB, pAB′, α′, β′)
     elseif conjB
-        _diagtensorcontract!(SV(C), SV(A), pA, conj(SV(B.diag)), pB, pAB, α, β)
+        _diagtensorcontract!(SV(C), SV(A), pA, conj(SV(B.diag)), pB, pAB′, α′, β′)
     else
-        _diagtensorcontract!(SV(C), SV(A), pA, SV(B.diag), pB, pAB, α, β)
+        _diagtensorcontract!(SV(C), SV(A), pA, SV(B.diag), pB, pAB′, α′, β′)
     end
     return C
 end
@@ -32,28 +39,32 @@ function tensorcontract!(
         α::Number, β::Number,
         ::StridedNative, allocator = DefaultAllocator()
     )
-    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
-    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
+    @nospecialize allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
+
+    argcheck_tensorcontract(C, A, pA, B, pB, pAB′)
+    dimcheck_tensorcontract(C, A, pA, B, pB, pAB′)
 
     rpA = reverse(pA)
     rpB = reverse(pB)
-    indCinoBA = let N₁ = numout(pA), N₂ = numin(pB)
-        map(n -> ifelse(n > N₁, n - N₁, n + N₂), linearize(pAB))
+    # note: `pAB` is only ever consumed through `linearize`/`invperm`, so the reversed
+    # permutation does not need to be repartitioned as `pAB` is
+    rpAB = let N₁ = numout(pA), N₂ = numin(pB)
+        map(n -> ifelse(n > N₁, n - N₁, n + N₂), pAB′)
     end
-    tpAB = trivialpermutation(pAB)
-    rpAB = (
-        TupleTools.getindices(indCinoBA, tpAB[1]),
-        TupleTools.getindices(indCinoBA, tpAB[2]),
-    )
 
     if conjA && conjB
-        _diagtensorcontract!(SV(C), conj(SV(B)), rpB, conj(SV(A.diag)), rpA, rpAB, α, β)
+        _diagtensorcontract!(SV(C), conj(SV(B)), rpB, conj(SV(A.diag)), rpA, rpAB, α′, β′)
     elseif conjA
-        _diagtensorcontract!(SV(C), SV(B), rpB, conj(SV(A.diag)), rpA, rpAB, α, β)
+        _diagtensorcontract!(SV(C), SV(B), rpB, conj(SV(A.diag)), rpA, rpAB, α′, β′)
     elseif conjB
-        _diagtensorcontract!(SV(C), conj(SV(B)), rpB, SV(A.diag), rpA, rpAB, α, β)
+        _diagtensorcontract!(SV(C), conj(SV(B)), rpB, SV(A.diag), rpA, rpAB, α′, β′)
     else
-        _diagtensorcontract!(SV(C), SV(B), rpB, SV(A.diag), rpA, rpAB, α, β)
+        _diagtensorcontract!(SV(C), SV(B), rpB, SV(A.diag), rpA, rpAB, α′, β′)
     end
     return C
 end
@@ -66,17 +77,24 @@ function tensorcontract!(
         α::Number, β::Number,
         ::StridedNative, allocator = DefaultAllocator()
     )
-    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
-    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
+    @nospecialize allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
+
+    argcheck_tensorcontract(C, A, pA, B, pB, pAB′)
+    dimcheck_tensorcontract(C, A, pA, B, pB, pAB′)
 
     if conjA && conjB
-        _diagdiagcontract!(SV(C), conj(SV(A.diag)), pA, conj(SV(B.diag)), pB, pAB, α, β)
+        _diagdiagcontract!(SV(C), conj(SV(A.diag)), pA, conj(SV(B.diag)), pB, pAB′, α′, β′)
     elseif conjA
-        _diagdiagcontract!(SV(C), conj(SV(A.diag)), pA, SV(B.diag), pB, pAB, α, β)
+        _diagdiagcontract!(SV(C), conj(SV(A.diag)), pA, SV(B.diag), pB, pAB′, α′, β′)
     elseif conjB
-        _diagdiagcontract!(SV(C), SV(A.diag), pA, conj(SV(B.diag)), pB, pAB, α, β)
+        _diagdiagcontract!(SV(C), SV(A.diag), pA, conj(SV(B.diag)), pB, pAB′, α′, β′)
     else
-        _diagdiagcontract!(SV(C), SV(A.diag), pA, SV(B.diag), pB, pAB, α, β)
+        _diagdiagcontract!(SV(C), SV(A.diag), pA, SV(B.diag), pB, pAB′, α′, β′)
     end
     return C
 end
@@ -89,21 +107,25 @@ function tensorcontract!(
         α::Number, β::Number,
         ::StridedNative, allocator = DefaultAllocator()
     )
-    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
-    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
+    @nospecialize allocator
 
-    A2 = StridedView(A.diag)
-    B2 = StridedView(B.diag)
-    C2 = StridedView(C.diag)
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
 
+    argcheck_tensorcontract(C, A, pA, B, pB, pAB′)
+    dimcheck_tensorcontract(C, A, pA, B, pB, pAB′)
+
+    C2 = SV(C.diag)
     if conjA && conjB
-        C2 .= C2 .* β .+ conj.(A2 .* B2) .* α
+        _diagdiagdiagcontract!(C2, conj(SV(A.diag)), conj(SV(B.diag)), α′, β′)
     elseif conjA
-        C2 .= C2 .* β .+ conj.(A2) .* B2 .* α
+        _diagdiagdiagcontract!(C2, conj(SV(A.diag)), SV(B.diag), α′, β′)
     elseif conjB
-        C2 .= C2 .* β .+ A2 .* conj.(B2) .* α
+        _diagdiagdiagcontract!(C2, SV(A.diag), conj(SV(B.diag)), α′, β′)
     else
-        C2 .= C2 .* β .+ A2 .* B2 .* α
+        _diagdiagdiagcontract!(C2, SV(A.diag), SV(B.diag), α′, β′)
     end
     return C
 end
@@ -112,7 +134,7 @@ function _diagtensorcontract!(
         C::StridedView,
         A::StridedView, pA::Index2Tuple,
         Bdiag::StridedView, pB::Index2Tuple,
-        pAB::Index2Tuple, α::Number, β::Number
+        pAB::IndexTuple, α::Number, β::Number
     )
     sizeA = i -> size(A, i)
     csizeA = sizeA.(pA[2])
@@ -122,7 +144,7 @@ function _diagtensorcontract!(
         totsize = (osizeA..., csizeA...)
         A2 = permutedims(A, linearize(pA))
         B2 = sreshape(Bdiag, ((one.(osizeA))..., csizeA...))
-        C2 = permutedims(C, invperm(linearize(pAB)))
+        C2 = permutedims(C, invperm(pAB))
 
     elseif numin(pB) == 0
         strideA = i -> stride(A, i)
@@ -130,7 +152,7 @@ function _diagtensorcontract!(
         totsize = (osizeA..., csizeA[1])
         A2 = StridedView(A.parent, totsize, newstrides, A.offset, A.op)
         B2 = sreshape(Bdiag, ((one.(osizeA))..., csizeA[1]))
-        C2 = permutedims(C, invperm(linearize(pAB)))
+        C2 = permutedims(C, invperm(pAB))
 
     else # numout(pB) == 2 # direct product
         scale!(C, β)
@@ -138,17 +160,28 @@ function _diagtensorcontract!(
         A2 = sreshape(permutedims(A, linearize(pA)), (osizeA..., 1))
         B2 = sreshape(Bdiag, ((one.(osizeA))..., length(Bdiag)))
 
-        C3 = permutedims(C, invperm(linearize(pAB)))
+        C3 = permutedims(C, invperm(pAB))
         sC = strides(C3)
         newstrides = (Base.front(Base.front(sC))..., sC[end - 1] + sC[end])
         totsize = (osizeA..., length(Bdiag))
         C2 = StridedView(C3.parent, totsize, newstrides, C3.offset, C3.op)
     end
 
-    op1 = Base.Fix2(scale, α) ∘ *
-    op2 = Base.Fix2(scale, β)
-    Strided._mapreducedim!(op1, +, op2, totsize, (C2, A2, B2))
+    Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), totsize, (C2, A2, B2))
+
+    return C
+end
 
+function _diagdiagdiagcontract!(
+        C::StridedView, Adiag::StridedView, Bdiag::StridedView, α::Number, β::Number
+    )
+    totsize = (length(C),)
+    # required: `β` was standardized, so `Zero()` no longer kills NaNs in uninitialized `C`
+    if iszero(β)
+        Strided._mapreducedim!(Scaler(α), nothing, nothing, totsize, (C, Adiag, Bdiag))
+    else
+        Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), totsize, (C, Adiag, Bdiag))
+    end
     return C
 end
 
@@ -156,7 +189,7 @@ function _diagdiagcontract!(
         C::StridedView,
         Adiag::StridedView, pA::Index2Tuple,
         Bdiag::StridedView, pB::Index2Tuple,
-        pAB::Index2Tuple, α::Number, β::Number
+        pAB::IndexTuple, α::Number, β::Number
     )
     if numin(pA) == 1 # matrix multiplication
         scale!(C, β)
@@ -181,16 +214,14 @@ function _diagdiagcontract!(
         A2 = sreshape(Adiag, (length(Adiag), 1))
         B2 = sreshape(Bdiag, (1, length(Adiag)))
 
-        C3 = permutedims(C, invperm(linearize(pAB)))
+        C3 = permutedims(C, invperm(pAB))
         strC = strides(C3)
         newstrides = (strC[1] + strC[2], strC[3] + strC[4])
         totsize = (length(A2), length(B2))
         C2 = StridedView(C3.parent, totsize, newstrides, C3.offset, C3.op)
     end
 
-    op1 = Base.Fix2(scale, α) ∘ *
-    op2 = Base.Fix2(scale, β)
-    Strided._mapreducedim!(op1, +, op2, totsize, (C2, A2, B2))
+    Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), totsize, (C2, A2, B2))
 
     return C
 end
diff --git a/src/implementation/strided.jl b/src/implementation/strided.jl
index c9669fd1..3dd4bf11 100644
--- a/src/implementation/strided.jl
+++ b/src/implementation/strided.jl
@@ -18,12 +18,20 @@ function tensoradd!(
         α::Number, β::Number,
         backend::StridedBackend, allocator = DefaultAllocator()
     )
+    @nospecialize backend allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    p = linearize(pA)
+
     # resolve conj flags and absorb into StridedView constructor to avoid type instabilities later on
     if conjA
-        stridedtensoradd!(SV(C), conj(SV(A)), pA, α, β, backend, allocator)
+        stridedtensoradd!(SV(C), conj(SV(A)), p, α′, β′)
     else
-        stridedtensoradd!(SV(C), SV(A), pA, α, β, backend, allocator)
+        stridedtensoradd!(SV(C), SV(A), p, α′, β′)
     end
+
     return C
 end
 
@@ -33,12 +41,20 @@ function tensortrace!(
         α::Number, β::Number,
         backend::StridedBackend, allocator = DefaultAllocator()
     )
+    @nospecialize backend allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    p′ = linearize(p)
+
     # resolve conj flags and absorb into StridedView constructor to avoid type instabilities later on
     if conjA
-        stridedtensortrace!(SV(C), conj(SV(A)), p, q, α, β, backend, allocator)
+        stridedtensortrace!(SV(C), conj(SV(A)), p′, q, α′, β′)
     else
-        stridedtensortrace!(SV(C), SV(A), p, q, α, β, backend, allocator)
+        stridedtensortrace!(SV(C), SV(A), p′, q, α′, β′)
     end
+
     return C
 end
 
@@ -48,26 +64,59 @@ function tensorcontract!(
         B::AbstractArray, pB::Index2Tuple, conjB::Bool,
         pAB::Index2Tuple,
         α::Number, β::Number,
-        backend::StridedBackend, allocator = DefaultAllocator()
+        backend::StridedBLAS, allocator = DefaultAllocator()
     )
+    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
+    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
+
+    (Base.mightalias(C, A) || Base.mightalias(C, B)) &&
+        throw(ArgumentError("output tensor must not be aliased with input tensor"))
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
+
     # resolve conj flags and absorb into StridedView constructor to avoid type instabilities later on
     if conjA && conjB
-        stridedtensorcontract!(
-            SV(C), conj(SV(A)), pA, conj(SV(B)), pB, pAB, α, β, backend, allocator
-        )
+        blas_contract!(SV(C), conj(SV(A)), pA, conj(SV(B)), pB, pAB′, α′, β′, backend, allocator)
     elseif conjA
-        stridedtensorcontract!(
-            SV(C), conj(SV(A)), pA, SV(B), pB, pAB, α, β, backend, allocator
-        )
+        blas_contract!(SV(C), conj(SV(A)), pA, SV(B), pB, pAB′, α′, β′, backend, allocator)
     elseif conjB
-        stridedtensorcontract!(
-            SV(C), SV(A), pA, conj(SV(B)), pB, pAB, α, β, backend, allocator
-        )
+        blas_contract!(SV(C), SV(A), pA, conj(SV(B)), pB, pAB′, α′, β′, backend, allocator)
     else
-        stridedtensorcontract!(
-            SV(C), SV(A), pA, SV(B), pB, pAB, α, β, backend, allocator
-        )
+        blas_contract!(SV(C), SV(A), pA, SV(B), pB, pAB′, α′, β′, backend, allocator)
     end
+
+    return C
+end
+
+function tensorcontract!(
+        C::AbstractArray,
+        A::AbstractArray, pA::Index2Tuple, conjA::Bool,
+        B::AbstractArray, pB::Index2Tuple, conjB::Bool,
+        pAB::Index2Tuple,
+        α::Number, β::Number,
+        backend::StridedNative, allocator = DefaultAllocator()
+    )
+    @nospecialize backend allocator
+
+    # standardize input types for compilation time
+    α′ = standardize_scalartype(C, α)
+    β′ = standardize_scalartype(C, β)
+    pAB′ = linearize(pAB)
+
+    # resolve conj flags and absorb into StridedView constructor to avoid type instabilities later on
+    if conjA && conjB
+        stridedtensorcontract!(SV(C), conj(SV(A)), pA, conj(SV(B)), pB, pAB′, α′, β′)
+    elseif conjA
+        stridedtensorcontract!(SV(C), conj(SV(A)), pA, SV(B), pB, pAB′, α′, β′)
+    elseif conjB
+        stridedtensorcontract!(SV(C), SV(A), pA, conj(SV(B)), pB, pAB′, α′, β′)
+    else
+        stridedtensorcontract!(SV(C), SV(A), pA, SV(B), pB, pAB′, α′, β′)
+    end
+
     return C
 end
 
@@ -83,40 +132,34 @@ end
 (s::Scaler)(x, y) = scale(x * y, s.α)
 
 function stridedtensoradd!(
-        C::StridedView,
-        A::StridedView, pA::Index2Tuple,
-        α::Number, β::Number,
-        ::StridedBackend, allocator = DefaultAllocator()
+        C::StridedView, A::StridedView, pA::IndexTuple, α::Number, β::Number,
     )
     argcheck_tensoradd(C, A, pA)
     dimcheck_tensoradd(C, A, pA)
-    if !istrivialpermutation(pA) && Base.mightalias(C, A)
+    !istrivialpermutation(pA) && Base.mightalias(C, A) &&
         throw(ArgumentError("output tensor must not be aliased with input tensor"))
+    Ap = permutedims(A, pA)
+    if iszero(β)
+        Strided._mapreducedim!(Scaler(α), nothing, nothing, size(C), (C, Ap))
+    else
+        Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), size(C), (C, Ap))
     end
-
-    A′ = permutedims(A, linearize(pA))
-    Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), size(C), (C, A′))
     return C
 end
 
 function stridedtensortrace!(
-        C::StridedView,
-        A::StridedView, p::Index2Tuple, q::Index2Tuple,
-        α::Number, β::Number,
-        ::StridedBackend, allocator = DefaultAllocator()
+        C::StridedView, A::StridedView, p::IndexTuple, q::Index2Tuple, α::Number, β::Number,
     )
     argcheck_tensortrace(C, A, p, q)
     dimcheck_tensortrace(C, A, p, q)
-
     Base.mightalias(C, A) &&
         throw(ArgumentError("output tensor must not be aliased with input tensor"))
-
-    sizeA = i -> size(A, i)
-    strideA = i -> stride(A, i)
-    tracesize = sizeA.(q[1])
-    newstrides = (strideA.(linearize(p))..., (strideA.(q[1]) .+ strideA.(q[2]))...)
-    newsize = (size(C)..., tracesize...)
-
+    newsize = linearize(size(C), TupleTools.getindices(size(A), q[1]))
+    stA = strides(A)
+    newstrides = linearize(
+        TupleTools.getindices(stA, p),
+        TupleTools.getindices(stA, q[1]) .+ TupleTools.getindices(stA, q[2])
+    )
     A′ = SV(A.parent, newsize, newstrides, A.offset, A.op)
     Strided._mapreducedim!(Scaler(α), Adder(), Scaler(β), newsize, (C, A′))
     return C
@@ -126,15 +169,10 @@ function stridedtensorcontract!(
         C::StridedView,
         A::StridedView, pA::Index2Tuple,
         B::StridedView, pB::Index2Tuple,
-        pAB::Index2Tuple,
+        pAB::IndexTuple,
         α::Number, β::Number,
-        backend::StridedBLAS, allocator = DefaultAllocator()
+        backend::StridedBLAS, allocator
     )
-    argcheck_tensorcontract(C, A, pA, B, pB, pAB)
-    dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
-
-    (Base.mightalias(C, A) || Base.mightalias(C, B)) &&
-        throw(ArgumentError("output tensor must not be aliased with input tensor"))
 
     blas_contract!(C, A, pA, B, pB, pAB, α, β, backend, allocator)
     return C
@@ -144,9 +182,8 @@ function stridedtensorcontract!(
         C::StridedView,
         A::StridedView, pA::Index2Tuple,
         B::StridedView, pB::Index2Tuple,
-        pAB::Index2Tuple,
-        α::Number, β::Number,
-        ::StridedNative, allocator = DefaultAllocator()
+        pAB::IndexTuple,
+        α::Number, β::Number
     )
     argcheck_tensorcontract(C, A, pA, B, pB, pAB)
     dimcheck_tensorcontract(C, A, pA, B, pB, pAB)
@@ -164,7 +201,7 @@ function stridedtensorcontract!(
         (one.(osizeA)..., osizeB..., csizeB...)
     )
     CS = sreshape(
-        permutedims(C, invperm(linearize(pAB))),
+        permutedims(C, invperm(pAB)),
         (osizeA..., osizeB..., one.(csizeA)...)
     )
     tsize = (osizeA..., osizeB..., csizeA...)
diff --git a/src/indices.jl b/src/indices.jl
index 5e5bfeb7..fdf8bf51 100644
--- a/src/indices.jl
+++ b/src/indices.jl
@@ -18,9 +18,12 @@ and `N₂` right indices.
 const Index2Tuple{N₁, N₂} = Tuple{IndexTuple{N₁}, IndexTuple{N₂}}
 
 linearize(p::Index2Tuple) = (p[1]..., p[2]...)
+linearize(p::IndexTuple) = p
+linearize(a::Tuple, b::Tuple) = (a..., b...)
 numout(p::Index2Tuple) = length(p[1])
 numin(p::Index2Tuple) = length(p[2])
 numind(p::Index2Tuple) = numout(p) + numin(p)
+numind(p::IndexTuple) = length(p)
 
 trivialpermutation(p::IndexTuple{N}) where {N} = ntuple(identity, Val(N))
 function trivialpermutation(p::Index2Tuple)
diff --git a/src/interface.jl b/src/interface.jl
index bdb5002f..85af1337 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -246,6 +246,23 @@ Obtain the type information of `C`, where `C` would be the output of
 """
 function tensorcontract_type end
 
+"""
+    standardize_scalartype(C, α::Number) -> α′
+
+Convert the scalar `α` into a standardized type suitable for combining with the output tensor `C`.
+This hook can be used to reduce the number of compiled specializations of some kernels,
+by ensuring that the scalars always arrive in a canonical type.
+
+!!! warning
+    Standardizing `Zero()` or `One()` discards their strong-zero and strong-one semantics, so a
+    kernel receiving a standardized scalar must either combine it through `VectorInterface`
+    (`scale`, `scale!`, `add`) or branch explicitly on `iszero(β)`. A plain `β * C` would
+    otherwise propagate `NaN`s from uninitialized output memory when `β == 0`.
+"""
+standardize_scalartype(C, α::Number) = _standardize_scalartype(scalartype(C), α)
+_standardize_scalartype(::Type{T}, α::Number) where {T <: Number} = convert(T, α)
+_standardize_scalartype(::Type, α::Number) = α      # ensure polynomials and symbolic types are left alone
+
 """
     tensoralloc(ttype, structure, [istemp=false, allocator])
 
diff --git a/src/precompile.jl b/src/precompile.jl
index 829e27e8..9ce13274 100644
--- a/src/precompile.jl
+++ b/src/precompile.jl
@@ -1,5 +1,5 @@
-using PrecompileTools: PrecompileTools
-using Preferences: @load_preference, load_preference
+using PrecompileTools: PrecompileTools, @setup_workload, @compile_workload
+using Preferences: @load_preference
 
 # Validate preferences input
 # --------------------------
@@ -49,92 +49,110 @@ const PRECOMPILE_CONTRACT_NDIMS = validate_contract_ndims(
     @load_preference("precompile_contract_ndims", [4, 2])
 )
 
-# Copy from PrecompileTools.workload_enabled but default to false
-function workload_enabled(mod::Module = @__MODULE__)
-    return try
-        if load_preference(PrecompileTools, "precompile_workloads", true)
-            return load_preference(mod, "precompile_workload", false)
-        else
-            return false
-        end
-    catch
-        false
-    end
+# Precompilation workload
+# ------------------------
+# The workload actually runs representative tensor operations so that PrecompileTools caches
+# the specializations that get compiled. Each operation family is factored into a reusable
+# `precompile_*` function that can also be called from a downstream package's own
+# `@compile_workload` to precompile for a different backend, allocator, or array type.
+#
+# `@compile_workload` is enabled by default and honors the standard `precompile_workload`
+# preference, which can be flipped to disable precompilation:
+#
+#     using TensorOperations, Preferences
+#     set_preferences!(TensorOperations, "precompile_workload" => false; force=true)
+
+# Tensor constructor used by the precompile workloads: build a rank-`N` tensor of scalar type
+# `T`. Downstream callers can add methods for other array types to reuse the `precompile_*`
+# functions below.
+precompile_maketensor(T, N) = zeros(T, ntuple(Returns(2), N))
+
+"""
+    precompile_tensoradd(T, N, backend, allocator)
+
+Run [`tensoradd!`](@ref) and [`tensoralloc_add`](@ref) for scalar type `T` and output rank `N`,
+using `backend` and `allocator`, so that their specializations are precompiled.
+"""
+function precompile_tensoradd(
+        T, N, backend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    C = precompile_maketensor(T, N)
+    A = precompile_maketensor(T, N)
+    pA = (ntuple(identity, N), ())
+
+    tensoradd!(C, A, pA, false, One(), Zero(), backend, allocator)
+    tensoradd!(C, A, pA, false, one(T), Zero(), backend, allocator)
+    tensoradd!(C, A, pA, false, one(T), zero(T), backend, allocator)
+
+    tensoralloc_add(T, A, pA, false, Val(true), allocator)
+    tensoralloc_add(T, A, pA, false, Val(false), allocator)
+    return nothing
 end
 
-# Using explicit precompile statements here instead of @compile_workload:
-# Actually running the precompilation through PrecompileTools leads to longer compile times
-# Keeping the workload_enabled functionality to have the option of disabling precompilation
-# in a compatible manner with the rest of the ecosystem
-if workload_enabled()
-    # tensoradd!
-    # ----------
-    for T in PRECOMPILE_ELTYPES
-        for N in 0:PRECOMPILE_ADD_NDIMS
-            C = Array{T, N}
-            A = Array{T, N}
-            pA = Index2Tuple{N, 0}
-
-            precompile(tensoradd!, (C, A, pA, Bool, One, Zero))
-            precompile(tensoradd!, (C, A, pA, Bool, T, Zero))
-            precompile(tensoradd!, (C, A, pA, Bool, T, T))
-
-            precompile(tensoralloc_add, (T, A, pA, Bool, Val{true}))
-            precompile(tensoralloc_add, (T, A, pA, Bool, Val{false}))
-        end
-    end
-
-    # tensortrace!
-    # ------------
-    for T in PRECOMPILE_ELTYPES
-        for N1 in 0:PRECOMPILE_TRACE_NDIMS[1], N2 in 0:PRECOMPILE_TRACE_NDIMS[2]
-            C = Array{T, N1}
-            A = Array{T, N1 + 2N2}
-            p = Index2Tuple{N1, 0}
-            q = Index2Tuple{N2, N2}
-
-            precompile(tensortrace!, (C, A, p, q, Bool, One, Zero))
-            precompile(tensortrace!, (C, A, p, q, Bool, T, Zero))
-            precompile(tensortrace!, (C, A, p, q, Bool, T, T))
+"""
+    precompile_tensortrace(T, (N1, N2), backend, allocator)
+
+Run [`tensortrace!`](@ref) for scalar type `T`, output rank `N1`, and `N2` traced index pairs,
+using `backend` and `allocator`, so that their specializations are precompiled.
+"""
+function precompile_tensortrace(
+        T, (N1, N2), backend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    C = precompile_maketensor(T, N1)
+    A = precompile_maketensor(T, N1 + 2N2)
+    p = (ntuple(identity, N1), ())
+    q = (ntuple(i -> N1 + i, N2), ntuple(i -> N1 + N2 + i, N2))
+
+    tensortrace!(C, A, p, q, false, One(), Zero(), backend, allocator)
+    tensortrace!(C, A, p, q, false, one(T), Zero(), backend, allocator)
+    tensortrace!(C, A, p, q, false, one(T), zero(T), backend, allocator)
+
+    # allocation re-uses tensoralloc_add
+    return nothing
+end
 
-            # allocation re-uses tensoralloc_add
-        end
-    end
+"""
+    precompile_tensorcontract(T, (N1, N2, N3), backend, allocator)
+
+Run [`tensorcontract!`](@ref) and [`tensoralloc_contract`](@ref) for scalar type `T`, with `N1`
+and `N3` free output indices on the two inputs and `N2` contracted indices, using `backend` and
+`allocator`, so that their specializations are precompiled.
+"""
+function precompile_tensorcontract(
+        T, (N1, N2, N3), backend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    NA = N1 + N2
+    NB = N2 + N3
+    NC = N1 + N3
+    C = precompile_maketensor(T, NC)
+    A = precompile_maketensor(T, NA)
+    B = precompile_maketensor(T, NB)
+    pA = (ntuple(identity, N1), ntuple(i -> N1 + i, N2))
+    pB = (ntuple(identity, N2), ntuple(i -> N2 + i, N3))
+    pAB = (ntuple(identity, NC), ())
+
+    tensorcontract!(C, A, pA, false, B, pB, false, pAB, One(), Zero(), backend, allocator)
+    tensorcontract!(C, A, pA, false, B, pB, false, pAB, one(T), Zero(), backend, allocator)
+    tensorcontract!(C, A, pA, false, B, pB, false, pAB, one(T), zero(T), backend, allocator)
+
+    tensoralloc_contract(T, A, pA, false, B, pB, false, pAB, Val(true), allocator)
+    tensoralloc_contract(T, A, pA, false, B, pB, false, pAB, Val(false), allocator)
+    return nothing
+end
 
-    # tensorcontract!
-    # ---------------
-    for T in PRECOMPILE_ELTYPES
-        for N1 in 0:PRECOMPILE_CONTRACT_NDIMS[1], N2 in 0:PRECOMPILE_CONTRACT_NDIMS[2],
-                N3 in 0:PRECOMPILE_CONTRACT_NDIMS[1]
-
-            NA = N1 + N2
-            NB = N2 + N3
-            NC = N1 + N3
-            C, A, B = Array{T, NC}, Array{T, NA}, Array{T, NB}
-            pA = Index2Tuple{N1, N2}
-            pB = Index2Tuple{N2, N3}
-            pAB = Index2Tuple{NC, 0}
-
-            precompile(tensorcontract!, (C, A, pA, Bool, B, pB, Bool, pAB, One, Zero))
-            precompile(tensorcontract!, (C, A, pA, Bool, B, pB, Bool, pAB, T, Zero))
-            precompile(tensorcontract!, (C, A, pA, Bool, B, pB, Bool, pAB, T, T))
-
-            precompile(tensoralloc_contract, (T, A, pA, Bool, B, pB, Bool, pAB, Val{true}))
-            precompile(tensoralloc_contract, (T, A, pA, Bool, B, pB, Bool, pAB, Val{false}))
+@setup_workload begin
+    @compile_workload begin
+        for T in PRECOMPILE_ELTYPES
+            for N in 0:PRECOMPILE_ADD_NDIMS
+                precompile_tensoradd(T, N)
+            end
+            for N1 in 0:PRECOMPILE_TRACE_NDIMS[1], N2 in 0:PRECOMPILE_TRACE_NDIMS[2]
+                precompile_tensortrace(T, (N1, N2))
+            end
+            for N1 in 0:PRECOMPILE_CONTRACT_NDIMS[1], N2 in 0:PRECOMPILE_CONTRACT_NDIMS[2],
+                    N3 in 0:PRECOMPILE_CONTRACT_NDIMS[1]
+                precompile_tensorcontract(T, (N1, N2, N3))
+            end
         end
     end
-else
-    @info """
-    TensorOperations can optionally be instructed to precompile several functions, which can be used to reduce the time to first execution (TTFX).
-    This is disabled by default as this can take a while on some machines, and is only relevant for contraction-heavy workloads.
-
-    To enable or disable precompilation, you can use the following script:
-
-    ```julia
-    using TensorOperations, Preferences
-    set_preferences!(TensorOperations, "precompile_workload" => true; force=true)
-    ```
-
-    This will create a `LocalPreferences.toml` file next to your current `Project.toml` file to store this setting in a persistent way.
-    """
 end
diff --git a/test/methods.jl b/test/methods.jl
index 22b8afec..501e9ace 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -252,3 +252,54 @@ end
         )
     end
 end
+
+# β = 0 must act as a strong zero: uninitialized output must not leak NaNs
+#-------------------------------------------------------------------------
+@testset "strong zero β with Diagonal ($T)" for T in (Float64, ComplexF64)
+    Zero = TensorOperations.VectorInterface.Zero
+    n = 4
+    d = randn(T, n)
+    S = Diagonal(d)
+    A = randn(T, n, n, n, n)
+    F = randn(T, n, n)
+    poison(dims...) = fill(T(NaN), dims...)
+
+    # C[i,j,k,l] = A[i,j,k,a] * S[a,l]
+    C1 = poison(n, n, n, n)
+    tensorcontract!(
+        C1, A, ((1, 2, 3), (4,)), false, S, ((1,), (2,)), false,
+        ((1, 2, 3, 4), ()), 1, Zero()
+    )
+    @test C1 ≈ [A[i, j, k, l] * d[l] for i in 1:n, j in 1:n, k in 1:n, l in 1:n]
+
+    # C[i,j] = A[i,j,a,b] * S[a,b]
+    C2 = poison(n, n)
+    tensorcontract!(
+        C2, A, ((1, 2), (3, 4)), false, S, ((1, 2), ()), false,
+        ((1, 2), ()), 1, Zero()
+    )
+    @test C2 ≈ [sum(A[i, j, a, a] * d[a] for a in 1:n) for i in 1:n, j in 1:n]
+
+    # C[a,c,b,d] = F[a,b] * S[c,d]
+    C3 = poison(n, n, n, n)
+    tensorcontract!(
+        C3, F, ((1, 2), ()), false, S, ((), (1, 2)), false,
+        ((1, 3, 2, 4), ()), 1, Zero()
+    )
+    @test C3 ≈ [F[a, b] * (c == e ? d[c] : zero(T)) for a in 1:n, c in 1:n, b in 1:n, e in 1:n]
+
+    # C[i,k] = S[i,j] * S[j,k], into a dense matrix and into a Diagonal
+    C4 = poison(n, n)
+    tensorcontract!(
+        C4, S, ((1,), (2,)), false, S, ((1,), (2,)), false,
+        ((1, 2), ()), 1, Zero()
+    )
+    @test C4 ≈ Array(S * S)
+
+    C5 = Diagonal(poison(n))
+    tensorcontract!(
+        C5, S, ((1,), (2,)), false, S, ((1,), (2,)), false,
+        ((1, 2), ()), 1, Zero()
+    )
+    @test C5 ≈ S * S
+end
diff --git a/test/tensor.jl b/test/tensor.jl
index 60b0b390..2ac3f71e 100644
--- a/test/tensor.jl
+++ b/test/tensor.jl
@@ -357,6 +357,10 @@ end
     S3 = similar(S)
     @tensor S3[i, k] = S[i, j] * S[j, k]
     @test S2 ≈ S3 ≈ S * S
+    # β = 0 must act as a strong zero: uninitialized output must not leak NaNs
+    S4 = Diagonal(fill(T(NaN), 5))
+    @tensor S4[i, k] = S[i, j] * S[j, k]
+    @test S4 ≈ S * S
     Str = @tensor S[i, j] * S[i, j]
     @test Str ≈ sum(S.diag .^ 2)