diff --git a/src/StarAlgebras.jl b/src/StarAlgebras.jl index c8c8757..edfbcd5 100644 --- a/src/StarAlgebras.jl +++ b/src/StarAlgebras.jl @@ -8,8 +8,8 @@ import LinearAlgebra import MutableArithmetics as MA -export StarAlgebra, AlgebraElement -export basis, coeffs, star +export StarAlgebra, AlgebraElement, Term +export basis, coeffs, star, coefficient, basis_element function star end @@ -32,6 +32,7 @@ include("mtables.jl") # Algebras and elts include("types.jl") +include("term.jl") include("algebra_elts.jl") include("star.jl") diff --git a/src/term.jl b/src/term.jl new file mode 100644 index 0000000..9d2a791 --- /dev/null +++ b/src/term.jl @@ -0,0 +1,139 @@ +# This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE +# Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat + +""" + struct Term{T,B} <: MA.AbstractMutable + coefficient::T + basis_element::B + end + +A representation of the multiplication between a `coefficient` and a `basis_element`. +This is the generic analog of a single term in an algebra: a scalar times a basis element. + +The `coefficient` does not need to be a `Number`. For instance, in multivariate +polynomial GCD computations, the coefficient can itself be a polynomial. +""" +struct Term{T,B} <: MA.AbstractMutable + coefficient::T + basis_element::B +end + +# Copy constructor (needed e.g. for array operations) +Term{T,B}(t::Term{T,B}) where {T,B} = t + +""" + coefficient(t::Term) + +Return the coefficient of the term `t`. +""" +coefficient(t::Term) = t.coefficient + +""" + basis_element(t::Term) + +Return the basis element of the term `t`. +""" +basis_element(t::Term) = t.basis_element + +Base.iszero(t::Term) = iszero(coefficient(t)) +Base.isone(t::Term) = isone(coefficient(t)) && isone(basis_element(t)) + +Base.zero(t::Term) = Term(zero(coefficient(t)), basis_element(t)) +Base.one(t::Term) = Term(one(coefficient(t)), one(basis_element(t))) +Base.one(::Type{Term{T,B}}) where {T,B} = Term(one(T), one(B)) + +function Base.:(==)(t1::Term, t2::Term) + c1 = coefficient(t1) + c2 = coefficient(t2) + if iszero(c1) + return iszero(c2) + else + return c1 == c2 && basis_element(t1) == basis_element(t2) + end +end +function Base.isequal(t1::Term, t2::Term) + c1 = coefficient(t1) + c2 = coefficient(t2) + if iszero(c1) + return iszero(c2) + else + return isequal(c1, c2) && isequal(basis_element(t1), basis_element(t2)) + end +end + +function Base.hash(t::Term, u::UInt) + if iszero(t) + return hash(0, u) + elseif isone(coefficient(t)) + return hash(basis_element(t), u) + else + return hash(basis_element(t), hash(coefficient(t), u)) + end +end + +function Base.convert(::Type{Term{T,B}}, t::Term) where {T,B} + return Term{T,B}(convert(T, coefficient(t)), convert(B, basis_element(t))) +end +Base.convert(::Type{Term{T,B}}, t::Term{T,B}) where {T,B} = t + +function Base.promote_rule( + ::Type{Term{S,B1}}, + ::Type{Term{T,B2}}, +) where {S,T,B1,B2} + return Term{promote_type(S, T),promote_type(B1, B2)} +end + +Base.:^(x::Term, p::Integer) = Base.power_by_squaring(x, p) + +Base.ndims(::Union{Type{<:Term},Term}) = 0 +Base.broadcastable(t::Term) = Ref(t) + +# `Base.power_by_squaring` calls `Base.copy` and we want +# `t^1` to be a mutable copy of `t` so `copy` needs to be +# the same as `mutable_copy`. +Base.copy(t::Term) = MA.mutable_copy(t) +function MA.mutable_copy(t::Term) + return Term( + MA.copy_if_mutable(coefficient(t)), + MA.copy_if_mutable(basis_element(t)), + ) +end + +function MA.mutability(::Type{Term{T,B}}) where {T,B} + if MA.mutability(T) isa MA.IsMutable && MA.mutability(B) isa MA.IsMutable + return MA.IsMutable() + else + return MA.IsNotMutable() + end +end + +# dot for Term to avoid recursive fallback in LinearAlgebra.dot +function LinearAlgebra.dot(t1::Term, t2::Term) + return coefficient(t1) * + coefficient(t2) * + (basis_element(t1) * basis_element(t2)) +end +function LinearAlgebra.dot(x, t::Term) + return x * t +end +function LinearAlgebra.dot(t::Term, x) + return star(t) * x +end + +function MA.operate_to!(t::Term, ::typeof(*), t1::Term, t2::Term) + MA.operate_to!(t.coefficient, *, coefficient(t1), coefficient(t2)) + MA.operate_to!(t.basis_element, *, basis_element(t1), basis_element(t2)) + return t +end + +function MA.operate!(::typeof(*), t1::Term, t2::Term) + MA.operate!(*, t1.coefficient, coefficient(t2)) + MA.operate!(*, t1.basis_element, basis_element(t2)) + return t1 +end + +function MA.operate!(::typeof(one), t::Term) + MA.operate!(one, t.coefficient) + MA.operate!(one, t.basis_element) + return t +end diff --git a/test/runtests.jl b/test/runtests.jl index 77d4dfd..209ff07 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,6 +46,7 @@ function test_vector_interface(basis, vector = collect(basis)) end @testset "StarAlgebras" begin + include("term.jl") include("basic.jl") # proof of concept using PermutationGroups diff --git a/test/term.jl b/test/term.jl new file mode 100644 index 0000000..18dd895 --- /dev/null +++ b/test/term.jl @@ -0,0 +1,246 @@ +# This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE +# Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat + +# Mutable wrapper for testing MA operations +mutable struct MutableInt <: Number + value::Int +end +Base.zero(::MutableInt) = MutableInt(0) +Base.zero(::Type{MutableInt}) = MutableInt(0) +Base.one(::MutableInt) = MutableInt(1) +Base.one(::Type{MutableInt}) = MutableInt(1) +Base.isone(x::MutableInt) = isone(x.value) +Base.iszero(x::MutableInt) = iszero(x.value) +Base.:(==)(a::MutableInt, b::MutableInt) = a.value == b.value +Base.:*(a::MutableInt, b::MutableInt) = MutableInt(a.value * b.value) +MA.mutability(::Type{MutableInt}) = MA.IsMutable() +function MA.operate_to!( + res::MutableInt, + ::typeof(*), + a::MutableInt, + b::MutableInt, +) + res.value = a.value * b.value + return res +end +function MA.operate!(::typeof(*), a::MutableInt, b::MutableInt) + a.value *= b.value + return a +end +function MA.operate!(::typeof(one), a::MutableInt) + a.value = 1 + return a +end +function MA.copy_if_mutable(a::MutableInt) + return MutableInt(a.value) +end +function MA.mutable_copy(a::MutableInt) + return MutableInt(a.value) +end + +# Minimal arithmetic on Term{Int,Monomial} using the bivariate example +# to test ^, dot, etc. within SA (no downstream needed). +function Base.:*(t1::Term{Int,Monomial}, t2::Term{Int,Monomial}) + return Term( + coefficient(t1) * coefficient(t2), + basis_element(t1) * basis_element(t2), + ) +end +function Base.:*(α::Int, t::Term{Int,Monomial}) + return Term(α * coefficient(t), basis_element(t)) +end +function Base.:*(t::Term{Int,Monomial}, α::Int) + return Term(coefficient(t) * α, basis_element(t)) +end +function Base.:*(α::Int, m::Monomial) + return Term(α, m) +end +SA.star(m::Monomial) = m +SA.star(t::Term{Int,Monomial}) = Term(coefficient(t), star(basis_element(t))) +Base.conj(t::Term{Int,Monomial}) = star(t) + +@testset "Term" begin + @testset "Accessors" begin + t = Term(3.0, Monomial((1, 2))) + @test coefficient(t) == 3.0 + @test basis_element(t) == Monomial((1, 2)) + end + + @testset "iszero / zero" begin + t = Term(3, Monomial((1, 0))) + @test !iszero(t) + z = zero(t) + @test iszero(z) + @test coefficient(z) == 0 + @test basis_element(z) == Monomial((1, 0)) + + @test iszero(Term(0.0, Monomial((0, 0)))) + @test !iszero(Term(1, Monomial((0, 0)))) + end + + @testset "copy / mutable_copy" begin + m = Monomial((2, 3)) + t = Term(5, m) + t2 = copy(t) + @test coefficient(t2) == 5 + @test basis_element(t2) == m + + t3 = MA.mutable_copy(t) + @test coefficient(t3) == coefficient(t) + @test basis_element(t3) == basis_element(t) + + # Test that copy of mutable fields produces independent copies + a = MutableInt(7) + b = MutableInt(3) + t_mut = Term(a, b) + t_copy = copy(t_mut) + @test coefficient(t_copy) == MutableInt(7) + @test basis_element(t_copy) == MutableInt(3) + # Mutating the copy should not affect the original + t_copy.coefficient.value = 99 + @test coefficient(t_mut) == MutableInt(7) + end + + @testset "copy constructor" begin + t = Term(3, Monomial((1, 0))) + T = typeof(t) + t2 = T(t) + @test t2 === t + end + + @testset "mutability" begin + # Int and Monomial (immutable struct) are both not mutable + @test MA.mutability(Term{Int,Monomial}) isa MA.IsNotMutable + # MutableInt on both sides gives mutable + @test MA.mutability(Term{MutableInt,MutableInt}) isa MA.IsMutable + # Mixed: one immutable makes the term not mutable + @test MA.mutability(Term{MutableInt,Int}) isa MA.IsNotMutable + @test MA.mutability(Term{Int,MutableInt}) isa MA.IsNotMutable + end + + @testset "operate_to!(*, ...)" begin + t1 = Term(MutableInt(3), MutableInt(2)) + t2 = Term(MutableInt(4), MutableInt(5)) + res = Term(MutableInt(0), MutableInt(0)) + MA.operate_to!(res, *, t1, t2) + @test coefficient(res) == MutableInt(12) + @test basis_element(res) == MutableInt(10) + # Original terms are unchanged + @test coefficient(t1) == MutableInt(3) + @test coefficient(t2) == MutableInt(4) + end + + @testset "operate!(*, ...)" begin + t1 = Term(MutableInt(3), MutableInt(2)) + t2 = Term(MutableInt(4), MutableInt(5)) + MA.operate!(*, t1, t2) + @test coefficient(t1) == MutableInt(12) + @test basis_element(t1) == MutableInt(10) + end + + @testset "operate!(one, ...)" begin + t = Term(MutableInt(5), MutableInt(7)) + MA.operate!(one, t) + @test coefficient(t) == MutableInt(1) + @test basis_element(t) == MutableInt(1) + end + + @testset "one / isone" begin + t = Term(3, Monomial((1, 0))) + @test !isone(t) + t1 = Term(1, Monomial((0, 0))) + @test isone(t1) + @test one(t1) == t1 + # one from type (needs one(::Type{B}) for the basis element) + o = one(Term{Int,MutableInt}) + @test isone(o) + @test coefficient(o) == 1 + @test basis_element(o) == MutableInt(1) + end + + @testset "== / isequal" begin + t1 = Term(3.0, Monomial((1, 0))) + t2 = Term(3.0, Monomial((1, 0))) + t3 = Term(4.0, Monomial((1, 0))) + t4 = Term(3.0, Monomial((0, 1))) + @test t1 == t2 + @test t1 != t3 + @test t1 != t4 + @test isequal(t1, t2) + @test !isequal(t1, t3) + # zero terms are equal regardless of basis element + @test Term(0.0, Monomial((1, 0))) == Term(0.0, Monomial((0, 1))) + # isequal for zero terms + @test isequal(Term(0.0, Monomial((1, 0))), Term(0.0, Monomial((0, 1)))) + end + + @testset "hash" begin + t1 = Term(3, Monomial((1, 0))) + t2 = Term(3, Monomial((1, 0))) + @test hash(t1) == hash(t2) + # hash of zero term + @test hash(Term(0, Monomial((1, 0)))) == hash(0) + # hash of term with coefficient 1 matches hash of basis_element + t3 = Term(1, Monomial((2, 0))) + @test hash(t3) == hash(Monomial((2, 0))) + end + + @testset "convert / promote_rule" begin + t = Term(3, Monomial((1, 0))) + t2 = convert(Term{Float64,Monomial}, t) + @test coefficient(t2) == 3.0 + @test basis_element(t2) == Monomial((1, 0)) + @test promote_type(Term{Int,Monomial}, Term{Float64,Monomial}) == + Term{Float64,Monomial} + # identity convert + t3 = convert(typeof(t), t) + @test t3 === t + end + + @testset "^ (power)" begin + t = Term(2, Monomial((1, 0))) + @test t^1 == t + t2 = t^2 + @test coefficient(t2) == 4 + @test basis_element(t2) == Monomial((2, 0)) + @test t^0 == one(t) + end + + @testset "broadcastable / ndims" begin + t = Term(2, Monomial((1, 0))) + @test ndims(t) == 0 + @test ndims(typeof(t)) == 0 + @test Base.broadcastable(t) isa Ref + end + + @testset "dot" begin + # dot(::Term, ::Term) computes coeff*coeff * (basis*basis) + t1 = Term(2, Monomial((1, 0))) + t2 = Term(3, Monomial((0, 1))) + d = LinearAlgebra.dot(t1, t2) + @test d == 6 * Monomial((1, 1)) + # dot(scalar, Term) + @test LinearAlgebra.dot(2, t1) == Term(4, Monomial((1, 0))) + # dot(Term, scalar) uses star(t) * x + @test LinearAlgebra.dot(t1, 3) == Term(6, Monomial((1, 0))) + end + + @testset "various coefficient types" begin + # Float64 coefficient + t = Term(2.5, Monomial((1, 1))) + @test coefficient(t) == 2.5 + @test !iszero(t) + + # Complex coefficient + tc = Term(1 + 2im, Monomial((0, 1))) + @test coefficient(tc) == 1 + 2im + @test !iszero(tc) + @test iszero(zero(tc)) + @test coefficient(zero(tc)) == 0 + 0im + + # Rational coefficient + tr = Term(3 // 4, Monomial((2, 0))) + @test coefficient(tr) == 3 // 4 + @test coefficient(zero(tr)) == 0 // 1 + end +end