improved Leibniz algebra and ∇,Δ operators

README.md

@@ -10,3 +10,32 @@
 [![Liberapay patrons](https://img.shields.io/liberapay/patrons/chakravala.svg)](https://liberapay.com/chakravala)
 
 Cross-compatibility of [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) with [Reduce.jl](https://github.com/chakravala/Reduce.jl) for multivariable differential operators and tensor field operations.
+
+```Julia
+julia> using Reduce, Leibniz, Grassmann
+Reduce (Free CSL version, revision 4980), 06-May-19 ...
+
+julia> V = V"3" # load Reduce 1st! otherwise slow
+⟨+++⟩
+
+julia> V(∇)
+∂₁v₁ + ∂₂v₂ + ∂₃v₃
+
+julia> V(∇^0), V(∇^2)
+(1v, (∂₁² + ∂₂² + ∂₃²)v)
+
+julia> V(∇^3)
+(∂₁³ + ∂₂²∂₁ + ∂₃²∂₁)v₁ + (∂₁²∂₂ + ∂₂³ + ∂₃²∂₂)v₂ + (∂₁²∂₃ + ∂₂²∂₃ + ∂₃³)v₃
+
+julia> V(∇^4)
+((∂₁² + ∂₂² + ∂₃²) ^ 2)v
+
+julia> ∇^2 == Δ
+true
+
+julia> ∇, Δ
+(∂ₖvₖ, ∂ₖ²v)
+
+```
+
+The package 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.

src/Leibniz.jl

@@ -6,43 +6,81 @@ module Leibniz
 using DirectSum, StaticArrays, Requires
 using LinearAlgebra, AbstractTensors, AbstractLattices
 import Base: *, ^, +, -, /, show
+import DirectSum: value
 
-export Differential, Partial, Derivation, d, ∂, ∇, Δ
+export Differential, Partial, Derivation, d, ∂, ∇, Δ, @operator
 
 abstract type Operator end
 
-struct Differential{D,O} <: Operator end
+struct Differential{D,O,T} <: Operator
+    v::T
+end
+
 Differential{D}() where D = Differential{D,1}()
+Differential{D,O}() where {D,O} = Differential{D,O}(true)
+Differential{D,O}(v::T) where {D,O,T} = Differential{D,O,T}(v)
+
+value(d::Differential{D,O,T} where {D,O}) where T = d.v
 
-show(io::IO,::Differential{D,0} where D) = print(io,1)
-show(io::IO,::Differential{D,1}) where D = print(io,'∂',DirectSum.subs[D])
-show(io::IO,::Differential{D,O}) where {D,O} = print(io,'∂',DirectSum.subs[D],DirectSum.sups[O])
+show(io::IO,d::Differential{D,0,Bool} where D) = print(io,value(d) ? 1 : -1)
+show(io::IO,d::Differential{D,0} where D) = print(io,value(d))
+show(io::IO,d::Differential{D,1,Bool}) where D = print(io,value(d) ? "" : "-",'∂',DirectSum.subs[D])
+show(io::IO,d::Differential{D,1}) where D = print(io,value(d),'∂',DirectSum.subs[D])
+show(io::IO,d::Differential{D,O,Bool}) where {D,O} = print(io,value(d) ? "" : "-",'∂',DirectSum.subs[D],DirectSum.sups[O])
+show(io::IO,d::Differential{D,O}) where {D,O} = print(io,value(d),'∂',DirectSum.subs[D],DirectSum.sups[O])
 
-+(d::Differential) = d
++(d::O) where O<:Operator = d
 +(r,d::O) where O<:Operator = d+r
++(a::Differential{D,O},b::Differential{D,O}) where {D,O} = (c=a.v+b.v; iszero(c) ? 0 : Differential{D,O}(c))
+-(a::Differential{D,O},b::Differential{D,O}) where {D,O} = (c=a.v-b.v; iszero(c) ? 0 : Differential{D,O}(c))
+-(d::Differential{D,O}) where {D,O} = Differential{D,O}(-d.v)
 *(d::Differential{D,0} where D,r) = r
 *(r,d::Differential) = d*r
-*(::Differential{D,1},::Differential{D,1}) where D = Differential{D,2}()
-*(::Differential{D,O1},::Differential{D,O2}) where {D,O1,O2} = Differential{D,O1+O2}()
-^(::Differential{D,O},o::T) where {D,O,T<:Integer} = Differential{D,O*o}()
+*(a::Differential{D,1},b::Differential{D,1}) where D = (c=a.v*b.v; iszero(c) ? 0 : Differential{D,2}(c))
+*(a::Differential{D,A},b::Differential{D,B}) where {D,A,B} = (c=a.v*b.v; iszero(c) ? 0 : Differential{D,A+B}(c))
+*(a::Differential{D,O,Bool},b::I) where {D,O,I<:Number} = isone(b) ? a : Differential{D,O,I}(value(a) ? b : -b)
+*(a::Differential{D,O,T},b::I) where {D,O,T,I<:Number} = isone(b) ? a : Differential{D,O}(value(a)*b)
+^(d::Differential{D,O},o::T) where {D,O,T<:Integer} = iszero(o) ? 1 : Differential{D,O*o}(value(d)^o)
+^(d::Differential{D,O,Bool},o::T) where {D,O,T<:Integer} = iszero(o) ? 1 : Differential{D,O*o}(value(d) ? true : iseven(o))
+
+struct Partial{D,T} <: Operator
+    v::T
+end
+
+Partial{D}() where D = Partial{D}(true)
+Partial{D}(v::T) where {D,T} = Partial{D,T}(v)
 
-struct Partial{D} <: Operator end
+value(d::Partial{D,T} where D) where T = d.v
 
-show(io::IO,::Partial{D}) where D = print(io,'∂',[DirectSum.subs[k] for k ∈ DirectSum.indices(UInt(D))]...)
-show(io::IO,::Partial{0}) = print(io,1)
-show(io::IO,::Partial{UInt(0)}) = print(io,1)
+show(io::IO,d::Partial{D,Bool}) where D = print(io,value(d) ? "" : "-",'∂',[DirectSum.subs[k] for k ∈ DirectSum.indices(UInt(D))]...)
+show(io::IO,d::Partial{D}) where D = print(io,value(d),'∂',[DirectSum.subs[k] for k ∈ DirectSum.indices(UInt(D))]...)
+show(io::IO,d::Partial{0,Bool}) = print(io,value(d) ? 1 : -1)
+show(io::IO,d::Partial{0}) = print(io,value(d))
+show(io::IO,d::Partial{UInt(0),Bool}) = print(io,value(d) ? 1 : -1)
+show(io::IO,d::Partial{UInt(0)}) = print(io,value(d))
 
 indexint(D) = DirectSum.bit2int(DirectSum.indexbits(max(D...),D))
 
 ∂(D::T) where T<:Integer = Differential{D}()
 ∂(D::T...) where T<:Integer = Partial{indexint(D)}()
 
-*(::Differential{D1,1},::Differential{D2,1}) where {D1,D2} = ∂(D1,D2)
+*(r,d::Partial) = d*r
+*(a::Differential{D1,1},b::Differential{D2,1}) where {D1,D2} = (c=a.v*b.v; iszero(c) ? 0 : Partial{indexint((D1,D2))}(c))
+*(a::Partial{D,Bool},b::I) where {D,I<:Number} = isone(b) ? a : Partial{D,I}(value(a) ? b : -b)
+*(a::Partial{D,T},b::I) where {D,T,I<:Number} = isone(b) ? a : Partial{D}(value(a)*b)
++(a::Partial{D},b::Partial{D}) where D = (c=a.v+b.v; iszero(c) ? 0 : Partial{indexint((D1,D2))}(c))
+-(a::Partial{D},b::Partial{D}) where D = (c=a.v-b.v; iszero(c) ? 0 : Partial{indexint((D1,D2))}(c))
+#-(d::Partial{D,Bool}) where {D,O} = Partial{D,Bool}(!value(d))
+-(d::Partial{D}) where D = Partial{D}(-value(d))
 
 struct OperatorExpr{T} <: Operator
     expr::T
 end
 
+macro operator(expr)
+    OperatorExpr(expr)
+end
+
 show(io::IO,d::OperatorExpr) = print(io,'(',d.expr,')')
 
 add(d,n) = OperatorExpr(Expr(:call,:+,d,n))
@@ -74,7 +112,7 @@ end
 +(d::Differential,n::Differential) = add(d,n)
 +(d::OperatorExpr,n) = plus(d,n)
 +(d::OperatorExpr,n::O) where O<:Operator = plus(d,n)
--(d::O) where O<:Operator = OperatorExpr(Expr(:call,:-,d))
+-(d::OperatorExpr) = OperatorExpr(Expr(:call,:-,d))
 
 #add(d,n) = OperatorExpr(Expr(:call,:+,d,n))
 
@@ -102,8 +140,16 @@ end
 
 *(d::OperatorExpr,n::Differential) = times(d,n)
 *(d::OperatorExpr,n::Partial) = times(d,n)
-*(n::Differential,d::OperatorExpr) = times(d,n)
-*(n::Partial,d::OperatorExpr) = times(d,n)
+*(d::OperatorExpr,n::OperatorExpr) = OperatorExpr(Expr(:call,:*,d,n))
+
+*(a::Differential,b::Differential{D,O}) where {D,O} = Differential{D,O}(a*OperatorExpr(b.v))
+*(a::Differential{D,O},b::OperatorExpr) where {D,O} = Differential{D,O}(value(a)*b.expr)
+*(a::Partial,b::Partial{D}) where D = Partial{D}(a*OperatorExpr(b.v))
+*(a::Partial{D},b::OperatorExpr) where {D,O} = Partial{D}(value(a)*b.expr)
+*(a::Differential,b::Partial{D}) where D = Partial{D}(a*OperatorExpr(b.v))
+*(a::Partial,b::Differential{D}) where D = Differential{D}(a*OperatorExpr(b.v))
+
+^(d::OperatorExpr,n::T) where T<:Integer = iszero(n) ? 1 : isone(n) ? d : OperatorExpr(Expr(:call,:^,d,n))
 
 ## generic
 
@@ -117,16 +163,16 @@ 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ₖ")
-show(io::IO,v::Derivation{T,O}) where {T,O} = print(io,v.v.λ,"∂ₖ",O==1 ? "" : DirectSum.sups[O],"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) ? "ₖ" : "")
 
--(v::Derivation{Bool,O}) where {T,O} = Derivation{Bool,O}(LinearAlgebra.UniformScaling{Bool}(!v.v.λ))
--(v::Derivation{T,O}) where {T,O} = Derivation{T,O}(LinearAlgebra.UniformScaling{T}(-v.v.λ))
+-(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.λ))
 
 function ^(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}(LinearAlgebra.UniformScaling{t}(x))
+    Derivation{t,O*n}(UniformScaling{t}(x))
 end
 
 for op ∈ (:+,:-,:*)
@@ -156,10 +202,8 @@ end
 ∇ = Derivation(LinearAlgebra.I)
 Δ = ∇^2
 
-include("grassmann.jl")
-
 function __init__()
-    #@require Grassmann="4df31cd9-4c27-5bea-88d0-e6a7146666d8" include("grassmann.jl")
+    @require Grassmann="4df31cd9-4c27-5bea-88d0-e6a7146666d8" include("grassmann.jl")
     @require Reduce="93e0c654-6965-5f22-aba9-9c1ae6b3c259" include("symbolic.jl")
     #@require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" nothing
 end

src/grassmann.jl

@@ -4,8 +4,9 @@
 
 using Grassmann
 
-function (V::Signature)(d::Derivation{T,O}) where {T,O,S}
-    sum([(∂(k)^O)*getbasis(V,1<<(k-1)) for k ∈ 1:ndims(V)])
+function (V::Signature{N})(d::Derivation{T,O}) where {N,T,O}
+    x = sum([(∂(k)^2) for k ∈ 1:N])^div(isodd(O) ? O-1 : O,2)
+    isodd(O) ? sum([(x*∂(k))*getbasis(V,1<<(k-1)) for k ∈ 1:N]) : x*getbasis(V,0)
 end
 
 ∂(ω::T) where T<:TensorAlgebra{V} where V = ω⋅V(∇)