Complex-valued variables

Project.toml

@@ -13,7 +13,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-MultivariatePolynomials = "0.5"
+MultivariatePolynomials = "0.5.3"
 MutableArithmetics = "1"
 Reexport = "1"
 julia = "1"

src/DynamicPolynomials.jl

@@ -72,11 +72,11 @@ function MP.similar_variable(
 ) where {V,M,S}
     return MP.similar_variable(P, S)
 end
-function MP.similar_variable(
-    ::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
-    s::Symbol,
-) where {V,M}
-    return Variable(string(s), V, M)
+function MP.similar_variable(p::PolyType{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,
 end
 
 include("promote.jl")

src/comp.jl

@@ -27,14 +27,19 @@ function (==)(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return x.variable_order.order.id == y.variable_order.order.id
+    return x.variable_order.order.id == y.variable_order.order.id &&
+           x.kind == y.kind
 end
 
 function Base.isless(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return isless(y.variable_order.order.id, x.variable_order.order.id)
+    if x.variable_order.order.id == y.variable_order.order.id
+        return isless(y.kind, x.kind)
+    else
+        return isless(y.variable_order.order.id, x.variable_order.order.id)
+    end
 end
 
 # Comparison of Monomial

src/mono.jl

@@ -113,3 +113,10 @@ function __add_variables!(
     mono.z[map] = tmp
     return
 end
+
+# for efficiency reasons
+function Base.conj(x::Monomial{V,M}) where {V<:Commutative,M}
+    cv = conj.(x.vars)
+    perm = sortperm(cv, rev = true)
+    return Monomial{V,M}(cv[perm], x.z[perm])
+end

src/monomial_vector.jl

@@ -91,6 +91,28 @@
 function MultivariatePolynomials.maxdegree(x::MonomialVector)
     return isempty(x) ? 0 : maximum(sum.(x.Z))
 end
+# Complex-valued degrees for monomial vectors
+for (fun, call, def, ret) in [
+    (:extdegree_complex, :extrema, (0, 0), :((min(v1, v2), max(v1, v2)))),
+    (:mindegree_complex, :minimum, 0, :(min(v1, v2))),
+    (:maxdegree_complex, :maximum, 0, :(max(v1, v2))),
+]
+    eval(quote
+        function MultivariatePolynomials.$fun(x::MonomialVector)
+            isempty(x) && return $def
+            vars = variables(x)
+            @assert(!any(isrealpart, vars) && !any(isimagpart, vars))
+            grouping = isconj.(vars)
+            v1 = $call(sum(z[grouping]) for z in x.Z)
+            grouping = map(!, grouping)
+            v2 = $call(sum(z[grouping]) for z in x.Z)
+            return $ret
+        end
+    end)
+end
+# faster complex-related functions
+Base.isreal(x::MonomialVector) = all(isreal, x.vars)
+Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z)
 
 MP.variables(m::Union{Monomial,MonomialVector}) = m.vars
 
@@ -121,7 +143,11 @@
 # TODO replace by MP function
 function _error_for_negative_degree(deg)
     if deg < 0
-        throw(ArgumentError("The degree should be a nonnegative number but the provided degree `$deg` is negative."))
+        throw(
+            ArgumentError(
+                "The degree should be a nonnegative number but the provided degree `$deg` is negative.",
+            ),
+        )
     end
 end
 

src/subs.jl

@@ -9,23 +9,73 @@
 
 function fillmap!(
     vals,
-    vars::Vector{<:Variable{<:NonCommutative}},
+    vars::Vector{<:Variable{C}},
     s::MP.Substitution,
-)
-    for j in eachindex(vars)
-        if vars[j] == s.first
-            vals[j] = s.second
+) where {C}
+    # We may assign a complex or real variable to its value (ordinary substitution).
+    # We may also assign a complex value to its conjugate, or just the real or imaginary parts
+    # Any combination of z, conj(z), real(z), imag(z), imag(conj(z)) can occur in either the polynomial or the substitution,
+    # and we must handle all of them correctly.
+    # Note: This may or may not work... Issues can arise if the substitutions contain the real and imaginary (or only one of
+    # those) of a variable separately whenever vals is not of the correct type:
+    # - Unless subs() originally had a polynomial-valued rhs, vals will be scalars, monomials, or terms. So when we try to
+    #   assign a polynomial to its value (which is necessary, as the one-step substitution of only the real or only the
+    #   imaginary part is incomplete), this is an impossible conversion.
+    # - The coefficients in vals might not be complex-valued; but to assign only a part of the variable, we necessarily need to
+    #   introduce an explicit imaginary coefficient to the value.
+    # Currently, we don't do anything to catch these errors.
+    if s.first.kind == REAL
+        for j in eachindex(vars)
+            if vars[j] == s.first
+                vals[j] = s.second
+                C == Commutative && break
+            end
+        end
+    else
+        for j in eachindex(vars)
+            if vars[j].variable_order.order.id ==
+               s.first.variable_order.order.id
+                if s.first.kind == COMPLEX || s.first.kind == CONJ
+                    value = s.first.kind == CONJ ? conj(s.second) : s.second
+                    if vars[j].kind == COMPLEX
+                        vals[j] = value
+                    elseif vars[j].kind == CONJ
+                        vals[j] = conj(value)
+                    elseif vars[j].kind == REAL_PART
+                        vals[j] = real(value)
+                    else
+                        vals[j] = imag(value)
+                    end
+                elseif s.first.kind == REAL_PART
+                    isreal(s.second) || error(
+                        "Cannot assign a complex value to the real part of an expression",
+                    )
+                    value = real(s.second) # just to make sure the type is correct
+                    if vars[j].kind == COMPLEX
+                        vals[j] = value + im * imag(vals[j])
+                    elseif vars[j].kind == CONJ
+                        vals[j] = value - im * imag(vals[j])
+                    elseif vars[j].kind == REAL_PART
+                        vals[j] = value
+                    end
+                    # else we know the real part but use the imaginary part; do nothing
+                else
+                    @assert(s.first.kind == IMAG_PART)
+                    isreal(s.second) || error(
+                        "Cannot assign a complex value to the imaginary part of an expression",
+                    )
+                    value = real(s.second) # just to make sure the type is correct
+                    if vars[j].kind == COMPLEX
+                        vals[j] = real(vals[j]) + im * value
+                    elseif vars[j].kind == CONJ
+                        vals[j] = real(vals[j]) - im * value
+                    elseif vars[j].kind == IMAG_PART
+                        vals[j] = value
+                    end
+                    # else we know the imaginary part but use the real part; do nothing
+                end
+            end
         end
-    end
-end
-function fillmap!(
-    vals,
-    vars::Vector{<:Variable{<:Commutative}},
-    s::MP.Substitution,
-)
-    j = findfirst(isequal(s.first), vars)
-    if j !== nothing
-        vals[j] = s.second
     end
 end
 function fillmap!(vals, vars, s::MP.AbstractMultiSubstitution)