Implement quadratic forms

src/StarAlgebras.jl

@@ -8,6 +8,8 @@ import MutableArithmetics as MA
 export StarAlgebra, AlgebraElement
 export basis, coeffs, star
 
+function star end
+
 # AbstractCoefficients
 ## abstract definitions
 include("coefficients.jl")
@@ -31,6 +33,7 @@ include("algebra_elts.jl")
 include("star.jl")
 
 include("arithmetic.jl")
+include("quadratic_form.jl")
 include("show.jl")
 
 # augmented basis implementation

src/algebra_elts.jl

@@ -37,5 +37,6 @@ function LinearAlgebra.dot(w::AbstractVector, b::AlgebraElement)
     return LinearAlgebra.dot(w, coeffs(b))
 end
 function LinearAlgebra.dot(a::AlgebraElement, w::AbstractVector)
+    @assert key_type(basis(parent(a))) <: Integer
     return LinearAlgebra.dot(coeffs(a), w)
 end

src/arithmetic.jl

@@ -60,10 +60,6 @@
     end
 end
 
-function Base.:*(args::Vararg{AlgebraElement,N}) where {N}
-    return MA.operate_to!(_preallocate_output(*, args...), *, args...)
-end
-
 Base.:^(a::AlgebraElement, p::Integer) = Base.power_by_squaring(a, p)
 
 # mutable API
@@ -130,34 +126,34 @@
     X::AlgebraElement,
     Y::AlgebraElement,
 )
-    @assert parent(res) === parent(X) === parent(Y)
+    @assert parent(res) == parent(X)
+    @assert parent(X) == parent(Y)
     MA.operate_to!(coeffs(res), -, coeffs(X), coeffs(Y))
     return res
 end
 
 function MA.operate_to!(
     res::AlgebraElement,
     ::typeof(*),
-    args::Vararg{AlgebraElement,N},
+    ABC::Vararg{AlgebraElement,N},
 ) where {N}
-    for arg in args
-        @assert parent(res) == parent(arg)
-    end
+    @assert allequal(parent, (res, ABC...))
     mstr = mstructure(basis(res))
-    MA.operate_to!(coeffs(res), mstr, coeffs.(args)...)
+    MA.operate_to!(coeffs(res), mstr, map(coeffs, ABC)..., true)
     return res
 end
 
 function MA.operate!(
-    ::UnsafeAddMul{typeof(*)},
+    mul::UnsafeAddMul,
     res::AlgebraElement,
-    args::Vararg{AlgebraElement,N},
+    ABC::Vararg{AlgebraElement,N},
 ) where {N}
-    for arg in args
-        @assert parent(res) == parent(arg)
-    end
-    mstr = mstructure(basis(res))
-    MA.operate!(UnsafeAddMul(mstr), coeffs(res), coeffs.(args)...)
+    MA.operate!(mul, coeffs(res), map(coeffs, ABC)..., true)
+    return res
+end
+
+function MA.operate!(add::UnsafeAdd, res::AlgebraElement, A::AlgebraElement)
+    MA.operate!(add, coeffs(res), coeffs(A))
     return res
 end
 

src/coefficients.jl

@@ -115,7 +115,7 @@ function LinearAlgebra.dot(w::AbstractVector, ac::AbstractCoefficients)
     if isempty(values(ac))
         return zero(MA.promote_sum_mul(eltype(w), value_type(ac)))
     else
-        return sum(w[i] * star(v) for (i, v) in nonzero_pairs(ac))
+        return sum(w[i] * v for (i, v) in nonzero_pairs(ac))
     end
 end
 

src/mstructures.jl

@@ -31,52 +31,74 @@ When the product is not representable faithfully,
 """
 abstract type MultiplicativeStructure end
 
-"""
-    struct UnsafeAddMul{M<:Union{typeof(*),MultiplicativeStructure}}
-        structure::M
-    end
-
-The value of `(op::UnsafeAddMul)(a, b, c)` is `a + structure(b, c)`
-where `a` is not expected to be canonicalized before the operation `+`
-and should not be expected to be canonicalized after either.
-"""
 struct UnsafeAddMul{M<:Union{typeof(*),MultiplicativeStructure}}
     structure::M
 end
 
-function MA.operate_to!(res, ms::MultiplicativeStructure, args::Vararg{Any,N}) where {N}
+"""
+    operate_to!(res, ms::MultiplicativeStructure, A, B[, α = true])
+Compute `α·A·B` storing the result in `res`. Return `res`.
+
+`res` is assumed to behave like `AbstractCoefficients` and not aliased with any
+other arguments.
+`A` and `B` are assumed to behave like `AbstractCoefficients`, while `α` will
+be treated as a scalar.
+
+Canonicalization of the result happens only once at the end of the operation.
+"""
+function MA.operate_to!(
+    res,
+    ms::MultiplicativeStructure,
+    args::Vararg{Any,N},
+) where {N}
     if any(Base.Fix1(===, res), args)
-        throw(ArgumentError("No alias allowed"))
+        throw(
+            ArgumentError(
+                "Aliasing arguments in multiplication is not supported",
+            ),
+        )
     end
     MA.operate!(zero, res)
-    MA.operate!(UnsafeAddMul(ms), res, args...)
+    res = MA.operate!(UnsafeAddMul(ms), res, args...)
     MA.operate!(canonical, res)
     return res
 end
 
-function MA.operate!(::UnsafeAddMul, res, c)
-    for (k, v) in nonzero_pairs(c)
+struct UnsafeAdd end
+
+function MA.operate!(::UnsafeAdd, res, b)
+    for (k, v) in nonzero_pairs(b)
         unsafe_push!(res, k, v)
     end
     return res
 end
 
-function MA.operate!(op::UnsafeAddMul, res, b, c, args::Vararg{Any, N}) where {N}
-    for (kb, vb) in nonzero_pairs(b)
-        for (kc, vc) in nonzero_pairs(c)
-            for (k, v) in nonzero_pairs(op.structure(kb, kc))
+function MA.operate!(op::UnsafeAddMul, res, A, B, α)
+    for (kA, vA) in nonzero_pairs(A)
+        for (kB, vB) in nonzero_pairs(B)
+            for (k, v) in nonzero_pairs(op.structure(kA, kB))
+                cfs = MA.@rewrite α * vA * vB * v
                 MA.operate!(
-                    op,
+                    UnsafeAdd(),
                     res,
-                    SparseCoefficients((_key(op.structure, k),), (vb * vc * v,)),
-                    args...,
+                    SparseCoefficients((_key(op.structure, k),), (cfs,)),
                 )
             end
         end
     end
     return res
 end
 
+function MA.operate!(op::UnsafeAddMul, res, A, B, C, α)
+    for (kA, vA) in nonzero_pairs(A)
+        for (kB, vB) in nonzero_pairs(B)
+            cfs = MA.@rewrite α * vA * vB
+            MA.operate!(op, res, op.structure(kA, kB), C, cfs)
+        end
+    end
+    return res
+end
+
 struct DiracMStructure{Op} <: MultiplicativeStructure
     op::Op
 end

src/quadratic_form.jl

@@ -0,0 +1,84 @@
+"""
+    QuadraticForm(Q)
+    QuadraticForm{star}(Q)
+A simple wrapper for representing a quadratic form.
+
+`QuadraticForm(Q)` represents an honest quadratic form, however in the context
+of `*`-algebras a more canonical object is a **sesquilinear form** which can be
+constructed as `QuadraticForm{star}(Q)`.
+
+Both objects are defined by matrix coefficients with respect to their bases `b`,
+i.e. their values at `i`-th and `j`-th basis elements:
+    `star.(b)·Q·b = ∑ᵢⱼ Q[i,j]·star(b[i])·b[j]`.
+
+# Necessary methods:
+ * `basis(Q)` - a **finite** basis `b` of the quadratic form;
+ * `Base.getindex(Q, i::T, j::T)` - the value at `i`-th and `j`-th basis elements
+   where `T` is the type of indicies of `basis(Q)`;
+ * `Base.eltype(Q)` - the type of `Q[i,j]·star(b[i])·b[j]`.
+"""
+struct QuadraticForm{T,involution}
+    Q::T
+    QuadraticForm(Q) = new{typeof(Q),identity}(Q)
+    QuadraticForm{star}(Q) = new{typeof(Q),star}(Q)
+end
+
+Base.eltype(qf::QuadraticForm) = eltype(qf.Q)
+basis(qf::QuadraticForm) = basis(qf.Q)
+Base.getindex(qf::QuadraticForm, i::T, j::T) where {T} = qf.Q[i, j]
+
+function MA.operate_to!(
+    res,
+    ms::MultiplicativeStructure,
+    Q::QuadraticForm,
+    args...,
+)
+    MA.operate!(zero, res)
+    MA.operate!(UnsafeAddMul(ms), res, Q, args...)
+    MA.operate!(canonical, res)
+    return res
+end
+
+function MA.operate!(
+    op::UnsafeAddMul,
+    res,
+    Q::QuadraticForm{T,ε},
+) where {T,ε}
+    for (i, b1) in pairs(basis(Q))
+        b1★ = ε(b1)
+        for (j, b2) in pairs(basis(Q))
+            MA.operate!(op, res, coeffs(b1★), coeffs(b2), Q[i, j])
+        end
+    end
+end
+
+function MA.operate!(
+    op::UnsafeAddMul,
+    res,
+    Q::QuadraticForm{T,ε},
+    p,
+) where {T,ε}
+    for (i, b1) in pairs(basis(Q))
+        b1★ = ε(b1)
+        for (j, b2) in pairs(basis(Q))
+            MA.operate!(op, res, coeffs(b1★), coeffs(b2), coeffs(p), Q[i, j])
+        end
+    end
+end
+
+"""
+    AlgebraElement(qf::QuadraticForm, A::AbstractStarAlgebra)
+Construct an algebra element in `A` representing quadratic form `qf`.
+
+!!! warning
+    It is assumed that all basis elements of `qf` belong to `A`, or at least
+    that `keys` of their coefficients can be found in the basis of `A`.
+"""
+function AlgebraElement(qf::QuadraticForm, A::AbstractStarAlgebra)
+    @assert all(b -> parent(b) == A, basis(qf))
+    res = zero(eltype(qf), A)
+    MA.operate_to!(coeffs(res), mstructure(A), qf)
+    return res
+end
+
+(A::AbstractStarAlgebra)(qf::QuadraticForm) = AlgebraElement(qf, A)

test/caching_allocations.jl

@@ -31,10 +31,10 @@ end
 
     @test Y isa AlgebraElement
 
-    @static if v"1.10" ≤ VERSION < v"1.11"
+    @static if v"1.10" ≤ VERSION
         star(Y)
         star(Y)
-        @test (@allocations star(Y)) ≤ 4
+        @test (@allocations star(Y)) ≤ 6
     end
 
     @test SA.supp(Y) == basis(fRG)[1:k]
@@ -48,7 +48,7 @@ end
     @test all(!iszero(mstr.table))
     @test_throws SA.UndefRefError all(!iszero, SA.mstructure(fRG).table)
 
-    @static if v"1.10" ≤ VERSION < v"1.11"
+    @static if v"1.10" ≤ VERSION
         @test (@allocations Y * Y) > k^2 - 2 * k
         @test Y * Y isa AlgebraElement
         @test (@allocations Y * Y) ≤ 26
@@ -63,15 +63,15 @@ end
     @test all(!iszero, SA.mstructure(fRG).table)
 
 
-    @static if v"1.10" ≤ VERSION < v"1.11"
+    @static if v"1.10" ≤ VERSION
         YY = deepcopy(Y)
         _alloc_test(YY, *, Y, Y, 25)
         _alloc_test(YY, +, Y, Y, 0)
         _alloc_test(YY, -, Y, Y, 0)
 
         # SparseArrays calls `Base.unalias` which allocates:
-        _alloc_test(YY, +, YY, Y, 2)
-        _alloc_test(YY, +, Y, YY, 2)
+        _alloc_test(YY, +, YY, Y, 4)
+        _alloc_test(YY, +, Y, YY, 4)
     end
 end
 

test/monoid_algebra.jl

@@ -127,17 +127,22 @@
             @test_throws ArgumentError MA.operate_to!(d, *, d, b)
 
             d = deepcopy(a)
-            @test @allocated(MA.operate_to!(d, *, 2, d)) == 0
-            @test d == 2a
-
+            MA.operate_to!(d, *, 2, a) # preallocates memory in coefficients
             @test @allocated(MA.operate_to!(d, *, 2, a)) == 0
             @test d == 2a
 
-            MA.operate!(zero, d)
-            MA.operate!(SA.UnsafeAddMul(*), d, a, b, b)
-            MA.operate!(SA.canonical, SA.coeffs(d))
-            @test a * b^2 == *(a, b, b)
-            @test d == *(a, b, b)
+            @test @allocated(MA.operate_to!(d, *, 2, d)) == 0
+            @test d == 4a
+
+            MA.operate_to!(d, *, a, b)
+            @test d == a * b
+            MA.operate!(SA.UnsafeAddMul(SA.mstructure(RG)), d, a, b)
+            @test d == 2 * a * b
+
+            MA.operate_to!(d, *, a, b, b)
+            @test d == a * b * b
+            MA.operate!(SA.UnsafeAddMul(SA.mstructure(RG)), d, a, b, b)
+            @test d == 2 * a * b * b
         end
     end
 end