From 82d0cc89f0bb43b2cb28a0c83f889e9637e5bf31 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 2 Dec 2024 00:42:12 +0100
Subject: [PATCH] allow efficient both quadratic and sesqui-linear forms

---
 src/StarAlgebras.jl    |  2 ++
 src/mstructures.jl     | 18 +++++++++++++-----
 src/types.jl           |  8 ++++++++
 test/quadratic_form.jl | 12 +++++++++---
 4 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/src/StarAlgebras.jl b/src/StarAlgebras.jl
index 7540279..7d2ef69 100644
--- a/src/StarAlgebras.jl
+++ b/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")
diff --git a/src/mstructures.jl b/src/mstructures.jl
index 1a7dfdf..67197d3 100644
--- a/src/mstructures.jl
+++ b/src/mstructures.jl
@@ -96,10 +96,16 @@ end
 
 """
     QuadraticForm(Q)
+    QuadraticForm{star}(Q)
 A simple wrapper for representing a quadratic form.
 
-A `QuadraticForm(Q)` represents actually a **sesquilinear form** with respect to
-its basis `b`, i.e. `star.(b)·Q·b = ∑ᵢⱼ Q[i,j]·star(b[i])·b[j]`.
+`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;
@@ -107,17 +113,19 @@ its basis `b`, i.e. `star.(b)·Q·b = ∑ᵢⱼ Q[i,j]·star(b[i])·b[j]`.
    where `T` is the type of indicies of `basis(Q)`;
  * `Base.eltype(Q)` - the type of `Q[i,j]`.
 """
-struct QuadraticForm{T}
+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!(op::UnsafeAddMul, res, Q::QuadraticForm)
+function MA.operate!(op::UnsafeAddMul, res, Q::QuadraticForm{T,ε}) where {T,ε}
     for (i, b1) in pairs(basis(Q))
-        b1★ = star(b1)
+        b1★ = ε(b1)
         for (j, b2) in pairs(basis(Q))
             MA.operate!(op, res, coeffs(b1★), coeffs(b2); cfs = Q[i, j])
         end
diff --git a/src/types.jl b/src/types.jl
index 3b57d88..bf59271 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -132,6 +132,14 @@ function AlgebraElement{T}(X::AlgebraElement) where {T}
     return AlgebraElement(w, parent(X))
 end
 
+"""
+    AlgebraElement(qf::QuadraticForm, A::AbstractStarAlgebra)
+Construct an algebra element in `A` representing quadratic form.
+
+!!! 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)
diff --git a/test/quadratic_form.jl b/test/quadratic_form.jl
index a1b2b66..d5f455e 100644
--- a/test/quadratic_form.jl
+++ b/test/quadratic_form.jl
@@ -29,7 +29,7 @@ Base.getindex(g::Gramm, i, j) = g.matrix[i, j]
     ]
     Q = SA.QuadraticForm(Gramm(m, gbasis))
     b = basis(Q)
-    @test A(Q) == π * b[1]' * b[1]
+    @test A(Q) == π * b[1] * b[1]
 
     m = [
         0 1//2 0 0
@@ -39,7 +39,7 @@ Base.getindex(g::Gramm, i, j) = g.matrix[i, j]
     ]
     Q = SA.QuadraticForm(Gramm(m, gbasis))
     b = basis(Q)
-    @test A(Q) == 1 // 2 * b[1]' * b[2]
+    @test A(Q) == 1 // 2 * b[1] * b[2]
 
     m = [
         0 0 0 0
@@ -49,7 +49,7 @@ Base.getindex(g::Gramm, i, j) = g.matrix[i, j]
     ]
     Q = SA.QuadraticForm(Gramm(m, gbasis))
     b = basis(Q)
-    @test A(Q) == π * b[2]' * b[2] - b[2]' * b[3]
+    @test A(Q) == π * b[2] * b[2] - b[2] * b[3]
 
     m = [
         0 0 0 0
@@ -59,9 +59,15 @@ Base.getindex(g::Gramm, i, j) = g.matrix[i, j]
     ]
     Q = SA.QuadraticForm(Gramm(m, gbasis))
     b = basis(Q)
+    @test A(Q) == π * b[3] * b[2] + b[4] * b[3]
+
+    Q = SA.QuadraticForm{SA.star}(Gramm(m, gbasis))
     @test A(Q) == π * b[3]' * b[2] + b[4]' * b[3]
 
     m = ones(Int, 4, 4)
     Q = SA.QuadraticForm(Gramm(m, gbasis))
+    @test A(Q) == sum(bi * bj for bi in gbasis for bj in gbasis)
+
+    Q = SA.QuadraticForm{SA.star}(Gramm(m, gbasis))
     @test A(Q) == sum(bi' * bj for bi in gbasis for bj in gbasis)
 end