diff --git a/docs/src/types.md b/docs/src/types.md index e80d42a5..08302dee 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -14,6 +14,14 @@ name_base_indices variable_union_type similar_variable @similar_variable +conj(::AbstractVariable) +real(::AbstractVariable) +imag(::AbstractVariable) +isreal(::AbstractVariable) +isrealpart +isimagpart +isconj +ordinary_variable ``` ## Monomials @@ -47,6 +55,8 @@ coefficient_type monomial constant_term zero_term +degree_complex +halfdegree ``` ## Polynomials @@ -75,6 +85,16 @@ monic map_coefficients map_coefficients! map_coefficients_to! +conj(::_APL) +real(::_APL) +imag(::_APL) +isreal(::_APL) +mindegree_complex +minhalfdegree +maxdegree_complex +maxhalfdegree +extdegree_complex +exthalfdegree ``` ## Rational Polynomial Function @@ -90,4 +110,8 @@ monomial_vector_type empty_monomial_vector sort_monomial_vector merge_monomial_vectors +conj(::AbstractVector{<:AbstractMonomial}) +real(::AbstractVector{<:AbstractMonomial}) +imag(::AbstractVector{<:AbstractMonomial}) +isreal(::AbstractVector{<:AbstractMonomial}) ``` diff --git a/src/MultivariatePolynomials.jl b/src/MultivariatePolynomials.jl index 303997fa..f3bc5f34 100644 --- a/src/MultivariatePolynomials.jl +++ b/src/MultivariatePolynomials.jl @@ -77,6 +77,7 @@ include("hash.jl") include("promote.jl") include("conversion.jl") +include("complex.jl") include("operators.jl") include("comparison.jl") diff --git a/src/complex.jl b/src/complex.jl new file mode 100644 index 00000000..7a3f4e1f --- /dev/null +++ b/src/complex.jl @@ -0,0 +1,381 @@ +export isrealpart, + isimagpart, + isconj, + ordinary_variable, + degree_complex, + halfdegree, + mindegree_complex, + minhalfdegree, + maxdegree_complex, + maxhalfdegree, + extdegree_complex, + exthalfdegree + +""" + isreal(x::AbstractVariable) + +Return whether a given variable was declared as a real-valued or complex-valued variable (also their conjugates are complex, +but their real and imaginary parts are not). +By default, all variables are real-valued. +""" +Base.isreal(::AbstractVariable) = true + +""" + isrealpart(x::AbstractVariable) + +Return whether the given variable is the real part of a complex-valued variable. + +See also [`isreal`](@ref), [`isimagpart`](@ref), [`isconj`](@ref). +""" +isrealpart(::AbstractVariable) = false + +""" + isimagpart(x::AbstractVariable) + +Return whether the given variable is the imaginary part of a complex-valued variable. + +See also [`isreal`](@ref), [`isrealpart`](@ref), [`isconj`](@ref). +""" +isimagpart(::AbstractVariable) = false + +""" + isconj(x::AbstractVariable) + +Return whether the given variable is obtained by conjugating a user-defined complex-valued variable. + +See also [`isreal`](@ref), [`isrealpart`](@ref), [`isimagpart`](@ref). +""" +isconj(::AbstractVariable) = false + +""" + ordinary_variable(x::Union{AbstractVariable, AbstractVector{<:AbstractVariable}}) + +Given some (complex-valued) variable that was transformed by conjugation, taking its real part, or taking its +imaginary part, return the original variable as it was defined by the user. + +See also [`conj`](@ref), [`real`](@ref), [`imag`](@ref). +""" +ordinary_variable(x::AbstractVariable) = MA.copy_if_mutable(x) + +""" + conj(x::AbstractVariable) + +Return the complex conjugate of a given variable if it was declared as a complex variable; else return the +variable unchanged. + +See also [`isreal`](@ref), [`isconj`](@ref). +""" +Base.conj(x::AbstractVariable) = MA.copy_if_mutable(x) +@doc """ + conj(x::AbstractPolynomialLike) + +Return the complex conjugate of `x` by applying conjugation to all coefficients and variables. +""" conj(::_APL) +@doc """ + conj(x::AbstractVector{<:AbstractMonomial}) + +Return the complex conjugate of `x` by applying conjugation to monomials. +""" conj(::AbstractVector{<:AbstractMonomial}) + +""" + real(x::AbstractVariable) + +Return the real part of a given variable if it was declared as a complex variable; else return the variable +unchanged. + +See also [`imag`](@ref). +""" +Base.real(x::AbstractVariable) = MA.copy_if_mutable(x) +@doc """ + real(x::AbstractPolynomialLike) + +Return the real part of `x` by applying [`real`](@ref) to all coefficients and variables; for this purpose, every +complex-valued variable is decomposed into its real- and imaginary parts. + +See also [`imag`](@ref). +""" real(::_APL) +@doc """ + real(x::AbstractVector{<:AbstractMonomial}) + +Return the real part of `x` by applying [`real`](@ref) to all monomials; for this purpose, every complex-valued variable is +decomposed into its real- and imaginary parts. Note that the result will no longer be a monomial vector. + +See also [`imag`](@ref). +""" real(::AbstractVector{<:AbstractMonomial}) + +""" + imag(x::AbstractVariable) + +Return the imaginary part of a given variable if it was declared as a complex variable; else return zero. + +See also [`isreal`](@ref), [`isimagpart`](@ref), [`real`](@ref). +""" +Base.imag(::AbstractVariable) = MA.Zero() +@doc """ + imag(x::AbstractPolynomialLike) + +Return the imaginary part of `x` by applying [`imag`](@ref) to all coefficients and variables; for this purpose, every +complex-valued variable is decomposed into its real- and imaginary parts. + +See also [`real`](@ref). +""" imag(::_APL) +@doc """ + imag(x::AbstractVector{<:AbstractMonomial}) + +Return the imaginary part of `x` by applying [`imag`](@ref) to all monomials; for this purpose, every complex-valued variable +is decomposed into its real- and imaginary parts. Note that the result will no longer be a monomial vector. + +See also [`real`](@ref). +""" imag(::AbstractVector{<:AbstractMonomial}) + +# extend to higher-level elements. We make all those type-stable (but we need convert, as the construction method may +# deliver simpler types than the inputs if they were deliberately casted, e.g., term to monomial) +""" + isreal(p::AbstractPolynomialLike) + +Returns `true` if and only if no single variable in `p` was declared as a complex variable (in the sense that [`isreal`](@ref) +applied on them would be `true`) and no coefficient is complex-valued. +""" +function Base.isreal(p::_APL) + for v in variables(p) + if !isreal(v) && !iszero(maxdegree(p, v)) + return false + end + end + return all(isreal, coefficients(p)) +end +""" + isreal(p::AbstractVector{<:AbstractMonomial}) + +Returns `true` if and only if every single monomial in `p` would is real-valued. +""" +Base.isreal(p::AbstractVector{<:AbstractMonomial}) = all(isreal, p) + +function Base.conj(x::M) where {M<:AbstractMonomial} + return isreal(x) ? x : + convert( + M, + reduce( + *, + conj(var)^exp for (var, exp) in powers(x); + init = constant_monomial(x), + ), + ) +end +function Base.conj(x::V) where {V<:AbstractVector{<:AbstractMonomial}} + return isreal(x) ? MA.copy_if_mutable(x) : monomial_vector(conj.(x)) +end +function Base.conj(x::T) where {T<:AbstractTerm} + return isreal(x) ? MA.copy_if_mutable(x) : + convert(T, conj(coefficient(x)) * conj(monomial(x))) +end +function Base.conj(x::P) where {P<:AbstractPolynomial} + return iszero(x) || isreal(x) ? MA.copy_if_mutable(x) : + convert(P, polynomial([conj(t) for t in x])) +end + +# Real and imaginary parts are harder to realize. The real part of a monomial can easily be a polynomial. +for fun in [:real, :imag] + eval( + quote + function Base.$fun(x::_APL) + # Note: x may also be a polynomial; we could handle this case separately by performing the replacement in each + # individual term, then adding them all up. This could potentially lower the overall memory requirement (in case + # the expansions of the individual terms simplify) at the expense of not being able to exploit optimizations that + # subst can do by knowing the full polynomial. + iszero(x) && + return zero(polynomial_type(x, real(coefficient_type(x)))) + # We replace every complex variable by its decomposition into real and imaginary part + subst_vars = filter(Base.:! ∘ isreal, variables(x)) + # To avoid a stack overflow on promote_type, we'll handle the empty case separately + full_version = + isempty(subst_vars) ? polynomial(x) : + subs( + x, + subst_vars => + [real(var) + 1im * imag(var) for var in subst_vars], + ) + # Now everything that is imaginary can be recognized by its coefficient + return convert( + polynomial_type( + full_version, + real(coefficient_type(full_version)), + ), + map_coefficients!($fun, full_version), + ) + end + function Base.$fun(x::AbstractVector{<:AbstractMonomial}) + return map(Base.$fun, x) + end + end, + ) +end + +# Also give complex-valued degree definitions. We choose not to overwrite degree, as this will lead to issues in monovecs +# and their sorting. So now there are two ways to calculate degrees: strictly by considering all variables independently, +# and also by looking at their complex structure. +""" + degree_complex(t::AbstractTermLike) + +Return the _total complex degree_ of the monomial of the term `t`, i.e., the maximum of the total degree of the declared +variables in `t` and the total degree of the conjugate variables in `t`. +To be well-defined, the monomial must not contain real parts or imaginary parts of variables. + + degree_complex(t::AbstractTermLike, v::AbstractVariable) + +Returns the exponent of the variable `v` or its conjugate in the monomial of the term `t`, whatever is larger. + +See also [`isconj`](@ref). +""" +function degree_complex(t::AbstractTermLike) + vars = variables(t) + @assert(!any(isrealpart, vars) && !any(isimagpart, vars)) + grouping = isconj.(vars) + exps = exponents(t) + return max(sum(exps[grouping]), sum(exps[map(!, grouping)])) +end +function degree_complex(t::AbstractTermLike, var::AbstractVariable) + return degree_complex(monomial(t), var) +end +degree_complex(v::AbstractVariable, var::AbstractVariable) = (v == var ? 1 : 0) +function degree_complex(m::AbstractMonomial, v::AbstractVariable) + deg = 0 + deg_c = 0 + c_v = conj(v) + for (var, exp) in powers(m) + @assert(!isrealpart(var) && !isimagpart(var)) + if var == v + deg += exp + elseif var == c_v + deg_c += exp + end + end + return max(deg, deg_c) +end + +""" + halfdegree(t::AbstractTermLike) + +Return the equivalent of `ceil(degree(t)/2)`` for real-valued terms or `degree_complex(t)` for terms with only complex +variables; however, respect any mixing between complex and real-valued variables. +""" +function halfdegree(t::AbstractTermLike) + realdeg = 0 + cpdeg = 0 + conjdeg = 0 + for (var, exp) in powers(t) + if isreal(var) + realdeg += exp + else + if isconj(var) + conjdeg += exp + else + @assert(!isrealpart(var) && !isimagpart(var)) + cpdeg += exp + end + end + end + return ((realdeg + 1) >> 1) + max(cpdeg, conjdeg) +end + +""" + mindegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Return the minimal total complex degree of the monomials of `p`, i.e., `minimum(degree_complex, terms(p))`. + + mindegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable) + +Return the minimal complex degree of the monomials of `p` in the variable `v`, i.e., `minimum(degree_complex.(terms(p), v))`. +""" +function mindegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) + return isempty(X) ? 0 : + minimum(t -> degree_complex(t, args...), X, init = 0) +end +function mindegree_complex(p::_APL, args...) + return mindegree_complex(terms(p), args...) +end + +""" + minhalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Return the minmal half degree of the monomials of `p`, i.e., `minimum(halfdegree, terms(p))` +""" +function minhalfdegree(X::AbstractVector{<:AbstractTermLike}, args...) + return isempty(X) ? 0 : minimum(halfdegree, X, init = 0) +end +minhalfdegree(p::_APL) = minhalfdegree(terms(p)) + +""" + maxdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Return the maximal total complex degree of the monomials of `p`, i.e., `maximum(degree_complex, terms(p))`. + + maxdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable) + +Return the maximal complex degree of the monomials of `p` in the variable `v`, i.e., `maximum(degree_complex.(terms(p), v))`. +""" +function maxdegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) + return mapreduce(t -> degree_complex(t, args...), max, X, init = 0) +end +function maxdegree_complex(p::_APL, args...) + return maxdegree_complex(terms(p), args...) +end + +""" + maxhalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Return the maximal half degree of the monomials of `p`, i.e., `maximum(halfdegree, terms(p))` +""" +function maxhalfdegree(X::AbstractVector{<:AbstractTermLike}, args...) + return isempty(X) ? 0 : maximum(halfdegree, X, init = 0) +end +maxhalfdegree(p::_APL) = maxhalfdegree(terms(p)) + +""" + extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Returns the extremal total complex degrees of the monomials of `p`, i.e., `(mindegree_complex(p), maxdegree_complex(p))`. + + extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable) + +Returns the extremal complex degrees of the monomials of `p` in the variable `v`, i.e., +`(mindegree_complex(p, v), maxdegree_complex(p, v))`. +""" +function extdegree_complex( + p::Union{_APL,AbstractVector{<:AbstractTermLike}}, + args..., +) + return (mindegree_complex(p, args...), maxdegree_complex(p, args...)) +end + +""" + exthalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}) + +Return the extremal half degree of the monomials of `p`, i.e., `(minhalfdegree(p), maxhalfdegree(p))` +""" +function exthalfdegree(p::Union{_APL,AbstractVector{<:AbstractTermLike}}) + return (minhalfdegree(p), maxhalfdegree(p)) +end + +function ordinary_variable(x::AbstractVector{<:AbstractVariable}) + # let's assume the number of elements in x is small, else a conversion to a dict (probably better OrderedDict) would + # be better + results = similar(x, 0) + sizehint!(results, length(x)) + j = 0 + for el in x + ov = ordinary_variable(el) + found = false + for i in 1:j + if results[i] == ov + found = true + break + end + end + if !found + push!(results, ov) + j += 1 + end + end + return results +end diff --git a/src/default_polynomial.jl b/src/default_polynomial.jl index 5cc6fb20..a7b45bf6 100644 --- a/src/default_polynomial.jl +++ b/src/default_polynomial.jl @@ -213,8 +213,6 @@ function Base.:(-)( return map_coefficients(J -> J.λ, p1, nonzero = true) - p2 end -LinearAlgebra.adjoint(x::Polynomial) = polynomial!(adjoint.(terms(x))) - function map_coefficients( f::F, p::Polynomial; diff --git a/src/default_term.jl b/src/default_term.jl index 92c9664a..34509bbc 100644 --- a/src/default_term.jl +++ b/src/default_term.jl @@ -31,8 +31,6 @@ end (t::Term)(s...) = substitute(Eval(), t, s) -LinearAlgebra.adjoint(t::Term) = Term(adjoint(coefficient(t)), monomial(t)) - function Base.convert(::Type{Term{T,M}}, m::AbstractMonomialLike) where {T,M} return Term(one(T), convert(M, m)) end diff --git a/src/operators.jl b/src/operators.jl index aee6b139..22dbf3d9 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -446,10 +446,10 @@ function _promote_terms( return promote(t1, t2) end -LinearAlgebra.adjoint(v::AbstractVariable) = v -LinearAlgebra.adjoint(m::AbstractMonomial) = m +LinearAlgebra.adjoint(v::AbstractVariable) = conj(v) +LinearAlgebra.adjoint(m::AbstractMonomial) = conj(m) function LinearAlgebra.adjoint(t::AbstractTerm) - return _term(LinearAlgebra.adjoint(coefficient(t)), monomial(t)) + return _term(adjoint(coefficient(t)), adjoint(monomial(t))) end function LinearAlgebra.adjoint(p::AbstractPolynomialLike) return polynomial(map(LinearAlgebra.adjoint, terms(p))) @@ -457,6 +457,13 @@ end function LinearAlgebra.adjoint(r::RationalPoly) return adjoint(numerator(r)) / adjoint(denominator(r)) end +LinearAlgebra.hermitian_type(::Type{T}) where {T<:AbstractPolynomialLike} = T +function LinearAlgebra.hermitian(v::AbstractPolynomialLike, ::Symbol) + isreal(v) || error( + "Complex-valued polynomials cannot be interpreted as hermitian scalars", + ) + return v +end LinearAlgebra.transpose(v::AbstractVariable) = v LinearAlgebra.transpose(m::AbstractMonomial) = m diff --git a/src/show.jl b/src/show.jl index 92df2050..74b5901a 100644 --- a/src/show.jl +++ b/src/show.jl @@ -27,12 +27,18 @@ Base.show(io::IO, p::TypesWithShow) = show(io, MIME"text/plain"(), p) # VARIABLES function _show(io::IO, mime::MIME, var::AbstractVariable) base, indices = name_base_indices(var) - if isempty(indices) - print(io, base) + if isconj(var) + for c in base + print(io, c, '\u0305') # displays as overbar (̄z) + end else print(io, base) + end + if !isempty(indices) print_subscript(io, mime, indices) end + isrealpart(var) && print(io, "ᵣ") # suffixes for real and imaginary part: zᵣ, zᵢ + return isimagpart(var) && print(io, "ᵢ") end function print_subscript(io::IO, ::MIME"text/print", index) diff --git a/test/commutativetests.jl b/test/commutativetests.jl index 74775473..adc29c22 100644 --- a/test/commutativetests.jl +++ b/test/commutativetests.jl @@ -9,6 +9,7 @@ include("polynomial.jl") include("det.jl") include("rational.jl") +isdefined(Mod, Symbol("@polycvar")) && include("complex.jl") include("promote.jl") include("hash.jl") diff --git a/test/complex.jl b/test/complex.jl new file mode 100644 index 00000000..709362e7 --- /dev/null +++ b/test/complex.jl @@ -0,0 +1,69 @@ +@testset "ComplexPoly" begin + Mod.@polyvar a + @test isreal(a) + @test !isrealpart(a) + @test !isimagpart(a) + @test !isconj(a) + @test ordinary_variable(a) == a + @test conj(a) == real(a) == a + @test_broken iszero(imag(a)) # no iszero for MA.Zero + @test isreal(a^3 + 5a^2 + 4a) + @test !isreal(a^3 + 5im * a^2 + 4a) + + Mod.@polycvar x y + @test !isreal(x) + @test !isrealpart(x) && + !isrealpart(conj(x)) && + isrealpart(real(x)) && + !isrealpart(imag(x)) + @test !isimagpart(x) && + !isimagpart(conj(x)) && + !isimagpart(real(x)) && + isimagpart(imag(x)) + @test !isconj(x) && isconj(conj(x)) && !isconj(real(x)) && !isconj(imag(x)) + @test ordinary_variable(conj(x)) == + ordinary_variable(real(x)) == + ordinary_variable(imag(x)) == + ordinary_variable(x) + @test conj(x) != x && conj(conj(x)) == x + @test real(x) == real(conj(x)) + @test imag(conj(x)) == -imag(x) + @test real(x * conj(x)) == real(x)^2 + imag(x)^2 + @test real(x + conj(x)) == 2real(x) + @test iszero(imag(x + conj(x))) + @test iszero(real(x - conj(x))) + @test imag(x - conj(x)) == 2imag(x) + @test !isreal(x^3 + 5x^2 + 4x) + @test isreal(a^3 + 5a^2 + 4a + x^0) + @test conj(7x^3 + 6y^2 - (3 + im)x + 8a) == + 7conj(x)^3 + 6conj(y)^2 - (3 - im) * conj(x) + 8a + @test conj(monomial_vector([x, y, x * y^3, y * x^4])) == + monomial_vector([conj(x), conj(y), conj(x * y^3), conj(y * x^4)]) + @test real(x * y^2 * a) == + ( + real(x) * real(y)^2 - imag(x) * 2 * real(y) * imag(y) - + real(x) * imag(y)^2 + ) * a + @test imag(x * y^2 * a) == + ( + imag(x) * real(y)^2 + real(x) * 2 * real(y) * imag(y) - + imag(x) * imag(y)^2 + ) * a + herm = [a x+y^2; conj(x)+conj(y)^2 x^2+conj(x)^2] + @test herm == herm' + + @test degree_complex(x * y^2 * conj(y)^3) == 3 + @test degree_complex(x * y^2 * conj(y)^3, x) == 1 + @test degree_complex(x * y^2 * conj(y)^3, y) == 3 + @test halfdegree(x * y^2 * conj(y^3)) == 3 + @test halfdegree(x * a^5 * conj(y)) == 4 + @test ordinary_variable([x, y, conj(x), a, real(x), imag(y)]) == [x, y, a] + + @test subs(4x + 8y^2 - 6x^3, [x, y] => [2 + 4im, 9 - im]) == + (1176 - 32im) * x^0 + @test subs(4x + 8y^2 - 6x^3, [x, conj(y)] => [2 + 4im, 9 - im]) == + (1176 + 256im) * x^0 + @test_broken subs(4x + 8y^2 - 6x^3, [x, real(y)] => [2 + 4im, 9]) + @test subs(4x + 8y^2 - 6x^3, [x, real(y)] => [2 + 4im, 9 + 0 * x^0]) == + 1184 + 112im + 144im * imag(y) - 8imag(y)^2 +end