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 @@ end
 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 @@ end
 # 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,63 @@ end
 
 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 follow the following rules:
+    # - If a single substitution rule determines the value of a real variable, just substitute it.
+    # - If a single substitution rule determines the value of a complex variable or its conjugate, substitute the appropriate
+    #   value whereever something related to this variable is found (i.e., the complex variable, the conjugate variable, or
+    #   its real or imaginary part)
+    # - If a single substitution rule determines the value of the real or imaginary part of a complex variable alone, then only
+    #   replace the real or imaginary parts if they occur explicitly. Don't do a partial substitution, i.e., `z` with the rule
+    #   `zᵣ => 1` is left alone and not changed into `1 + im*zᵢ`. Even if both the real and imaginary parts are substituted as
+    #   two individual rules (which we don't know of in this method), `z` will not be replaced.
+    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 == REAL_PART
+                        vals[j] = value
+                    end
+                    # don't do a partial substitution
+                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 == IMAG_PART
+                        vals[j] = value
+                    end
+                    # don't do a partial substitution
+                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)