diff --git a/src/Grassmann.jl b/src/Grassmann.jl
index a62bfbe..1fb79cc 100644
--- a/src/Grassmann.jl
+++ b/src/Grassmann.jl
@@ -26,12 +26,12 @@ import AbstractTensors: Values, Variables, FixedVector, clifford, hodge, wedge,
 export ⊕, ℝ, @V_str, @S_str, @D_str, Manifold, Submanifold, Signature, DiagonalForm, value
 export @basis, @basis_str, @dualbasis, @dualbasis_str, @mixedbasis, @mixedbasis_str, Λ
 export ℝ0, ℝ1, ℝ2, ℝ3, ℝ4, ℝ5, ℝ6, ℝ7, ℝ8, ℝ9, mdims, tangent, metric, antimetric
-export hodge, wedge, vee, complement, dot, antidot
+export hodge, wedge, vee, complement, dot, antidot, istangent
 
 import Base: @pure, ==, isapprox
 import Base: print, show, getindex, setindex!, promote_rule, convert, adjoint
 import DirectSum: V0, ⊕, generate, basis, getalgebra, getbasis, dual, Zero, One, Zero, One
-import Leibniz: hasinf, hasorigin, dyadmode, value, pre, vsn, metric, mdims
+import Leibniz: hasinf, hasorigin, dyadmode, value, pre, vsn, metric, mdims, gdims
 import Leibniz: bit2int, indexbits, indices, diffvars, diffmask
 import Leibniz: symmetricmask, indexstring, indexsymbol, combo, digits_fast
 import DirectSum: antimetric
@@ -117,8 +117,8 @@ end
 d(ω::T) where T<:TensorAlgebra = Manifold(ω)(∇)∧ω
 δ(ω::T) where T<:TensorAlgebra = -∂(ω)
 
-function boundary_rank(t,d=gdims(t))
-    out = gdims(∂(t))
+function boundary_rank(t,d=count_gdims(t))
+    out = count_gdims(∂(t))
     out[1] = 0
     for k ∈ 2:length(out)-1
         @inbounds out[k] = min(out[k],d[k+1])
@@ -127,7 +127,7 @@ function boundary_rank(t,d=gdims(t))
 end
 
 function boundary_null(t)
-    d = gdims(t)
+    d = count_gdims(t)
     r = boundary_rank(t,d)
     l = length(d)
     out = zeros(Variables{l,Int})
@@ -143,7 +143,7 @@ end
 Compute the Betti numbers.
 """
 function betti(t::T) where T<:TensorAlgebra
-    d = gdims(t)
+    d = count_gdims(t)
     r = boundary_rank(t,d)
     l = length(d)-1
     out = zeros(Variables{l,Int})
@@ -569,8 +569,7 @@ function __init__()
         DyadicChain(m::StaticArrays.SMatrix{N,N}) where N = Chain{Submanifold(N),1}(m)
         Chain{V,G}(m::StaticArrays.SMatrix{N,N}) where {V,G,N} = Chain{V,G}(Chain{V,G}.(getindex.(Ref(m),:,StaticArrays.SVector{N}(1:N))))
         Chain{V,G,Chain{W,G}}(m::StaticArrays.SMatrix{M,N}) where {V,W,G,M,N} = Chain{V,G}(Chain{W,G}.(getindex.(Ref(m),:,StaticArrays.SVector{N}(1:N))))
-        Base.exp(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(exp(StaticArrays.SMatrix(A)))
-        Base.log(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(log(StaticArrays.SMatrix(A)))
+        #Base.log(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(log(StaticArrays.SMatrix(A)))
         LinearAlgebra.eigvals(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain(Values{binomial(mdims(V),G)}(LinearAlgebra.eigvals(StaticArrays.SMatrix(A))))
         LinearAlgebra.eigvecs(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain(Chain.(Values{binomial(mdims(A),G)}.(getindex.(Ref(LinearAlgebra.eigvecs(StaticArrays.SMatrix(A))),:,list(1,binomial(mdims(A),G))))))
         function LinearAlgebra.eigen(A::Chain{V,G,<:Chain{V,G}}) where {V,G}
diff --git a/src/composite.jl b/src/composite.jl
index dc81447..8e9a8bd 100644
--- a/src/composite.jl
+++ b/src/composite.jl
@@ -216,6 +216,109 @@ function Base.exp(t::T,::Val{hint}) where T<:TensorGraded{V} where {V,hint}
     end
 end
 
+@inline Base.expm1(A::DyadicChain) = exp(A)-I
+@inline Base.exp(A::DyadicChain{V,G,T,1}) where {V,G,T} = Chain{V,G}(Values(Chain{V,G}(exp(A[1][1]))))
+
+@inline function Base.exp(A::DyadicChain{V,G,<:Real,2}) where {V,G}
+    T = typeof(exp(zero(valuetype(A))))
+    @inbounds a = A[1][1]
+    @inbounds c = A[1][2]
+    @inbounds b = A[2][1]
+    @inbounds d = A[2][2]
+    v = (a-d)^2 + 4*b*c
+    if v > 0
+        z = sqrt(v)
+        z1 = cosh(z / 2)
+        z2 = sinh(z / 2) / z
+    elseif v < 0
+        z = sqrt(-v)
+        z1 = cos(z / 2)
+        z2 = sin(z / 2) / z
+    else # if v == 0
+        z1 = T(1.0)
+        z2 = T(0.5)
+    end
+    r = exp((a + d) / 2)
+    m11 = r * (z1 + (a - d) * z2)
+    m12 = r * 2b * z2
+    m21 = r * 2c * z2
+    m22 = r * (z1 - (a - d) * z2)
+    Chain{V,G}(Chain{V,G}(m11, m21), Chain{V,G}(m12, m22))
+end
+
+@inline function Base.exp(A::DyadicChain{V,G,<:Complex,2}) where {V,G}
+    T = typeof(exp(zero(valuetype(A))))
+    @inbounds a = A[1][1]
+    @inbounds c = A[1][2]
+    @inbounds b = A[2][1]
+    @inbounds d = A[2][2]
+    z = sqrt((a - d)*(a - d) + 4*b*c )
+    e = expm1((a + d - z) / 2)
+    f = expm1((a + d + z) / 2)
+    ϵ = eps()
+    g = abs2(z) < ϵ^2 ? exp((a + d) / 2) * (1 + z^2 / 24) : (f - e) / z
+    m11 = (g * (a - d) + f + e) / 2 + 1
+    m12 = g * b
+    m21 = g * c
+    m22 = (-g * (a - d) + f + e) / 2 + 1
+    Chain{V,G}(Chain{V,G}(m11, m21), Chain{V,G}(m12, m22))
+end
+
+# Adapted from implementation in Base; algorithm from
+# Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
+function Base.exp(_A::DyadicChain{W,G,T,N}) where {W,G,T,N}
+    S = typeof((zero(T)*zero(T) + zero(T)*zero(T))/one(T))
+    A = Chain{W,G}(map.(S,value(_A)))
+    # omitted: matrix balancing, i.e., LAPACK.gebal!
+    nA = maximum(sum.(value.(map.(abs,value(A)))))
+    # marginally more performant than norm(A, 1)
+    ## For sufficiently small nA, use lower order Padé-Approximations
+    if (nA <= 2.1)
+        A2 = A*A
+        if nA > 0.95
+            U = S(8821612800)*I+A2*(S(302702400)*I+A2*(S(2162160)*I+A2*(S(3960)*I+A2)))
+            U = A*U
+            V = S(17643225600)*I+A2*(S(2075673600)*I+A2*(S(30270240)*I+A2*(S(110880)*I+S(90)*A2)))
+        elseif nA > 0.25
+            U = S(8648640)*I+A2*(S(277200)*I+A2*(S(1512)*I+A2))
+            U = A*U
+            V = S(17297280)*I+A2*(S(1995840)*I+A2*(S(25200)*I+S(56)*A2))
+        elseif nA > 0.015
+            U = S(15120)*I+A2*(S(420)*I+A2)
+            U = A*U
+            V = S(30240)*I+A2*(S(3360)*I+S(30)*A2)
+        else
+            U = S(60)*I+A2
+            U = A*U
+            V = S(120)*I+S(12)*A2
+        end
+        expA = (V - U) \ (V + U)
+    else
+        s  = log2(nA/5.4)               # power of 2 later reversed by squaring
+        if s > 0
+            si = ceil(Int,s)
+            A = A / S(2^si)
+        end
+        A2 = A*A
+        A4 = A2*A2
+        A6 = A2*A4
+        U = A6*(A6 + S(16380)*A4 + S(40840800)*A2) +
+            (S(33522128640)*A6 + S(10559470521600)*A4 + S(1187353796428800)*A2) +
+            S(32382376266240000)*I
+        U = A*U
+        V = A6*(S(182)*A6 + S(960960)*A4 + S(1323241920)*A2) +
+            (S(670442572800)*A6 + S(129060195264000)*A4 + S(7771770303897600)*A2) +
+            S(64764752532480000)*I
+        expA = (V - U) \ (V + U)
+        if s > 0            # squaring to reverse dividing by power of 2
+            for t=1:si
+                expA = expA*expA
+            end
+        end
+    end
+    expA
+end
+
 qlog(b::PseudoCouple,x::Int=10000) = qlog(multispin(b),x)
 qlog(b::AntiSpinor,x::Int=10000) = qlog(Multivector(b),x)
 function qlog(w::T,x::Int=10000) where T<:TensorAlgebra
@@ -804,6 +907,7 @@ function inv_approx(t::Chain{M,1,<:Chain{V,1}}) where {M,V}
     mdims(M) < mdims(V) ? (inv(tt⋅t))⋅tt : tt⋅inv(t⋅tt)
 end
 
+Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1,<:Chain}) where {M,W,V} = Chain{V,1}(t.\value(v))
 Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1}) where {M,W,V} = value(t)\v
 Base.in(v::Chain{V,1},t::Chain{W,1,<:Chain{V,1}}) where {V,W} = v ∈ value(t)
 Base.inv(t::Chain{V,1,<:Chain{W,1}}) where {W,V} = inv(value(t))
diff --git a/src/multivectors.jl b/src/multivectors.jl
index 82c3050..7887595 100644
--- a/src/multivectors.jl
+++ b/src/multivectors.jl
@@ -37,6 +37,18 @@ end
 import Leibniz: Fields, parval, mixed, mvecs, svecs, spinsum, spinsum_set
 parsym = (Symbol,parval...)
 
+"""
+    compactio(::Bool)
+
+Toggles compact numbers for extended `Grassmann` elements (enabled by default)
+"""
+compact = ( () -> begin
+        io=true
+        return (tf=io)->(io≠tf && (io=tf); return io)
+    end)()
+
+compactio(io::IO) = compact() ? IOContext(io, :compact => true) : io
+
 function showterm(io::IO,V,B::UInt,i::T,compact=get(io,:compact,false)) where T
     if !(|(broadcast(<:,T,parsym)...)) && signbit(i) && !isnan(i)
         print(io, compact ? "-" : " - ")
@@ -84,17 +96,24 @@ Chain(v::Chain{V,G,𝕂}) where {V,G,𝕂} = v
 
 DyadicProduct{V,W,G,T,N} = Chain{V,G,Chain{W,G,T,N},N}
 DyadicChain{V,G,T,N} = DyadicProduct{V,V,G,T,N}
+#TriadicProduct{V,W,G,T,N} = Chain{V,G,DyadicChain{W,G,T,N},N}
+#DyadicChain{V,G,T,N} = Chain{V,G,Chain{V,G,T,N},N}
+#TriadicChain{V,G,T,N,M} = Chain{V,G,DyadicChain{V,1,T,M},N}
 
 Base.Matrix(m::Chain{V,G,<:TensorGraded{W,G}}) where {V,W,G} = hcat(value.(Chain.(value(m)))...)
 Base.Matrix(m::Chain{V,G,<:Chain{W,G}}) where {V,W,G} = hcat(value.(value(m))...)
-DyadicChain(m::Matrix) = Chain{Submanifold(size(m)[1]),1}(m)
+Chain{V}(m::Matrix) where V = Chain{V,1}(m)
 function Chain{V,G}(m::Matrix) where {V,G}
     N = size(m)[2]
     Chain{V,G,Chain{N≠mdims(V) ? Submanifold(N) : V,G}}(m)
 end
 Chain{V,G,Chain{W,G}}(m::Matrix) where {V,W,G} = Chain{V,G}(Chain{W,G}.(getindex.(Ref(m),:,list(1,size(m)[2]))))
+DyadicChain{V,G}(m::Matrix) where {V,G} = Chain{V,G}(Chain{V,G}.(getindex.(Ref(m),:,list(1,size(m)[2]))))
+DyadicChain{V}(m::Matrix) where V = DyadicChain{V,1}(m)
+DyadicChain(m::Matrix) = DyadicChain{Submanifold(size(m)[1])}(m)
+Base.log(A::DyadicChain{V,G}) where {V,G} = DyadicChain{V,G}(log(Matrix(A)))
 
-export Chain, DyadicProduct, DyadicChain
+export Chain, DyadicChain, TriadicChain, DyadicProduct
 getindex(m::Chain,i::Int) = m.v[i]
 getindex(m::Chain,i::UnitRange{Int}) = m.v[i]
 getindex(m::Chain,i::T) where T<:AbstractVector = m.v[i]
@@ -118,6 +137,7 @@ Base.transpose(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = _transpose(value(t),V)
 
 function show(io::IO, m::Chain{V,G,T}) where {V,G,T}
     ib,compact = indexbasis(mdims(V),G),get(io,:compact,false)
+    io = compactio(io)
     @inbounds Leibniz.showvalue(io,V,ib[1],m.v[1])
     for k ∈ list(2,binomial(mdims(V),G))
         @inbounds showterm(io,V,ib[k],m.v[k],compact)
@@ -370,6 +390,7 @@ getindex(m::Multivector,i::T) where T<:AbstractVector = m.v[i]
 setindex!(m::Multivector{V,T} where V,k::T,i::Int,j::Int) where T = (m[i][j] = k)
 Base.firstindex(m::Multivector) = 0
 Base.lastindex(m::Multivector{V,T} where T) where V = mdims(V)
+@pure Base.length(m::Multivector{V}) where V = 1<<mdims(V)
 
 (m::Multivector{V,T})(g::Val{G}) where {V,T,G} = Chain{V,G,T}(m[g])
 (m::Multivector{V,T})(g::Int,i) where {V,T} = m(Val(g),i)
@@ -388,6 +409,7 @@ end
 
 function show(io::IO, m::Multivector{V,T}) where {V,T}
     N,compact,bases = mdims(V),get(io,:compact,false),true
+    io = compactio(io)
     bs = binomsum_set(N)
     @inbounds print(io,m.v[1])
     for i ∈ list(2,N+1)
@@ -492,6 +514,7 @@ for pinor ∈ (:Spinor,:AntiSpinor)
         #setindex!(m::$pinor{V,T} where V,k::T,i::Int,j::Int) where T = (m[i][j] = k)
         Base.firstindex(m::$pinor) = 0
         Base.lastindex(m::$pinor{V,T} where T) where V = mdims(V)
+        @pure Base.length(m::$pinor{V}) where V = 1<<(mdims(V)-1)
         grade(m::$pinor,g::Val) = m(g)
         (m::$pinor{V,T})(g::Int,i::Int) where {T,V} = m(Val(g),i)
         Multivector(t::$pinor{V}) where V = Multivector{V}(t)
@@ -643,6 +666,7 @@ end
 
 function Base.show(io::IO, m::Spinor{V,T}) where {V,T}
     N,compact = mdims(V),get(io,:compact,false)
+    io = compactio(io)
     bs = spinsum_set(N)
     @inbounds print(io,m.v[1])
     for i ∈ evens(2,N)
@@ -654,6 +678,7 @@ function Base.show(io::IO, m::Spinor{V,T}) where {V,T}
 end
 function Base.show(io::IO, m::AntiSpinor{V,T}) where {V,T}
     N,compact = mdims(V),get(io,:compact,false)
+    io = compactio(io)
     bs = antisum_set(N)
     for i ∈ evens(1,N)
         ib = indexbasis(N,i)
@@ -866,6 +891,7 @@ for couple ∈ (:Couple,:PseudoCouple)
         Spinor{V}(val::$couple{V,B,𝕂}) where {V,B,𝕂} = Spinor{V,𝕂}(val)
         Base.zero(::$couple{V,B,T}) where {V,B,T} = $couple{V,B}(zero(Complex{T}))
         Base.zero(::Type{$couple{V,B,T}}) where {V,B,T} = $couple{V,B}(zero(Complex{T}))
+        @pure Base.length(m::$couple) = 2
     end
 end
 
@@ -1007,6 +1033,7 @@ getindex(P::Dyadic,i::Int,j::Int) = P.x[i]*P.y[j]
 show(io::IO,P::Dyadic) = print(io,"(",P.x,")⊗(",P.y,")")
 
 DyadicChain(P::Dyadic{V}) where V = DyadicProduct{V}(P)
+#DyadicChain{V}(P::Dyadic{V}) where V = outer(P.x,P.y)
 DyadicChain{V}(P::Dyadic{V}) where V = DyadicProduct{V}(p)
 DyadicProduct(P::Dyadic{V}) where V = DyadicProduct{V}(P)
 DyadicProduct{V}(P::Dyadic{V}) where V = outer(P.x,P.y)
@@ -1184,27 +1211,27 @@ end
 @pure maxpseudograde(t::Type{<:AntiSpinor{V}}) where V = mdims(V)-1
 @pure nextmaxpseudograde(t) = maxpseudograde(t)-nextgrade(t)
 
-Leibniz.gdims(t::Tuple{Vector{<:Chain},Vector{Int}}) = gdims(t[1][findall(x->!iszero(x),t[2])])
-function Leibniz.gdims(t::Vector{<:Chain})
+Leibniz.count_gdims(t::Tuple{Vector{<:Chain},Vector{Int}}) = count_gdims(t[1][findall(x->!iszero(x),t[2])])
+function Leibniz.count_gdims(t::Vector{<:Chain})
     out = @inbounds zeros(Variables{mdims(Manifold(points(t)))+1,Int})
     @inbounds out[mdims(Manifold(t))+1] = length(t)
     return out
 end
-function Leibniz.gdims(t::Values{N,<:Vector}) where N
+function Leibniz.count_gdims(t::Values{N,<:Vector}) where N
     out = @inbounds zeros(Variables{mdims(points(t[1]))+1,Int})
     for i ∈ list(1,N)
         @inbounds out[mdims(Manifold(t[i]))+1] = length(t[i])
     end
     return out
 end
-function Leibniz.gdims(t::Values{N,<:Tuple}) where N
+function Leibniz.count_gdims(t::Values{N,<:Tuple}) where N
     out = @inbounds zeros(Variables{mdims(points(t[1][1]))+1,Int})
     for i ∈ list(1,N)
         @inbounds out[mdims(Manifold(t[i][1]))+1] = length(t[i][1])
     end
     return out
 end
-function Leibniz.gdims(t::Multivector{V}) where V
+function Leibniz.count_gdims(t::Multivector{V}) where V
     N = mdims(V)
     out = zeros(Variables{N+1,Int})
     bs = binomsum_set(N)
@@ -1217,8 +1244,8 @@ function Leibniz.gdims(t::Multivector{V}) where V
     return out
 end
 
-Leibniz.χ(t::Values{N,<:Vector}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
-Leibniz.χ(t::Values{N,<:Tuple}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
+Leibniz.χ(t::Values{N,<:Vector}) where N = (B=count_gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
+Leibniz.χ(t::Values{N,<:Tuple}) where N = (B=count_gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
 
 ## Adjoint
 
diff --git a/src/products.jl b/src/products.jl
index f4716cc..4de00e0 100644
--- a/src/products.jl
+++ b/src/products.jl
@@ -714,8 +714,8 @@ contraction(a::Proj{V,<:Chain{V,1,<:TensorNested}},b::TensorGraded{V,0}) where V
 #contraction(a::Chain{W,1,<:Proj{V}},b::Chain{V,1}) where {W,V} = Chain{W,1}(value(a).⋅b)
 contraction(a::Chain{W,1,<:Dyadic{V}},b::Chain{V,1}) where {W,V} = Chain{W,1}(value(a).⋅Ref(b))
 contraction(a::Proj{W,<:Chain{W,1,<:TensorNested{V}}},b::Chain{V,1}) where {W,V} = a.v:b
-contraction(a::Chain{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(value.(Ref(a).⋅value(b)))
-contraction(a::Chain{W,L,<:Chain},b::Chain{V,G,<:Chain{W,L}}) where {W,L,G,V} = Chain{V,G}(value.(Ref(a).⋅value(b)))
+contraction(a::Chain{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(value.(Single.(Ref(a).⋅value(b))))
+contraction(a::Chain{W,L,<:Chain,N},b::Chain{V,G,<:Chain{W,L},M}) where {W,L,G,V,N,M} = Chain{V,G}(value.(Ref(a).⋅value(b)))
 contraction(a::Multivector{W,<:Multivector},b::Multivector{V,<:Multivector{W}}) where {W,V} = Multivector{V}(column(Ref(a).⋅value(b)))
 Base.:(:)(a::Chain{V,1,<:Chain},b::Chain{V,1,<:Chain}) where V = sum(value(a).⋅value(b))
 Base.:(:)(a::Chain{W,1,<:Dyadic{V}},b::Chain{V,1}) where {W,V} = sum(value(a).⋅Ref(b))
@@ -725,7 +725,8 @@ contraction(a::Submanifold{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(
 contraction(a::Single{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(column(Ref(a).⋅value(b)))
 contraction(x::Chain{V,G,<:Chain},y::Single{V,G}) where {V,G} = value(y)*x[bladeindex(mdims(V),UInt(basis(y)))]
 contraction(x::Chain{V,G,<:Chain},y::Submanifold{V,G}) where {V,G} = x[bladeindex(mdims(V),UInt(y))]
-contraction(a::Chain{V,L,<:Chain{V,G}},b::Chain{V,G,<:Chain{V}}) where {V,G,L} = Chain{V,G}(matmul(value(a),value(b)))
+contraction(a::Chain{V,L,<:Chain{V,G},N},b::Chain{V,G,<:Chain{V},M}) where {V,G,L,N,M} = Chain{V,G}(contraction.(Ref(a),value(b)))
+contraction(x::Chain{V,L,<:Chain{V,G},N},y::Chain{V,G,<:Chain{V,L},N}) where {L,N,V,G} = Chain{V,G}(contraction.(Ref(x),value(y)))
 contraction(x::Chain{W,L,<:Chain{V,G},N},y::Chain{V,G,T,N}) where {W,L,N,V,G,T} = Chain{V,G}(matmul(value(x),value(y)))
 contraction(x::Chain{W,L,<:Multivector{V},N},y::Chain{V,G,T,N}) where {W,L,N,V,G,T} = Multivector{V}(matmul(value(x),value(y)))
 contraction(x::Multivector{W,<:Chain{V,G},N},y::Multivector{V,T,N}) where {W,N,V,G,T} = Chain{V,G}(matmul(value(x),value(y)))