diff --git a/appveyor.yml b/.appveyor.yml
similarity index 96%
rename from appveyor.yml
rename to .appveyor.yml
index 52b5405..3c9ce83 100644
--- a/appveyor.yml
+++ b/.appveyor.yml
@@ -4,6 +4,8 @@ environment:
   - julia_version: 1.1
   - julia_version: 1.2
   - julia_version: 1.3
+  - julia_version: 1.4
+  - julia_version: 1.5
   - julia_version: nightly
 
 platform:
diff --git a/.travis.yml b/.travis.yml
index e0d4a06..f653184 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -9,6 +9,8 @@ julia:
   - 1.1
   - 1.2
   - 1.3
+  - 1.4
+  - 1.5
   - nightly
 matrix:
   allow_failures:
diff --git a/Project.toml b/Project.toml
index ebdf1e7..c15bee3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,19 +1,17 @@
 name = "Leibniz"
 uuid = "edad4870-8a01-11e9-2d75-8f02e448fc59"
 authors = ["Michael Reed"]
-version = "0.0.5"
+version = "0.1.0"
 
 [deps]
 AbstractTensors = "a8e43f4a-99b7-5565-8bf1-0165161caaea"
-DirectSum = "22fd7b30-a8c0-5bf2-aabe-97783860d07c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 
 [compat]
 julia = "1"
-DirectSum = "0"
-AbstractTensors = "0"
-StaticArrays = "0"
+Combinatorics = "1"
+AbstractTensors = "0.5"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index ff8fa1d..0b78426 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # Leibniz.jl
 
-*Operator algebras for multivariate differentiable Julia expressions*
+*Bit entanglements for tensor algebra derivations and hypergraphs*
 
 [![Build Status](https://travis-ci.org/chakravala/Leibniz.jl.svg?branch=master)](https://travis-ci.org/chakravala/Leibniz.jl)
 [![Build status](https://ci.appveyor.com/api/projects/status/xb03dyfvhni6vrj5?svg=true)](https://ci.appveyor.com/project/chakravala/leibniz-jl)
@@ -9,11 +9,26 @@
 [![Gitter](https://badges.gitter.im/Grassmann-jl/community.svg)](https://gitter.im/Grassmann-jl/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge)
 [![Liberapay patrons](https://img.shields.io/liberapay/patrons/chakravala.svg)](https://liberapay.com/chakravala)
 
-Compatibility of [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) for multivariable differential operators and tensor field operations.
+Although intended for compatibility use with the [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) package for multivariable differential operators and tensor field operations, `Leibniz` can be used independently.
+
+### Extended dual index printing with full alphanumeric characters #62'
+
+To help provide a commonly shared and readable indexing to the user, some print methods are provided:
+```julia
+julia> Leibniz.printindices(stdout,Leibniz.indices(UInt(2^62-1)),false,"v")
+v₁₂₃₄₅₆₇₈₉₀abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ
+
+julia> Leibniz.printindices(stdout,Leibniz.indices(UInt(2^62-1)),false,"w")
+w¹²³⁴⁵⁶⁷⁸⁹⁰ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz
+```
+An application of this is in `Grassmann` and `DirectSum`, where dual indexing is used.
+
+# Derivation
+
+Generates the tensor algebra of multivariable symmetric Leibniz differentials and interfaces `using Reduce, Grassmann` to provide the `∇,Δ` vector field operators, enabling  mixed-symmetry tensors with arbitrary multivariate `Grassmann` manifolds.
 
 ```Julia
 julia> using Leibniz, Grassmann
-Reduce (Free CSL version, revision 4980), 06-May-19 ...
 
 julia> V = tangent(ℝ^3,4,3)
 ⟨+++⟩
@@ -37,6 +52,4 @@ julia> ∇, Δ
 (∂ₖvₖ, ∂ₖ²v)
 ```
 
-Generates the tensor algebra of multivariable symmetric Leibniz differentials and interfaces `using Reduce, Grassmann` to provide the `∇,Δ` vector field operators, enabling  mixed-symmetry tensors with arbitrary multivariate `Grassmann` manifolds.
-
 This is an initial undocumented pre-release registration for testing with other packages.
diff --git a/src/Leibniz.jl b/src/Leibniz.jl
index 21aabfe..c1578e9 100644
--- a/src/Leibniz.jl
+++ b/src/Leibniz.jl
@@ -1,171 +1,90 @@
 module Leibniz
 
-#   This file is part of Leibniz.jl. It is licensed under the GPL license
+#   This file is part of Leibniz.jl. It is licensed under the AGPL license
 #   Leibniz Copyright (C) 2019 Michael Reed
 
-using DirectSum, StaticArrays #, Requires
 using LinearAlgebra, AbstractTensors
-import Base: *, ^, +, -, /, \, show, zero
-import DirectSum: value, V0, mixedmode, pre, diffmode
+export Manifold, Differential, Derivation, d, ∂, ∇, Δ
+import Base: getindex, convert, @pure, +, *, ∪, ∩, ⊆, ⊇, ==, show, zero
+import LinearAlgebra: det, rank
 
-export Differential, Monomial, Derivation, d, ∂, ∇, Δ, @operator
+## Manifold{N}
 
-abstract type Operator{V} end #<: TensorAlgebra{V} end
+import AbstractTensors: TensorAlgebra, Manifold, TensorGraded, scalar, isscalar, involute
+import AbstractTensors: vector, isvector, bivector, isbivector, volume, isvolume, ⋆, mdims
+import AbstractTensors: value, valuetype, interop, interform, even, odd, isnull, norm
+import AbstractTensors: TupleVector, Values, Variables, FixedVector, SVector, MVector
+import AbstractTensors: basis, complementleft, complementlefthodge, unit, involute, clifford
+abstract type TensorTerm{V,G} <: TensorGraded{V,G} end
 
-abstract type Polynomial{V,G} <: Operator{V} end
+## utilities
 
-+(d::O) where O<:Operator = d
-+(r,d::O) where O<:Operator = d+r
+include("utilities.jl")
 
-struct Monomial{V,G,D,O,T} <: Polynomial{V,G}
-    v::T
-end
+#="""
+    floatprecision(s)
 
-Monomial{V,G,D}() where {V,G,D} = Monomial{V,G,D}(true)
-Monomial{V,G,D}(v::T) where {V,G,D,T} = Monomial{V,G,D,1,T}(v)
-Monomial{V,G,D,O}() where {V,G,D,O,T} = Monomial{V,G,D,O}(true)
-Monomial{V,G,D,O}(v::T) where {V,G,D,O,T} = Monomial{V,G,D,O,T}(v)
-
-zero(::Monomial) = Monomial{V0,0,0,0}()
-
-value(d::Monomial{V,G,D,T} where {V,G,D}) where T = d.v
-
-sups(O) = O ≠ 1 ? DirectSum.sups[O] : ""
-
-show(io::IO,d::Monomial{V,G,D,O,Bool} where G) where {V,D,O} = print(io,value(d) ? "" : "-",pre[mixedmode(V)>0 ? 4 : 3],[DirectSum.subs[k] for k ∈ DirectSum.shift_indices(V,UInt(D))]...,sups(O))
-show(io::IO,d::Monomial{V,G,D,O} where G) where {V,D,O} = print(io,value(d),pre[mixedmode(V)>0 ? 4 : 3],[DirectSum.subs[k] for k ∈ DirectSum.shift_indices(V,UInt(D))]...,sups(O))
-show(io::IO,d::Monomial{V,G,D,0,Bool} where {V,G,D}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,G,0,O,Bool} where {V,G,O}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,0,D,O,Bool} where {V,D,O}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,G,D,0} where {V,G,D}) = print(io,value(d))
-show(io::IO,d::Monomial{V,G,0} where {V,G}) = print(io,value(d))
-show(io::IO,d::Monomial{V,0} where V) = print(io,value(d))
-show(io::IO,d::Monomial{V,G,D,UInt(0),Bool} where {V,G,D}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,G,UInt(0),O,Bool} where {V,G,O}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,UInt(0),D,O,Bool} where {V,D,O}) = print(io,value(d) ? 1 : -1)
-show(io::IO,d::Monomial{V,G,D,UInt(0)} where {V,G,D}) = print(io,value(d))
-show(io::IO,d::Monomial{V,G,UInt(0)} where {V,G}) = print(io,value(d))
-show(io::IO,d::Monomial{V,UInt(0)} where V) = print(io,value(d))
-
-indexint(D) = DirectSum.bit2int(DirectSum.indexbits(max(D...),D))
-
-#∂(D::T...) where T<:Integer = Monomial{V0,length(D),indexint(D)}()
-#∂(V::S,D::T...) where {S<:Manifold,T<:Integer} = Monomial{V,length(D),indexint(D)}()
-
-*(r,d::Monomial) = d*r
-*(d::Monomial{V,G,D,0} where {V,G,D},r) = r
-*(d::Monomial{V,G,0,O} where {V,G,O},r) = r
-*(d::Monomial{V,0,D,O} where {V,D,O},r) = r
-function *(a::Monomial{V,1,D,O1},b::Monomial{V,1,D,O2}) where {V,D,O1,O2}
-    O = O1+O2
-    O > diffmode(V) && (return 0)
-    c = a.v*b.v
-    iszero(c) ? 0 : Monomial{V,2,D,O1+O2}(c)
-end
-function *(a::Monomial{V,1,D,1},b::Monomial{V,1,D,1}) where {V,D,O1,O2}
-    2 > diffmode(V) && (return 0)
-    c = a.v*b.v
-    iszero(c) ? 0 : Monomial{V,2,D,2}(c)
-end
-function *(a::Monomial{V,1,D1,1},b::Monomial{V,1,D2,1}) where {V,D1,D2}
-    2 > diffmode(V) && (return 0)
-    c = a.v*b.v
-    iszero(c) ? 0 : Monomial{V,2,D1|D2}(c)
-end
-*(a::Monomial{V,G,D,O,Bool},b::I) where {V,G,D,O,I<:Number} = isone(b) ? a : Monomial{V,G,D,O,I}(value(a) ? b : -b)
-*(a::Monomial{V,G,D,O,T},b::I) where {V,G,D,O,T,I<:Number} = isone(b) ? a : Monomial{V,G,D,O}(value(a)*b)
-+(a::Monomial{V,G,D,O},b::Monomial{V,G,D,O}) where {V,G,D,O} = (c=a.v+b.v; iszero(c) ? 0 : Monomial{V,G,D,O}(c))
--(a::Monomial{V,G,D,O},b::Monomial{V,G,D,O}) where {V,G,D,O} = (c=a.v-b.v; iszero(c) ? 0 : Monomial{V,G,D,O}(c))
-#-(d::Monomial{V,G,D,O,Bool}) where {V,G,D,O} = Monomial{V,G,D,O,Bool}(!value(d))
--(d::Monomial{V,G,D,O}) where {V,G,D,O} = Monomial{V,G,D,O}(-value(d))
-function ^(d::Monomial{V,G,D,O},o::T) where {V,G,D,O,T<:Integer}
-    Oo = O*o
-    GOo = G+Oo
-    GOo > diffmode(V) && (return 0)
-    iszero(o) ? 1 : Monomial{V,GOo,D,Oo}(value(d)^o)
-end
-function ^(d::Monomial{V,G,D,O,Bool},o::T) where {V,G,D,O,T<:Integer}
-    Oo = O*o
-    GOo = G+Oo
-    GOo > diffmode(V) && (return 0)
-    iszero(o) ? 1 : Monomial{V,GOo,D,Oo}(value(d) ? true : iseven(o))
-end
+Set float precision for display Float64 coefficents.
 
-struct OperatorExpr{T} <: Operator{T}
-    expr::T
-end
+Float coefficients `f` are printed as `@sprintf(s,f)`.
 
-macro operator(expr)
-    OperatorExpr(expr)
-end
+If `s == ""` (default), then `@sprintf` is not used.
+"""
+const floatprecision = ( () -> begin
+        gs::String = ""
+        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
+    end)()
+export floatprecision
+
+macro fprintf()
+    s = floatprecision()
+    isempty(s) ? :(m.v) : :(Printf.@sprintf($s,m.v))
+end=#
+
+# symbolic print types
 
-show(io::IO,d::OperatorExpr) = print(io,'(',d.expr,')')
-
-add(d,n) = OperatorExpr(Expr(:call,:+,d,n))
-
-function plus(d::OperatorExpr{T},n) where T
-    iszero(n) && (return d)
-    if T == Expr
-        if d.expr.head == :call
-            if d.expr.args[1] == :+
-                return OperatorExpr(Expr(:call,:+,push!(copy(d.expr.args[2:end]),n)...))
-            elseif d.expr.args[1] == :- && length(d.expr.args) == 2 && d.expr.args[2] == n
-                return 0
-            else
-                return OperatorExpr(Expr(:call,:+,d.expr,n))
-            end
-        else
-            throw(error("Operator expression not implemented"))
-        end
-    else
-        OperatorExpr(d.expr+n)
+parval = (Expr,Complex,Rational,TensorAlgebra)
+
+# number fields
+
+const Fields = (Real,Complex)
+const Field = Fields[1]
+const ExprField = Union{Expr,Symbol}
+
+extend_field(Field=Field) = (global parval = (parval...,Field))
+
+for T ∈ Fields
+    @eval begin
+        Base.:(==)(a::T,b::TensorTerm{V,G} where V) where {T<:$T,G} = G==0 ? a==value(b) : 0==a==value(b)
+        Base.:(==)(a::TensorTerm{V,G} where V,b::T) where {T<:$T,G} = G==0 ? value(a)==b : 0==value(a)==b
     end
 end
 
-+(d::Monomial,n::T) where T<:Number = iszero(n) ? d : add(d,n)
-+(d::Monomial,n::Monomial) = add(d,n)
-+(d::OperatorExpr,n) = plus(d,n)
-+(d::OperatorExpr,n::O) where O<:Operator = plus(d,n)
--(d::OperatorExpr) = OperatorExpr(Expr(:call,:-,d))
-
-#add(d,n) = OperatorExpr(Expr(:call,:+,d,n))
-
-function times(d::OperatorExpr{T},n) where T
-    iszero(n) && (return 0)
-    isone(n) && (return d)
-    if T == Expr
-        if d.expr.head == :call
-            if d.expr.args[1] ∈ (:+,:-)
-                return OperatorExpr(Expr(:call,:+,(d.expr.args[2:end] .* Ref(n))...))
-            elseif d.expr.args[1] == :*
-                (d.expr.args[2]*n)*d.expr.args[3] + d.expr.args[2]*(d.expr.args[3]*n)
-            elseif d.expr.args[1] == :/
-                (d.expr.args[2]*n)/d.expr.args[3] - (d.expr.args[2]*(d.expr.args[3]*n))/(d.expr.args[3]^2)
-            else
-                return OperatorExpr(Expr(:call,:*,d.expr,n))
-            end
-        else
-            throw(error("Operator expression not implemented"))
-        end
-    else
-        OperatorExpr(d.expr*n)
+Base.:(==)(a::TensorTerm,b::TensorTerm) = 0 == value(a) == value(b)
+
+for T ∈ (Fields...,Symbol,Expr)
+    @eval begin
+        Base.isapprox(a::S,b::T) where {S<:TensorAlgebra,T<:$T} = Base.isapprox(a,Simplex{Manifold(a)}(b))
+        Base.isapprox(a::S,b::T) where {S<:$T,T<:TensorAlgebra} = Base.isapprox(b,a)
     end
 end
 
-*(d::OperatorExpr,n::Monomial) = times(d,n)
-*(d::OperatorExpr,n::OperatorExpr) = OperatorExpr(Expr(:call,:*,d,n))
-*(n::T,d::OperatorExpr) where T<:Number = OperatorExpr(DirectSum.∏(n,d.expr))
+## fundamentals
+
+"""
+    getbasis(V::Manifold,v)
 
-*(a::Monomial{V,1},b::Monomial{V,1,D,O}) where {V,D,O} = Monomial{V,1,D,O}(a*OperatorExpr(b.v))
-*(a::Monomial,b::Monomial{V,G,D,O}) where {V,G,D,O} = Monomial{V,G,D,O}(a*OperatorExpr(b.v))
-*(a::Monomial{V,G,D,O},b::OperatorExpr) where {V,G,D,O} = Monomial{V,G,D,O}(value(a)*b.expr)
+Fetch a specific `SubManifold{G,V}` element from an optimal `SubAlgebra{V}` selection.
+"""
+@inline getbasis(V,b) = getbasis(V,UInt(b))
 
-^(d::OperatorExpr,n::T) where T<:Integer = iszero(n) ? 1 : isone(n) ? d : OperatorExpr(Expr(:call,:^,d,n))
+Base.one(V::T) where T<:TensorGraded = one(Manifold(V))
+Base.zero(V::T) where T<:TensorGraded = zero(Manifold(V))
 
-## generic
+@pure g_one(::Type{T}) where T = one(T)
+@pure g_zero(::Type{T}) where T = zero(T)
 
-Base.signbit(::O) where O<:Operator = false
-Base.abs(d::O) where O<:Operator = d
+## Derivation
 
 struct Derivation{T,O}
     v::UniformScaling{T}
@@ -174,19 +93,19 @@ end
 Derivation{T}(v::UniformScaling{T}) where T = Derivation{T,1}(v)
 Derivation(v::UniformScaling{T}) where T = Derivation{T}(v)
 
-show(io::IO,v::Derivation{Bool,O}) where O = print(io,(v.v.λ ? "" : "-"),"∂ₖ",O==1 ? "" : DirectSum.sups[O],"v",isodd(O) ? "ₖ" : "")
-show(io::IO,v::Derivation{T,O}) where {T,O} = print(io,v.v.λ,"∂ₖ",O==1 ? "" : DirectSum.sups[O],"v",isodd(O) ? "ₖ" : "")
+show(io::IO,v::Derivation{Bool,O}) where O = print(io,(v.v.λ ? "" : "-"),"∂ₖ",O==1 ? "" : AbstractTensors.sups[O],"v",isodd(O) ? "ₖ" : "")
+show(io::IO,v::Derivation{T,O}) where {T,O} = print(io,v.v.λ,"∂ₖ",O==1 ? "" : AbstractTensors.sups[O],"v",isodd(O) ? "ₖ" : "")
 
--(v::Derivation{Bool,O}) where {T,O} = Derivation{Bool,O}(UniformScaling{Bool}(!v.v.λ))
--(v::Derivation{T,O}) where {T,O} = Derivation{T,O}(UniformScaling{T}(-v.v.λ))
+Base.:-(v::Derivation{Bool,O}) where {T,O} = Derivation{Bool,O}(UniformScaling{Bool}(!v.v.λ))
+Base.:-(v::Derivation{T,O}) where {T,O} = Derivation{T,O}(UniformScaling{T}(-v.v.λ))
 
-function ^(v::Derivation{T,O},n::S) where {T,O,S<:Integer}
+function Base.:^(v::Derivation{T,O},n::S) where {T,O,S<:Integer}
     x = T<:Bool ? (isodd(n) ? v.v.λ : true ) : v.v.λ^n
     t = typeof(x)
     Derivation{t,O*n}(UniformScaling{t}(x))
 end
 
-for op ∈ (:+,:-,:*)
+for op ∈ (:(Base.:+),:(Base.:-),:(Base.:*))
     @eval begin
         $op(a::Derivation{A,O},b::Derivation{B,O}) where {A,B,O} = Derivation{promote_type(A,B),O}($op(a.v,b.v))
         $op(a::Derivation{A,O},b::B) where {A,B<:Number,O} = Derivation{promote_type(A,B),O}($op(a.v,b))
@@ -196,16 +115,16 @@ end
 
 unitype(::UniformScaling{T}) where T = T
 
-/(a::Derivation{A,O},b::Derivation{B,O}) where {A,B,O} = (x=a.v/b.v; Derivation{unitype(x),O}(x))
-/(a::Derivation{A,O},b::B) where {A,B<:Number,O} = (x=a.v/b; Derivation{unitype(x),O}(x))
-#/(a::A,b::Derivation{B,O}) where {A<:Number,B,O} = (x=a/b.v; Derivation{typeof(x),O}(x))
-\(a::Derivation{A,O},b::Derivation{B,O}) where {A,B,O} = (x=a.v\b.v; Derivation{unitype(x),O}(x))
-\(a::A,b::Derivation{B,O}) where {A<:Number,B,O} = (x=a\b.v; Derivation{unitype(x),O}(x))
+Base.:/(a::Derivation{A,O},b::Derivation{B,O}) where {A,B,O} = (x=Base.:/(a.v,b.v); Derivation{unitype(x),O}(x))
+Base.:/(a::Derivation{A,O},b::B) where {A,B<:Number,O} = (x=Base.:/(a.v,b); Derivation{unitype(x),O}(x))
+#Base.:/(a::A,b::Derivation{B,O}) where {A<:Number,B,O} = (x=Base.:/(a,b.v); Derivation{typeof(x),O}(x))
+Base.:\(a::Derivation{A,O},b::Derivation{B,O}) where {A,B,O} = (x=a.v\b.v; Derivation{unitype(x),O}(x))
+Base.:\(a::A,b::Derivation{B,O}) where {A<:Number,B,O} = (x=a\b.v; Derivation{unitype(x),O}(x))
 
 import AbstractTensors: ∧, ∨
 import LinearAlgebra: dot, cross
 
-for op ∈ (:+,:-,:*,:/,:\,:∧,:∨,:dot,:cross)
+for op ∈ (:(Base.:+),:(Base.:-),:(Base.:*),:(Base.:/),:(Base.:\),:∧,:∨,:dot,:cross)
     @eval begin
         $op(a::Derivation,b::B) where B<:TensorAlgebra = $op(Manifold(b)(a),b)
         $op(a::A,b::Derivation) where A<:TensorAlgebra = $op(a,Manifold(a)(b))
@@ -216,6 +135,16 @@ const ∇ = Derivation(LinearAlgebra.I)
 const Δ = ∇^2
 
 function d end
+function ∂ end
+
+include("generic.jl")
+include("operations.jl")
+include("indices.jl")
+
+bladeindex(cache_limit,one(UInt))
+basisindex(cache_limit,one(UInt))
+
+indexbasis(Int((sparse_limit+cache_limit)/2),1)
 
 #=function __init__()
     @require Reduce="93e0c654-6965-5f22-aba9-9c1ae6b3c259" include("symbolic.jl")
diff --git a/src/generic.jl b/src/generic.jl
new file mode 100644
index 0000000..333ec96
--- /dev/null
+++ b/src/generic.jl
@@ -0,0 +1,193 @@
+
+#   This file is part of Leibniz.jl. It is licensed under the AGPL license
+#   Leibniz Copyright (C) 2019 Michael Reed
+
+export bits, basis, grade, order, options, metric, polymode, dyadmode, diffmode, diffvars
+export valuetype, value, hasinf, hasorigin, isorigin, norm, indices, tangent, isbasis, ≅
+
+@pure grade(V::M) where M<:Manifold{N} where N = N-(isdyadic(V) ? 2 : 1)*diffvars(V)
+@pure grade(m::T) where T<:Real = 0
+@pure order(m) = 0
+@pure order(V::M) where M<:Manifold = diffvars(V)
+@pure options(::Int) = 0
+@pure metric(::Int) = zero(UInt)
+@pure metric(V::M,b::UInt) where M<:Manifold = PROD(V[indices(b)])
+@pure polymode(::Int) = true
+@pure dyadmode(::Int) = 0
+@pure diffmode(::Int) = 0
+@pure diffvars(::Int) = 0
+for mode ∈ (:options,:polymode,:dyadmode,:diffmode,:diffvars)
+    @eval @pure $mode(t::T) where T<:TensorAlgebra = $mode(Manifold(t))
+end
+
+@pure ≅(a,b) = grade(a) == grade(b) && order(a) == order(b) && diffmode(a) == diffmode(b)
+
+export isdyadic, isdual, istangent
+const mixedmode = dyadmode
+@pure isdyadic(t::T) where T<:TensorAlgebra = dyadmode(Manifold(t))<0
+@pure isdual(t::T) where T<:TensorAlgebra = dyadmode(Manifold(t))>0
+@pure istangent(t::T) where T<:TensorAlgebra = diffvars(Manifold(t))≠0
+
+@pure isbasis(x) = false
+@pure isdyadic(::Int) = false
+@pure isdual(::Int) = false
+@pure istangent(::Int) = false
+
+@inline value_diff(m::T) where T<:TensorTerm = (v=value(m);istensor(v) ? v : m)
+
+for T ∈ (:T,:(Type{T}))
+    @eval @pure isbasis(::$T) where T<:TensorAlgebra = false
+end
+@pure UInt(m::T) where T<:TensorTerm = UInt(basis(m))
+
+@pure hasconformal(V) = hasinf(V) && hasorigin(V)
+@pure hasinf(::Int) = false
+@pure hasinf(t::M) where M<:Manifold = hasinf(Manifold(t))
+#@pure hasinf(m::T) where T<:TensorAlgebra = hasinf(Manifold(m))
+@pure hasorigin(::Int) = false
+@pure hasorigin(t::M) where M<:Manifold = hasorigin(Manifold(t))
+#@pure hasorigin(m::T) where T<:TensorAlgebra = hasorigin(Manifold(m))
+
+@pure hasorigin(V,B::UInt) = hasinf(V) ? (Bits(2)&B)==Bits(2) : isodd(B)
+
+@pure function hasinf(V,A::UInt,B::UInt)
+    hasconformal(V) && (isodd(A) || isodd(B))
+end
+@pure function hasorigin(V,A::UInt,B::UInt)
+    hasconformal(V) && (hasorigin(V,A) || hasorigin(V,B))
+end
+
+@pure function hasinf2(V,A::UInt,B::UInt)
+    hasconformal(V) && isodd(A) && isodd(B)
+end
+@pure function hasorigin2(V,A::UInt,B::UInt)
+    hasconformal(V) && hasorigin(V,A) && hasorigin(V,B)
+end
+
+@pure function hasorigininf(V,A::UInt,B::UInt)
+    hasconformal(V) && hasorigin(V,A) && isodd(B) && !hasorigin(V,B) && !isodd(A)
+end
+@pure function hasinforigin(V,A::UInt,B::UInt)
+    hasconformal(V) && isodd(A) && hasorigin(V,B) && !isodd(B) && !hasorigin(V,A)
+end
+
+@pure function hasinf2origin(V,A::UInt,B::UInt)
+    hasinf2(V,A,B) && hasorigin(V,A,B)
+end
+@pure function hasorigin2inf(V,A::UInt,B::UInt)
+    hasorigin2(V,A,B) && hasinf(V,A,B)
+end
+
+@pure diffmask(::Int) = zero(UInt)
+@pure function diffmask(V::M) where M<:Manifold
+    d = diffvars(V)
+    if isdyadic(V)
+        v = ((one(UInt)<<d)-1)<<(mdims(V)-2d)
+        w = ((one(UInt)<<d)-1)<<(mdims(V)-d)
+        return d<0 ? (typemax(UInt)-v,typemax(UInt)-w) : (v,w)
+    end
+    v = ((one(UInt)<<d)-1)<<(mdims(V)-d)
+    d<0 ? typemax(UInt)-v : v
+end
+
+@pure function symmetricsplit(V,a)
+    sm,dm = symmetricmask(V,a),diffmask(V)
+    isdyadic(V) ? (sm&dm[1],sm&dm[2]) : sm
+end
+
+@pure function symmetricmask(V,a)
+    d = diffmask(V)
+    a&(isdyadic(V) ? |(d...) : d)
+end
+
+@pure function symmetricmask(V,a,b)
+    d = diffmask(V)
+    D = isdyadic(V) ? |(d...) : d
+    aD,bD = (a&D),(b&D)
+    return a&~D, b&~D, aD|bD, aD&bD
+end
+
+@pure function diffcheck(V,A::UInt,B::UInt)
+    d,db = diffvars(V),diffmask(V)
+    v = isdyadic(V) ? db[1]|db[2] : db
+    hi = hasinf2(V,A,B) && !hasorigin(V,A,B)
+    ho = hasorigin2(V,A,B) && !hasinf(V,A,B)
+    (hi || ho) || (d≠0 && count_ones(A&v)+count_ones(B&v)>diffmode(V))
+end
+
+## functors
+
+@pure loworder(N::Int) = N
+function supermanifold end
+
+## adjoint parities
+
+@pure parityreverse(G) = isodd(Int((G-1)*G/2))
+@pure parityinvolute(G) = isodd(G)
+@pure parityclifford(G) = parityreverse(G)⊻parityinvolute(G)
+const parityconj = parityreverse
+
+## reverse
+
+import Base: reverse, ~
+export involute, clifford
+
+@pure grade_basis(V::Int,B) = B&(one(UInt)<<V-1)
+@pure grade_basis(V,B) = B&(one(UInt)<<grade(V)-1)
+@pure grade(V,B) = count_ones(grade_basis(V,B))
+
+@doc """
+    ~(ω::TensorAlgebra)
+
+Reverse of an element: ~ω = (-1)^(grade(ω)*(grade(ω)-1)/2)*ω
+""" Base.conj
+#reverse(a::UniformScaling{Bool}) = UniformScaling(!a.λ)
+#reverse(a::UniformScaling{T}) where T<:Field = UniformScaling(-a.λ)
+
+"""
+    reverse(ω::TensorAlgebra)
+
+Reverse of an element: ~ω = (-1)^(grade(ω)*(grade(ω)-1)/2)*ω
+"""
+@inline Base.:~(b::TensorAlgebra) = Base.conj(b)
+#@inline ~(b::UniformScaling) = reverse(b)
+
+@doc """
+    involute(ω::TensorAlgebra)
+
+Involute of an element: ~ω = (-1)^grade(ω)*ω
+""" involute
+
+@doc """
+    clifford(ω::TensorAlgebra)
+
+Clifford conjugate of an element: clifford(ω) = involute(conj(ω))
+""" clifford
+
+odd(t::T) where T<:TensorGraded{V,G} where {V,G} = parityinvolute(G) ? t : zero(V)
+even(t::T) where T<:TensorGraded{V,G} where {V,G} = parityinvolute(G) ? zero(V) : t
+
+"""
+    imag(ω::TensorAlgebra)
+
+The `imag` part `(ω-(~ω))/2` is defined by `abs2(imag(ω)) == -(imag(ω)^2)`.
+"""
+Base.imag(t::T) where T<:TensorGraded{V,G} where {V,G} = parityreverse(G) ? t : zero(V)
+
+"""
+real(ω::TensorAlgebra)
+
+The `real` part `(ω+(~ω))/2` is defined by `abs2(real(ω)) == real(ω)^2`.
+"""
+Base.real(t::T) where T<:TensorGraded{V,G} where {V,G} = parityreverse(G) ? zero(V) : t
+
+Base.isfinite(b::T) where T<:TensorTerm = isfinite(value(b))
+
+# comparison (special case for scalars)
+
+Base.isless(a::T,b::S) where {T<:TensorTerm{V,0},S<:TensorTerm{W,0}} where {V,W} = isless(value(a),value(b))
+Base.isless(a::T,b) where T<:TensorTerm{V,0} where V = isless(value(a),b)
+Base.isless(a,b::T) where T<:TensorTerm{V,0} where V = isless(a,value(b))
+Base.:<=(x::T,y::S) where {T<:TensorTerm{V,0},S<:TensorTerm{W,0}} where {V,W} = isless(x,y) | (x == y)
+Base.:<=(x::T,y) where T<:TensorTerm{V,0} where V = isless(x,y) | (x == y)
+Base.:<=(x,y::T) where T<:TensorTerm{V,0} where V = isless(x,y) | (x == y)
diff --git a/src/indices.jl b/src/indices.jl
new file mode 100644
index 0000000..409880a
--- /dev/null
+++ b/src/indices.jl
@@ -0,0 +1,237 @@
+
+#   This file is part of Leibniz.jl. It is licensed under the AGPL license
+#   Leibniz Copyright (C) 2019 Michael Reed
+
+# alpha-numeric digits
+const vio = ('∞','∅')
+const digs = "1234567890"
+const low_case,upp_case = "abcdefghijklmnopqrstuvwxyz","ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+const low_greek,upp_greek = "αβγδϵζηθικλμνξοπρστυφχψω","ΑΒΓΔΕΖΗΘΙΚΛΜΝΞΟΡΣΤΥΦΨΩ"
+const alphanumv = digs*low_case*upp_case #*low_greek*upp_greek
+const alphanumw = digs*upp_case*low_case #*upp_greek*low_greek
+
+# subscript index
+const subs = Dict{Int,Char}(
+   -1 => vio[1],
+    0 => vio[2],
+    1 => '₁',
+    2 => '₂',
+    3 => '₃',
+    4 => '₄',
+    5 => '₅',
+    6 => '₆',
+    7 => '₇',
+    8 => '₈',
+    9 => '₉',
+    10 => '₀',
+    [j=>alphanumv[j] for j ∈ 11:36]...
+)
+
+# superscript index
+const sups = Dict{Int,Char}(
+   -1 => vio[1],
+    0 => vio[2],
+    1 => '¹',
+    2 => '²',
+    3 => '³',
+    4 => '⁴',
+    5 => '⁵',
+    6 => '⁶',
+    7 => '⁷',
+    8 => '⁸',
+    9 => '⁹',
+    10 => '⁰',
+    [j=>alphanumw[j] for j ∈ 11:36]...
+)
+
+# vector and co-vector prefix
+const pre = ("v","w","∂","ϵ")
+const PRE = ("X","x","Y","y")
+
+# vector space and dual-space symbols
+const vsn = (:V,:VV,:W)
+const VSN = (:Χ,:ΧΧ,:Υ) # \Chi,\Upsilon
+
+# converts indices into BitArray of length N
+@inline function indexbits(N::I,indices::T) where {I<:Integer,T<:SVTI}
+    out = falses(N)
+    for k ∈ indices
+        out[k] = true
+    end
+    return out
+end
+
+# index sets
+index_limit = 20
+const digits_fast_cache = Vector{SVector}[]
+const digits_fast_extra = Dict{UInt,SVector}[]
+@pure digits_fast_calc(b,N) = SVector{N+1,Int}(digits(b,base=2,pad=N+1))
+@pure function digits_fast(b,N)
+    if N>index_limit
+        n = N-index_limit
+        for k ∈ length(digits_fast_extra)+1:n
+            push!(digits_fast_extra,Dict{UInt,SVector{k+1,Int}}())
+        end
+        !haskey(digits_fast_extra[n],b) && push!(digits_fast_extra[n],b=>digits_fast_calc(b,N))
+        @inbounds digits_fast_extra[n][b]
+    else
+        for k ∈ length(digits_fast_cache)+1:min(N,index_limit)
+            push!(digits_fast_cache,[digits_fast_calc(d,k) for d ∈ 0:1<<(k+1)-1])
+            GC.gc()
+        end
+        @inbounds digits_fast_cache[N][b+1]
+    end
+end
+
+const indices_cache = Dict{UInt,Vector{Int}}()
+indices(b::Bits) = findall(digits(b,base=2).==1)
+function indices_calc(b::UInt,N::Int)
+    d = digits_fast(b,N)
+    l = length(d)
+    a = Int[]
+    for i ∈ 1:l
+        d[i] == 1 && push!(a,i)
+    end
+    return a
+end
+function indices(b::UInt,N::Int)
+    !haskey(indices_cache,b) && push!(indices_cache,b=>indices_calc(b,N))
+    return @inbounds indices_cache[b]
+end
+
+shift_indices(V::T,b::UInt) where T<:Manifold = shift_indices(supermanifold(V),b)
+function shift_indices!(s::T,set::Vector{Int}) where T<:Manifold
+    M = supermanifold(s)
+    if !isempty(set)
+        k = 1
+        hasinf(M) && set[1] == 1 && (set[1] = -1; k += 1)
+        shift = hasinf(M) + hasorigin(M)
+        hasorigin(M) && length(set)>=k && set[k]==shift && (set[k]=0;k+=1)
+        shift > 0 && (set[k:end] .-= shift)
+    end
+    return set
+end
+
+shift_indices(V::Int,b::UInt) = indices(b,V)
+shift_indices!(s::Int,set::Vector{Int}) = set
+
+# printing of indices
+@inline function printindex(i,l::Bool=false,e::String=pre[1],pre=pre)
+    t = i>36; j = t ? i-26 : i
+    (l&&(0<j≤10)) ? j : ((e∉pre[[1,3]])⊻t ? sups[j] : subs[j])
+ end
+ @inline printindices(io::IO,b::UInt,l::Bool=false,e::String=pre[1],pre::NTuple{4,String}=pre) = printindices(io,indices(b),l,e,pre)
+@inline printindices(io::IO,b::VTI,l::Bool=false,e::String=pre[1],pre::NTuple{4,String}=pre) = print(io,e,[printindex(i,l,e,pre) for i ∈ b]...)
+@inline printindices(io::IO,a::VTI,b::VTI,l::Bool=false,e::String=pre[1],f::String=pre[2]) = printindices(io,a,b,Int[],Int[],l,e,f)
+@inline function printindices(io::IO,a::VTI,b::VTI,c::VTI,d::VTI,l::Bool=false,e::String=pre[1],f::String=pre[2],g::String=pre[3],h::String=pre[4])
+    A,B,C,D = isempty(a),!isempty(b),!isempty(c),!isempty(d)
+    PRE = (e,f,g,h)
+    C && printindices(io,c,l,g,PRE)
+    D && printindices(io,d,l,h,PRE)
+    !((B || C || D) && A) && printindices(io,a,l,e,PRE)
+    B && printindices(io,b,l,f,PRE)
+end
+
+printindices(io::IO,V::Int,e::UInt,label::Bool=false) = printlabel(io,V,e,label,pre...)
+@inline function printlabel(io::IO,V::Int,e::UInt,label::Bool,vec,cov,duo,dif)
+    printindices(io,shift_indices(V,e),label,vec)
+    return io
+end
+
+@inline function printlabel(io::IO,V::T,e::UInt,label::Bool,vec,cov,duo,dif) where T<:Manifold
+    M = supermanifold(V)
+    N,D,C,db = mdims(M),diffvars(M),dyadmode(V),diffmask(V)
+    if C < 0
+        es = e & (~(db[1]|db[2]))
+        n = Int((N-2D)/2)
+        eps = shift_indices(V,e & db[1]).-(N-2D-hasinf(M)-hasorigin(M))
+        par = shift_indices(V,e & db[2]).-(N-D-hasinf(M)-hasorigin(M))
+        printindices(io,shift_indices(V,es & Bits(2^n-1)),shift_indices(V,es>>n),eps,par,label,vec,cov,duo,dif)
+    else
+        es = e & (~db)
+        eps = shift_indices(V,e & db).-(N-D-hasinf(M)-hasorigin(M))
+        if !isempty(eps)
+            printindices(io,shift_indices(V,es),Int[],C>0 ? Int[] : eps,C>0 ? eps : Int[],label,C>0 ? cov : vec,cov,C>0 ? dif : duo,dif)
+        else
+            printindices(io,shift_indices(V,es),label,C>0 ? string(cov) : vec)
+        end
+    end
+    return io
+end
+
+@inline printlabel(V::T,e::UInt,label::Bool,vec,cov,duo,dif) where T<:Manifold = printlabel(IOBuffer(),V,e,label,vec,cov,duo,dif) |> take! |> String
+
+function indexstring(V::M,D) where M<:Manifold
+    io = IOBuffer()
+    printlabel(io,V,D,true,PRE...)
+    String(take!(io))
+end
+
+indexsymbol(V::M,D) where M<:Manifold = Symbol(indexstring(V,D))
+
+indexsplit(B,N) = [UInt(1)<<(k-1) for k ∈ indices(B,N)]
+
+indexparity!(ind::SVector{N,Int}) where N = indexparity!(MVector(ind))
+function indexparity!(ind::MVector{N,Int}) where N
+    k = 1
+    t = false
+    while k < length(ind)
+        if ind[k] > ind[k+1]
+            ind[SVector(k,k+1)] = ind[SVector(k+1,k)]
+            t = !t
+            k ≠ 1 && (k -= 1)
+        else
+            k += 1
+        end
+    end
+    return t, ind
+end
+function indexparity!(ind::Vector{Int},s)
+    k = 1
+    t = false
+    while k < length(ind)
+        if ind[k] == ind[k+1]
+            ind[k] == 1 && hasinf(s) && (return t, ind, true)
+            s[ind[k]] && (t = !t)
+            deleteat!(ind,[k,k+1])
+        elseif ind[k] > ind[k+1]
+            ind[[k,k+1]] = ind[[k+1,k]]
+            t = !t
+            k ≠ 1 && (k -= 1)
+        else
+            k += 1
+        end
+    end
+    return t, ind, false
+end
+
+@noinline function indexparity(V::T,v::Symbol)::Tuple{Bool,Vector,T,Bool} where T
+    vs = string(v)
+    vt = vs[1:1]≠pre[1]
+    Z=match(Regex("([$(pre[1])]([0-9a-vx-zA-VX-Z]+))?([$(pre[2])]([0-9a-zA-Z]+))?"),vs)
+    ef = String[]
+    for k ∈ (2,4)
+        Z[k] ≠ nothing && push!(ef,Z[k])
+    end
+    length(ef) == 0 && (return false,Int[],V,true)
+    let W = V,fs=false
+        C = dyadmode(V)
+        X = C≥0 && mdims(V)<4sizeof(UInt)+1
+        X && (W = T<:Int ? 2V : (C>0 ? V'⊕V : V⊕V'))
+        V2 = (vt ⊻ (vt ? C≠0 : C>0)) ? V' : V
+        L = length(ef) > 1
+        M = X ? Int(mdims(W)/2) : mdims(W)
+        m = ((!L) && vt && (C<0)) ? M : 0
+        chars = (L || (Z[2] ≠ nothing)) ? alphanumv : alphanumw
+        (es,e,et) = indexparity!([findfirst(isequal(ef[1][k]),chars) for k∈1:length(ef[1])].+m,C<0 ? V : V2)
+        et && (return false,Int[],V,true)
+        w,d = if L
+            (fs,f,ft) = indexparity!([findfirst(isequal(ef[2][k]),alphanumw) for k∈1:length(ef[2])].+M,W)
+            ft && (return false,Int[],V,true)
+            W,[e;f]
+        else
+            V2,e
+        end
+        return es⊻fs, d, w, false
+    end
+end
diff --git a/src/operations.jl b/src/operations.jl
new file mode 100644
index 0000000..93cd181
--- /dev/null
+++ b/src/operations.jl
@@ -0,0 +1,165 @@
+
+#   This file is part of Leibniz.jl. It is licensed under the AGPL license
+#   Leibniz Copyright (C) 2019 Michael Reed
+
+export ⊕, χ, gdims
+
+"""
+    χ(::TensorAlgebra)
+
+Compute the Euler characteristic χ = ∑ₚ(-1)ᵖbₚ.
+"""
+χ(t::T) where T<:TensorAlgebra = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)]))
+χ(t::T) where T<:TensorTerm = χ(Manifold(t),bits(basis(t)),t)
+@inline χ(V,b::UInt,t) = iszero(t) ? 0 : isodd(count_ones(symmetricmask(V,b,b)[1])) ? 1 : -1
+
+function gdims(t::T) where T<:TensorTerm
+    B,N = bits(basis(t)),mdims(t)
+    g = count_ones(symmetricmask(Manifold(t),B,B)[1])
+    MVector{N+1,Int}([g==G ? abs(χ(t)) : 0 for G ∈ 0:N])
+end
+function gdims(t::T) where T<:TensorGraded{V,G} where {V,G}
+    N = mdims(V)
+    out = zeros(MVector{N+1,Int})
+    ib = indexbasis(N,G)
+    for k ∈ 1:length(ib)
+        @inbounds t[k] ≠ 0 && (out[count_ones(symmetricmask(V,ib[k],ib[k])[1])+1] += 1)
+    end
+    return out
+end
+
+## set theory ∪,∩,⊆,⊇
+
+@pure ∪(x::T) where T<:Manifold = x
+@pure ∪(a::A,b::B,c::C...) where {A<:Manifold,B<:Manifold,C<:Manifold} = ∪(a∪b,c...)
+
+@pure ∩(x::T) where T<:Manifold = x
+@pure ∩(a::A,b::B,c::C...) where {A<:Manifold,B<:Manifold,C<:Manifold} = ∩(a∩b,c...)
+
+## conversions
+
+@pure function mixed(V::M,ibk::UInt) where M<:Manifold
+    N,D,VC = mdims(V),diffvars(V),isdual(V)
+    return if D≠0
+        A,B = ibk&(UInt(1)<<(N-D)-1),ibk&diffmask(V)
+        VC ? (A<<(N-D))|(B<<N) : A|(B<<(N-D))
+    else
+        VC ? ibk<<N : ibk
+    end
+end
+
+@pure function combine(v::M,w::T,iak::UInt,ibk::UInt) where {T<:Manifold,M<:Manifold}
+    (isdual(v) ≠ isdual(w)) && throw(error("$v and $w incompatible"))
+    V,W = supermanifold(v),supermanifold(w)
+    return if istangent(V)||istangent(W)
+        gras1,gras2 = iak&(UInt(1)<<grade(V)-1),ibk&(UInt(1)<<grade(W)-1)
+        diffs = (iak&diffmask(W))|(ibk&diffmask(W))
+        gras1|(gras2<<grade(V))|(diffs<<mdims(W)) # A|(B<<(N-D))
+    else
+        iak|(ibk<<mdims(V)) # ibk
+    end
+end
+
+## complement parity
+
+@pure parityrighthodge(V::Int,B::Int,G,N=nothing) = isodd(V)⊻parityright(V,B,G,N)
+@pure paritylefthodge(V::Int,B::Int,G,N) = (isodd(G) && iseven(N)) ⊻ parityrighthodge(V,B,G,N)
+@pure parityright(V::Int,B::Int,G,N=nothing) = isodd(B+Int((G+1)*G/2))
+@pure parityleft(V::Int,B::Int,G,N) = (isodd(G) && iseven(N)) ⊻ parityright(V,B,G,N)
+
+for side ∈ (:left,:right)
+    p = Symbol(:parity,side)
+    pg = Symbol(p,:hodge)
+    pn = Symbol(p,:null)
+    pnp = Symbol(pn,:pre)
+    @eval begin
+        @pure $p(V::UInt,B::UInt,N::Int) = $p(0,sum(indices(B,N)),count_ones(B),N)
+        @pure $pg(V::UInt,B::UInt,N::Int) = $pg(count_ones(V&B),sum(indices(B,N)),count_ones(B),N)
+        @inline function $pn(V,B,v)
+            if hasconformal(V) && count_ones(B&UInt(3)) == 1
+                isodd(B) ? (2v) : (v/2)
+            else
+                v
+            end
+        end
+        @inline function $pnp(V,B,v)
+            if hasconformal(V) && count_ones(B&UInt(3)) == 1
+                isodd(B) ? Expr(:call,:*,2,v) : Expr(:call,:/,v,2)
+            else
+                v
+            end
+        end
+        @pure function $pg(V::Int,B::UInt,G=count_ones(B))
+            ind = indices(B&(UInt(1)<<(mdims(V)-diffvars(V))-1),mdims(V))
+            c = hasconformal(V) && (B&UInt(3) == UInt(2))
+            $p(0,sum(ind),G,mdims(V)-diffvars(V))⊻c ? -1 : 1
+        end
+    end
+end
+
+@pure function complement(N::Int,B::UInt,D::Int=0,P::Int=0)::UInt
+    UP,ND = UInt(1)<<(P==1 ? 0 : P)-1, N-D
+    C = ((~B)&(UP⊻(UInt(1)<<ND-1)))|(B&(UP⊻((UInt(1)<<D-1)<<ND)))
+    count_ones(C&UP)≠1 ? C⊻UP : C
+end
+
+# Hodge star ★
+
+const complementrighthodge = ⋆
+const complementright = !
+
+## complement
+
+export complementleft, complementright, ⋆, complementlefthodge, complementrighthodge
+
+@doc """
+    complementrighthodge(ω::TensorAlgebra)
+
+Grassmann-Poincare-Hodge complement: ⋆ω = ω∗I
+""" complementrighthodge
+
+@doc """
+    complementlefthodge(ω::TensorAlgebra)
+
+Grassmann-Poincare left complement: ⋆'ω = I∗'ω
+""" complementlefthodge
+
+@doc """
+    complementright(::TensorAlgebra)
+
+Non-metric variant of Grassmann-Poincare-Hodge complement.
+""" complementright
+
+@doc """
+    complementleft(::TensorAlgebra)
+
+Non-metric variant Grassmann-Poincare left complement.
+""" complementleft
+
+# QR compatibility
+
+@inline function LinearAlgebra.reflectorApply!(x::AbstractVector, τ::TensorAlgebra, A::StridedMatrix)
+    @assert !LinearAlgebra.has_offset_axes(x)
+    m, n = size(A)
+    if length(x) != m
+        throw(DimensionMismatch("reflector has length $(length(x)), which must match the dimension of matrix A, $m"))
+    end
+    @inbounds begin
+        for j = 1:n
+            #dot
+            vAj = A[1,j]
+            for i = 2:m
+                vAj += Base.conj(x[i])*A[i,j]
+            end
+
+            vAj = conj(τ)*vAj
+
+            #ger
+            A[1, j] -= vAj
+            for i = 2:m
+                A[i,j] -= x[i]*vAj
+            end
+        end
+    end
+    return A
+end
diff --git a/src/symbolic.jl b/src/symbolic.jl
deleted file mode 100644
index 3010244..0000000
--- a/src/symbolic.jl
+++ /dev/null
@@ -1,16 +0,0 @@
-
-#   This file is part of Leibniz.jl. It is licensed under the GPL license
-#   Leibniz Copyright (C) 2019 Michael Reed
-
-using Reduce
-
-Reduce.RExpr(::Monomial{V,G,D,0} where {V,G}) where D = RExpr("1")
-Reduce.RExpr(::Monomial{V,G,D,1} where {V,G}) where D = RExpr("d_$D")
-Reduce.RExpr(::Monomial{V,G,D,O} where {V,G}) where {D,O} = RExpr("d_$D^$O")
-
-for RE ∈ (RExpr,Expr,Symbol)
-    @eval begin
-        *(d::Monomial{V,1,D,1} where V,r::$RE) where D = Algebra.df(r,RExpr("v$D"))
-        *(d::Monomial{V,O,D,O} where V,r::$RE) where {D,O} = Algebra.df(r,RExpr("v$D"),RExpr("$O"))
-    end
-end
diff --git a/src/utilities.jl b/src/utilities.jl
new file mode 100644
index 0000000..86bb955
--- /dev/null
+++ b/src/utilities.jl
@@ -0,0 +1,262 @@
+
+#   This file is part of Leibniz.jl. It is licensed under the AGPL license
+#   Leibniz Copyright (C) 2019 Michael Reed
+
+import AbstractTensors: conj, inv, PROD, SUM, -, /
+import AbstractTensors: sqrt, abs, exp, expm1, log, log1p, sin, cos, sinh, cosh, ^
+
+Bits,bits = UInt,UInt
+
+const VTI = Union{Vector{Int},Tuple,NTuple}
+const SVTI = Union{Vector{Int},Tuple,NTuple,SVector}
+
+bit2int(b::BitArray{1}) = parse(UInt,join(reverse([t ? '1' : '0' for t ∈ b])),base=2)
+
+AbstractTensors.:-(x::Values) = Base.:-(x)
+AbstractTensors.:-(x::Values{N,Any} where N) = broadcast(-,x)
+@inline AbstractTensors.norm(z::Values{N,Any} where N) = sqrt(SUM(z.^2...))
+
+@pure binomial(N,G) = Base.binomial(N,G)
+@pure binomial_set(N) = Values(Int[binomial(N,g) for g ∈ 0:N]...)
+@pure intlog(M::Integer) = Int(log2(M))
+@pure promote_type(t...) = Base.promote_type(t...)
+@pure mvec(N,G,t) = Variables{binomial(N,G),t}
+@pure mvec(N,t) = Variables{2^N,t}
+@pure svec(N,G,t) = FixedVector{binomial(N,G),t}
+@pure svec(N,t) = FixedVector{1<<N,t}
+
+## constructor
+
+@inline assign_expr!(e,x::Vector{Any},v::Symbol,expr) = v ∈ e && push!(x,Expr(:(=),v,expr))
+
+function insert_expr(e,vec=:mvec,T=:(valuetype(a)),S=:(valuetype(b)),L=:(2^N);mv=0)
+    x = Any[] # Any[:(sigcheck(sig(a),sig(b)))]
+    assign_expr!(e,x,:N,:(mdims(V)))
+    assign_expr!(e,x,:M,:(Int(N/2)))
+    assign_expr!(e,x,:t,vec≠:mvec ? :Any : :(promote_type($T,$S)))
+    assign_expr!(e,x,:out,mv≠0 ? :(t=Any;convert(svec(N,Any),out)) : :(zeros($vec(N,t))))
+    assign_expr!(e,x,:r,:(binomsum(N,G)))
+    assign_expr!(e,x,:bng,:(binomial(N,G)))
+    assign_expr!(e,x,:bnl,:(binomial(N,L)))
+    assign_expr!(e,x,:ib,:(indexbasis(N,G)))
+    assign_expr!(e,x,:bs,:(binomsum_set(N)))
+    assign_expr!(e,x,:bn,:(binomial_set(N)))
+    assign_expr!(e,x,:df,:(dualform(V)))
+    assign_expr!(e,x,:di,:(dualindex(V)))
+    assign_expr!(e,x,:D,:(diffvars(V)))
+    assign_expr!(e,x,:μ,:(istangent(V)))
+    assign_expr!(e,x,:P,:(hasinf(V)+hasorigin(V)))
+    return x
+end
+
+## cache
+
+export binomsum, bladeindex, basisindex, indexbasis, lowerbits, expandbits
+
+const algebra_limit = 8
+const sparse_limit = 22
+const cache_limit = 12
+const fill_limit = 0.5
+
+import Combinatorics: combinations
+
+combo_calc(k,g) = collect(combinations(1:k,g))
+const combo_cache = Vector{Vector{Vector{Int}}}[]
+const combo_extra = Vector{Vector{Vector{Int}}}[]
+function combo(n::Int,g::Int)::Vector{Vector{Int}}
+    if g == 0
+        [Int[]]
+    elseif n>sparse_limit
+        N=n-sparse_limit
+        for k ∈ length(combo_extra)+1:N
+            push!(combo_extra,Vector{Vector{Int}}[])
+        end
+        @inbounds for k ∈ length(combo_extra[N])+1:g
+            @inbounds push!(combo_extra[N],Vector{Int}[])
+        end
+        @inbounds isempty(combo_extra[N][g]) && (combo_extra[N][g]=combo_calc(n,g))
+        @inbounds combo_extra[N][g]
+    else
+        for k ∈ length(combo_cache)+1:min(n,sparse_limit)
+            push!(combo_cache,([combo_calc(k,q) for q ∈ 1:k]))
+        end
+        @inbounds combo_cache[n][g]
+    end
+end
+
+binomsum_calc(n) = SVector{n+2,Int}([0;cumsum([binomial(n,q) for q=0:n])])
+const binomsum_cache = (SVector{N,Int} where N)[SVector(0),SVector(0,1)]
+const binomsum_extra = (SVector{N,Int} where N)[]
+@pure function binomsum(n::Int, i::Int)::Int
+    if n>sparse_limit
+        N=n-sparse_limit
+        for k ∈ length(binomsum_extra)+1:N
+            push!(binomsum_extra,SVector{0,Int}())
+        end
+        @inbounds isempty(binomsum_extra[N]) && (binomsum_extra[N]=binomsum_calc(n))
+        @inbounds binomsum_extra[N][i+1]
+    else
+        for k=length(binomsum_cache):n+1
+            push!(binomsum_cache,binomsum_calc(k))
+        end
+        @inbounds binomsum_cache[n+1][i+1]
+    end
+end
+@pure function binomsum_set(n::Int)::(SVector{N,Int} where N)
+    if n>sparse_limit
+        N=n-sparse_limit
+        for k ∈ length(binomsum_extra)+1:N
+            push!(binomsum_extra,SVector{0,Int}())
+        end
+        @inbounds isempty(binomsum_extra[N]) && (binomsum_extra[N]=binomsum_calc(n))
+        @inbounds binomsum_extra[N]
+    else
+        for k=length(binomsum_cache):n+1
+            push!(binomsum_cache,binomsum_calc(k))
+        end
+        @inbounds binomsum_cache[n+1]
+    end
+end
+
+@pure function bladeindex_calc(d,k)
+    H = indices(UInt(d),k)
+    findfirst(x->x==H,combo(k,count_ones(d)))
+end
+const bladeindex_cache = Vector{Int}[]
+const bladeindex_extra = Vector{Int}[]
+@pure function bladeindex(n::Int,s::UInt)::Int
+    if s == 0
+        1
+    elseif n>(index_limit)
+        bladeindex_calc(s,n)
+    elseif n>cache_limit
+        N = n-cache_limit
+        for k ∈ length(bladeindex_extra)+1:N
+            push!(bladeindex_extra,Int[])
+        end
+        @inbounds isempty(bladeindex_extra[N]) && (bladeindex_extra[N]=-ones(Int,1<<n-1))
+        @inbounds signbit(bladeindex_extra[N][s]) && (bladeindex_extra[N][s]=bladeindex_calc(s,n))
+        @inbounds bladeindex_extra[N][s]
+    else
+        j = length(bladeindex_cache)+1
+        for k ∈ j:min(n,cache_limit)
+            push!(bladeindex_cache,(bladeindex_calc.(1:1<<k-1,k)))
+            GC.gc()
+        end
+        @inbounds bladeindex_cache[n][s]
+    end
+end
+
+@inline basisindex_calc(d,k) = binomsum(k,count_ones(UInt(d)))+bladeindex(k,UInt(d))
+const basisindex_cache = Vector{Int}[]
+const basisindex_extra = Vector{Int}[]
+@pure function basisindex(n::Int,s::UInt)::Int
+    if s == 0
+        1
+    elseif n>(index_limit)
+        basisindex_calc(s,n)
+    elseif n>cache_limit
+        N = n-cache_limit
+        for k ∈ length(basisindex_extra)+1:N
+            push!(basisindex_extra,Int[])
+        end
+        @inbounds isempty(basisindex_extra[N]) && (basisindex_extra[N]=-ones(Int,1<<n-1))
+        @inbounds signbit(basisindex_extra[N][s]) && (basisindex_extra[N][s]=basisindex_calc(s,n))
+        @inbounds basisindex_extra[N][s]
+    else
+        j = length(basisindex_cache)+1
+        for k ∈ j:min(n,cache_limit)
+            push!(basisindex_cache,[basisindex_calc(d,k) for d ∈ 1:1<<k-1])
+            GC.gc()
+        end
+        @inbounds basisindex_cache[n][s]
+    end
+end
+
+index2int(k,c) = bit2int(indexbits(k,c))
+indexbasis_calc(k,G) = index2int.(k,combo(k,G))
+const indexbasis_cache = Vector{Vector{UInt}}[]
+const indexbasis_extra = Vector{Vector{UInt}}[]
+@pure function indexbasis(n::Int,g::Int)::Vector{UInt}
+    if n>sparse_limit
+        N = n-sparse_limit
+        for k ∈ length(indexbasis_extra)+1:N
+            push!(indexbasis_extra,Vector{UInt}[])
+        end
+        @inbounds for k ∈ length(indexbasis_extra[N])+1:g
+            @inbounds push!(indexbasis_extra[N],UInt[])
+        end
+        @inbounds if isempty(indexbasis_extra[N][g])
+            @inbounds indexbasis_extra[N][g] = indexbasis_calc(n,g)
+        end
+        @inbounds indexbasis_extra[N][g]
+    else
+        for k ∈ length(indexbasis_cache)+1:n
+            push!(indexbasis_cache,indexbasis_calc.(k,1:k))
+        end
+        @inbounds g>0 ? indexbasis_cache[n][g] : [zero(UInt)]
+    end
+end
+@pure indexbasis(N) = vcat(indexbasis(N,0),indexbasis_set(N)...)
+@pure indexbasis_set(N) = SVector(((N≠0 && N<sparse_limit) ? @inbounds(indexbasis_cache[N]) : Vector{Bits}[indexbasis(N,g) for g ∈ 0:N])...)
+
+# SubManifold
+
+const lowerbits_cache = Vector{Vector{UInt}}[]
+const lowerbits_extra = Dict{UInt,Dict{UInt,UInt}}[]
+@pure lowerbits_calc(N,S,B,k=indices(S,N)) = bit2int(indexbits(N,findall(x->x∈k,indices(B,N))))
+@pure function lowerbits(N,S,B)
+    if N>cache_limit
+        n = N-cache_limit
+        for k ∈ length(lowerbits_extra)+1:n
+            push!(lowerbits_extra,Dict{UInt,Dict{UInt,UInt}}())
+        end
+        @inbounds !haskey(lowerbits_extra[n],S) && push!(lowerbits_extra[n],S=>Dict{UInt,UInt}())
+        @inbounds !haskey(lowerbits_extra[n][S],B) && push!(lowerbits_extra[n][S],B=>lowerbits_calc(N,S,B))
+        @inbounds lowerbits_extra[n][S][B]
+    else
+        for k ∈ length(lowerbits_cache)+1:min(N,cache_limit)
+            push!(lowerbits_cache,Vector{Int}[])
+        end
+        for s ∈ length(lowerbits_cache[N])+1:S
+            k = indices(S,N)
+            push!(lowerbits_cache[N],[lowerbits_calc(N,s,d,k) for d ∈ UInt(0):UInt(1)<<(N+1)-1])
+        end
+        @inbounds lowerbits_cache[N][S][B+1]
+    end
+end
+
+const expandbits_cache = Dict{UInt,Dict{UInt,UInt}}[]
+@pure expandbits_calc(N,S,B) = bit2int(indexbits(N,indices(S,N)[indices(B,N)]))
+@pure function expandbits(N,S,B)
+    for k ∈ length(expandbits_cache)+1:N
+        push!(expandbits_cache,Dict{UInt,Dict{UInt,UInt}}())
+    end
+    @inbounds !haskey(expandbits_cache[N],S) && push!(expandbits_cache[N],S=>Dict{UInt,UInt}())
+    @inbounds !haskey(expandbits_cache[N][S],B) && push!(expandbits_cache[N][S],B=>expandbits_calc(N,S,B))
+    @inbounds expandbits_cache[N][S][B]
+end
+
+#=const expandbits_cache = Vector{Vector{UInt}}[]
+const expandbits_extra = Dict{UInt,Dict{UInt,UInt}}[]
+@pure expandbits_calc(N,S,B,k=indices(S,N)) = bit2int(indexbits(N,k[indices(B,N)]))
+@pure function expandbits(N,S,B)
+    if N>cache_limit
+        n = N-cache_limit
+        for k ∈ length(expandbits_extra)+1:n
+            push!(expandbits_extra,Dict{UInt,Dict{UInt,UInt}}())
+        end
+        @inbounds !haskey(expandbits_extra[n],S) && push!(expandbits_extra[n],S=>Dict{UInt,UInt}())
+        @inbounds !haskey(expandbits_extra[n][S],B) && push!(expandbits_extra[n][S],B=>expandbits_calc(N,S,B))
+        @inbounds expandbits_extra[n][S][B]
+    else
+        for k ∈ length(expandbits_cache)+1:min(N,cache_limit)
+            push!(expandbits_cache,Vector{Int}[])
+        end
+        for s ∈ length(expandbits_cache[N])+1:S
+            k = indices(S,N)
+            push!(expandbits_cache[N],[expandbits_calc(N,s,d,k) for d ∈ UInt(0):UInt(1)<<(N+1)-1])
+        end
+        @inbounds expandbits_cache[N][S][B+1]
+    end
+end=#