Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 58 additions & 56 deletions src/DynamicPolynomials.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,85 +8,87 @@ import MutableArithmetics as MA
import StarAlgebras as SA

include("var.jl")
#const CommutativeVariable{O,M} = Variable{Commutative{O},M}
#const NonCommutativeVariable{O,M} = Variable{NonCommutative{O},M}
include("mono.jl")
const DMonomialLike{V,M} = Union{Monomial{V,M},Variable{V,M}}
MA.mutability(::Type{<:Monomial{<:Commutative}}) = MA.IsMutable()
MA.mutability(::Type{<:Monomial{<:NonCommutative}}) = MA.IsNotMutable()
const _Term{V,M,T} = MP.Term{T,Monomial{V,M}}
function __add_variables!(t::_Term, allvars, map)
return __add_variables!(MP.monomial(t), allvars, map)
end
include("monomial_vector.jl")
include("poly.jl")
MA.mutability(::Type{<:Polynomial}) = MA.IsMutable()
const TermPoly{V,M,T} = Union{_Term{V,M,T},Polynomial{V,M,T}}
const PolyType{V,M} =
Union{Polynomial{V,M},_Term{V,M},Monomial{V,M},Variable{V,M}}

# Monomials are now MP.Polynomial{MP.Monomial, V, E}
# No separate Monomial struct — variables and exponents live in the basis element.

# The monomial type for DP variables
const DMonomialLike{V} = Union{MP.Polynomial{MP.Monomial,V},Variable{V}}
# Convenience alias for the specific monomial type with DP variables
const DPMonomial{V,M} = MP.Polynomial{MP.Monomial,Vector{Variable{V,M}},Vector{Int}}

function MP.variable_union_type(
::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
::Union{Variable{V,M},Type{<:Variable{V,M}}},
) where {V,M}
return Variable{V,M}
end
MP.constant_monomial(::Type{<:PolyType{V,M}}) where {V,M} = Monomial{V,M}()
function MP.constant_monomial(p::PolyType)
return Monomial(copy(MP.variables(p)), zeros(Int, nvariables(p)))
end
MP.monomial_type(::Type{<:PolyType{V,M}}) where {V,M} = Monomial{V,M}
MP.monomial_type(::PolyType{V,M}) where {V,M} = Monomial{V,M}
MP.ordering(p::PolyType) = MP.ordering(MP.variable_union_type(p))
#function MP.constant_monomial(::Type{Monomial{V,M}}, vars=Variable{V,M}[]) where {V,M}
# return Monomial{V,M}(vars, zeros(Int, length(vars)))
#end
function MP.term_type(
::Union{TermPoly{V,M,T},Type{<:TermPoly{V,M,T}}},
) where {V,M,T}
return _Term{V,M,T}
function MP.variable_union_type(
::Union{DPMonomial{V,M},Type{<:DPMonomial{V,M}}},
) where {V,M}
return Variable{V,M}
end

MP.constant_monomial(::Type{DPMonomial{V,M}}) where {V,M} = MP.Polynomial(
MP.Variables{MP.Monomial}(Variable{V,M}[]),
Int[],
)
MP.monomial_type(::Type{<:DPMonomial{V,M}}) where {V,M} = DPMonomial{V,M}
MP.monomial_type(::DPMonomial{V,M}) where {V,M} = DPMonomial{V,M}
MP.monomial_type(::Type{<:Variable{V,M}}) where {V,M} = DPMonomial{V,M}
MP.monomial_type(::Variable{V,M}) where {V,M} = DPMonomial{V,M}
# MP.ordering for Variable is in var.jl

function MP.term_type(
::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
::Union{Variable{V,M},Type{<:Variable{V,M}}},
::Type{T},
) where {V,M,T}
return _Term{V,M,T}
# Create a term via the convenience constructor
# term_type needs to return a constructible type...
# For now just return SA.Term{T} since concrete type depends on algebra
return SA.Term{T}
end
MP.term_type(::Type{Polynomial{V,M}}) where {V,M} = _Term{V,M}
MP.polynomial_type(::Type{_Term{V,M}}) where {V,M} = Polynomial{V,M}
MP.polynomial_type(::Type{_Term{V,M,T}}) where {T,V,M} = Polynomial{V,M,T}
function MP.polynomial_type(
::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
::Type{T},
) where {V,M,T}
return Polynomial{V,M,T}
end
MP.variables(p::AbstractArray{<:PolyType}) = mergevars(MP.variables.(p))[1]
function MP.nvariables(
p::Union{PolyType,MonomialVector,AbstractArray{<:PolyType}},
)

MP.variables(p::AbstractArray{<:Variable}) = mergevars(MP.variables.(p))[1]
function MP.nvariables(p::Union{Variable,AbstractArray{<:Variable}})
return length(MP.variables(p))
end
function MP.similar_variable(
P::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
P::Union{Variable{V,M},Type{<:Variable{V,M}}},
::Type{Val{S}},
) where {V,M,S}
return MP.similar_variable(P, S)
end
function MP.similar_variable(p::PolyType{V,M}, s::Symbol) where {V,M}
function MP.similar_variable(p::Variable{V,M}, s::Symbol) where {V,M}
return Variable(string(s), V, M, isreal(p) ? REAL : COMPLEX)
end
function MP.similar_variable(::Type{<:PolyType{V,M}}, s::Symbol) where {V,M}
return Variable(string(s), V, M, REAL) # we cannot infer this from the type
function MP.similar_variable(::Type{<:Variable{V,M}}, s::Symbol) where {V,M}
return Variable(string(s), V, M, REAL)
end

include("promote.jl")
# Create monomial from variable: Variable → Polynomial{Monomial,...}
function Base.convert(::Type{DPMonomial{V,M}}, x::Variable{V,M}) where {V,M}
return MP.Polynomial{MP.Monomial}(x)
end

# monomial(vars, exps) constructs a Polynomial{Monomial,...}
function MP.monomial(vars::Vector{Variable{V,M}}, z::Vector{Int}) where {V,M}
@assert !iscomm(V) || issorted(vars, rev = true)
return MP.Polynomial(MP.Variables{MP.Monomial}(vars), z)
end

# exponents for variables is in var.jl

include("operators.jl")
include("comp.jl")
include("promote.jl")

include("anti_diff.jl")
include("diff.jl")
include("subs.jl")
# Variable power → monomial
Base.:(^)(x::Variable{V,M}, i::Int) where {V,M} = MP.Polynomial(
MP.Variables{MP.Monomial}([x]),
[i],
)

include("div.jl")
# Variable + Variable → uses term + term → AlgebraElement
Base.:(+)(x::Variable, y::Variable) = MP.term(x) + MP.term(y)
Base.:(-)(x::Variable, y::Variable) = MP.term(x) - MP.term(y)

end # module
82 changes: 1 addition & 81 deletions src/comp.jl
Original file line number Diff line number Diff line change
@@ -1,12 +1,10 @@
import Base.==

# TODO This should be in Base with T instead of Variable{V,M}.
# See https://github.com/blegat/MultivariatePolynomials.jl/issues/3
# Comparison of Variable vectors
function Base.:(==)(x::Vector{Variable{V,M}}, y::Vector{Variable{V,M}}) where {V,M}
if length(x) != length(y)
false
else
#for (xi, yi) in zip(x, y)
for i in 1:length(x)
if x[i] != y[i]
return false
Expand All @@ -17,7 +15,6 @@ function Base.:(==)(x::Vector{Variable{V,M}}, y::Vector{Variable{V,M}}) where {V
end

# Comparison of Variable

const AnyCommutative{O} = Union{Commutative{O},NonCommutative{O}}

function Base.cmp(
Expand All @@ -31,81 +28,4 @@ function Base.cmp(
end
end

# TODO remove
Base.:(==)(x::Variable, y::Variable) = iszero(cmp(x, y))
Base.:(==)(x::Monomial, y::Monomial) = iszero(cmp(x, y))

# Comparison of MonomialVector
function (==)(x::MonomialVector{V,M}, y::MonomialVector{V,M}) where {V,M}
if length(x.Z) != length(y.Z)
return false
end
allvars, maps = mergevars([MP.variables(x), MP.variables(y)])
# Should be sorted in the same order since the non-common
# polyvar should have exponent 0
for (a, b) in zip(x.Z, y.Z)
A = zeros(Int, length(allvars))
B = zeros(Int, length(allvars))
A[maps[1]] = a
B[maps[2]] = b
if A != B
return false
end
end
return true
end
(==)(mv::AbstractVector, x::MonomialVector) = monomial_vector(mv) == x
(==)(x::MonomialVector, mv::AbstractVector) = x == monomial_vector(mv)

# Comparison of Term
function _compare(p::Polynomial{V,M}, q::Polynomial{V,M}, comparator) where {V,M}
# terms should be sorted and without zeros
if MP.nterms(p) != MP.nterms(q)
return false
end
for i in eachindex(p.a)
if !comparator(p.x[i], q.x[i])
# There should not be zero terms
@assert p.a[i] != 0
@assert q.a[i] != 0
return false
end
if !comparator(p.a[i], q.a[i])
return false
end
end
return true
end

(==)(p::Polynomial{V, M}, q::Polynomial{V, M}) where {V, M} = _compare(p, q, (==))
Base.isequal(p::Polynomial{V, M}, q::Polynomial{V, M}) where {V, M} = _compare(p, q, isequal)

function Base.isapprox(
p::Polynomial{V,M,S},
q::Polynomial{V,M,T};
rtol::Real = Base.rtoldefault(S, T, 0),
atol::Real = 0,
ztol::Real = iszero(atol) ? Base.rtoldefault(S, T, 0) : atol,
) where {V,M,S,T}
i = j = 1
while i <= length(p.x) || j <= length(q.x)
if i > length(p.x) || (j <= length(q.x) && q.x[j] < p.x[i])
if !isapproxzero(q.a[j], ztol = ztol)
return false
end
j += 1
elseif j > length(q.x) || p.x[i] < q.x[j]
if !isapproxzero(p.a[i], ztol = ztol)
return false
end
i += 1
else
if !isapprox(p.a[i], q.a[j], rtol = rtol, atol = atol)
return false
end
i += 1
j += 1
end
end
return true
end
18 changes: 0 additions & 18 deletions src/diff.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,21 +8,3 @@ function MP.differentiate(m::Monomial{V,M}, x::Variable{V,M}) where {V,M}
m.z[i] * Monomial(MP.variables(m), z)
end
end

function MP.differentiate(p::Polynomial{V,M,T}, x::Variable{V,M}) where {V,M,T}
# grlex order preserved
i = something(findfirst(isequal(x), MP.variables(p)), 0)
S = typeof(zero(T) * 0)
if iszero(i)
zero(Polynomial{V,M,S})
else
keep = findall(z -> z[i] > 0, p.x.Z)
Z = copy.(p.x.Z[keep])
a = Vector{S}(undef, length(keep))
for j in 1:length(Z)
a[j] = p.a[keep[j]] * Z[j][i]
Z[j][i] -= 1
end
Polynomial(a, MonomialVector(MP.variables(p), Z))
end
end
Loading
Loading