projekter/CPPolynomials.jl Secret

## CPPolynomials.jl
using MultivariatePolynomials
import DynamicPolynomials

# Here, we extend MultivariatePolynomials/DynamicPolynomials so that elementary complex-numbers related operations also work on polynomials.
# For this, we introduce special variables (though their type does not change, only their behavior under printing, substitution and some complex operations).
# Note that this requires a patch in DynamicPolynomials: a second internal constructor for PolyVar must be provided that allows to specify a custom id.
# We also define a new macro for DynamicPolynomials, @polycvar, which allows to directly define complex-valued variables.
# (The same could be done for the noncommuting case, which is left out here.)

# Prototyping: generic functions on all variables. All generic variables are real by default.
@eval MultivariatePolynomials begin
    iscomplex(x::AbstractVariable) = false;
    isrealpart(x::AbstractVariable) = false;
    isimagpart(x::AbstractVariable) = false;
    isconj(x::AbstractVariable) = false;
    ordvar(x::AbstractVariable) = x;
    Base.conj(x::AbstractVariable) = x;
    Base.real(x::AbstractVariable) = x;
    Base.imag(x::AbstractVariable) = 0;

    # 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)
    Base.conj(x::M) where {M<:AbstractMonomial} = convert(M, reduce(*, conj(var)^exp for (var, exp) in zip(variables(x), exponents(x)); init=constantmonomial(x)));
    Base.conj(x::V) where {V<:AbstractVector{<:AbstractMonomial}} = monovec(conj.(x));
    Base.conj(x::T) where {T<:AbstractTerm} = convert(T, conj(coefficient(x)) * conj(monomial(x)));
    Base.conj(x::P) where {P<:AbstractPolynomial} = convert(P, sum(conj(t) for t in x));
    LinearAlgebra.adjoint(x::AbstractVariable) = conj(x);
    LinearAlgebra.adjoint(x::AbstractMonomial) = conj(x);
    LinearAlgebra.adjoint(x::AbstractTerm) = _term(LinearAlgebra.adjoint(coefficient(x)), LinearAlgebra.adjoint(monomial(x)));

    realCoefficientTypePolynomialLike(::AbstractPolynomialLike{R}) where {R<:Real} = R;
    realCoefficientTypePolynomialLike(::AbstractPolynomialLike{Complex{R}}) where {R<:Real} = R;
    # 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::Union{AbstractMonomial,AbstractTerm})
                    # We replace every complex variable by its decomposition into real and imaginary part
                    substs = [var => real(var) + 1im*imag(var) for var in variables(x) if iscomplex(var)];
                    # subs throws an error if it doesn't receive at least one substitution
                    fullVersion = length(substs) > 0 ? subs(x, substs...) : polynomial(x);
                    # Now everything that is imaginary can be recognized by its coefficient
                    return convert(polynomialtype(fullVersion, realCoefficientTypePolynomialLike(fullVersion)),
                        mapcoefficients!($fun, fullVersion));
                end;
                Base.$fun(x::AbstractVector{<:AbstractMonomial}) = map(Base.$fun, x);
                Base.$fun(x::AbstractPolynomial{T}) where {T} = sum($fun(t) for t in x);
            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.
    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::AbstractMonomial, v::AbstractVariable)
        deg = 0;
        degC = 0;
        cV = conj(v);
        for (var, exp) in powers(m)
            @assert(!isrealpart(var) && !isimagpart(var));
            if var == v
                deg += exp;
            elseif var == cV
                degC += exp;
            end;
        end;
        return max(deg, degC);
    end;

    mindegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
        isempty(X) ? 0 : minimum(t -> degree_complex(t, args...), X);
    mindegree_complex(p::AbstractPolynomialLike, args...) = mindegree_complex(terms(p), args...);

    maxdegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
        mapreduce(t -> degree_complex(t, args...), max, X, init=0);
    maxdegree_complex(p::AbstractPolynomialLike, args...) = maxdegree_complex(terms(p), args...);

    extdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, args...) =
        (mindegree_complex(p, args...), maxdegree_complex(p, args...));

    function ordvar(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 = ordvar(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;

    function _show(io::IO, mime::MIME, var::AbstractVariable)
        base, indices = name_base_indices(var);
        if isconj(var)
            for c in base
                print(io, c, '\u0305');
            end;
        else
            print(io, base);
        end;
        if !isempty(indices)
            print_subscript(io, mime, indices);
        end;
        isrealpart(var) && print(io, "ᵣ");
        isimagpart(var) && print(io, "ᵢ");
    end;

    function _show(io::IO, mime::MIME"text/print", var::AbstractVariable)
        if isconj(var)
            for c in name(var)
                print(io, c, '\u0305');
            end;
        else
            print(io, name(var));
        end;
        isrealpart(var) && print(io, "ᵣ");
        isimagpart(var) && print(io, "ᵢ");
    end;
end;

using MultivariatePolynomials: iscomplex, isrealpart, isimagpart, isconj, ordvar, degree_complex,
    mindegree_complex, maxdegree_complex, extdegree_complex

# In the special case for DynamicPolynomials, we define complex-valued variables.
# We hijack the ordinary PolyVar for this and use the three most significant bits to encode
# which behavior our variable represents. Indeed, this will run us into trouble if we create
# more than 2^61 variables...
@eval DynamicPolynomials begin
    # Make sure the following inner constructor for PolyVar is inserted into the source:
    # PolyVar{C}(id::Int, name::AbstractString) where {C} = new(id, name);

    # We choose the following custom encoding
    # - look at the id of the PolyVar
    # - ordinary real variable: id ≥ 0
    # - complex variable: id < 0
    # - real part of complex variable: id & 0x4000000000000000 is set
    # - imaginary part of complex variable: id & 0x2000000000000000 is set
    # - conjugate variable: id < 0, but id & 0x6000000000000000 is unset
    const IsRealPart = Int(0x4000000000000000);
    const IsImagPart = Int(0x2000000000000000);
    const IsNotConj = IsRealPart | IsImagPart;

    using MultivariatePolynomials: iscomplex, isrealpart, isimagpart, isconj, ordvar, degree_complex,
        mindegree_complex, maxdegree_complex, extdegree_complex

    MP.iscomplex(x::PolyVar{C}) where {C} = x.id < 0;
    MP.isrealpart(x::PolyVar{C}) where {C} = x.id ≥ 0 && (x.id & IsRealPart) != 0;
    MP.isimagpart(x::PolyVar{C}) where {C} = x.id ≥ 0 && (x.id & IsImagPart) != 0;
    MP.isconj(x::PolyVar{C}) where {C} = x.id < 0 && (x.id & IsNotConj) == 0;
    function MP.ordvar(x::PolyVar{C}) where {C}
        if x.id ≥ 0
            if x.id & IsNotConj == 0
                return x; # this never was a complex variable
            else
                return PolyVar{C}(-(x.id & ~IsNotConj), x.name);
            end;
        else
            return PolyVar{C}(x.id | IsNotConj, x.name);
        end;
    end;

    function MP.similarvariable(p::PolyVar{C}, ::Type{Val{V}}) where {C,V}
        var = PolyVar{C}(string(V), Val(:complex));
        if isconj(p)
            return conj(var);
        elseif isrealpart(p)
            return real(var);
        elseif isimagpart(p)
            return imag(var);
        else
            return var;
        end;
    end;
    MP.similarvariable(p::PolyVar{C}, V::Symbol) where {C} = similarvariable(p, Val{V});

    function Base.conj(x::PolyVar{C}) where {C}
        !iscomplex(x) && return x;
        if isconj(x)
            return PolyVar{C}(x.id | IsNotConj, x.name);
        else
            return PolyVar{C}(x.id & ~IsNotConj, x.name);
        end;
    end;
    # for efficiency reasons
    function Base.conj(x::Monomial{true})
        cv = conj.(x.vars);
        perm = sortperm(cv, rev=true);
        return Monomial{true}(cv[perm], x.z[perm]);
    end;
    # we don't define complex-valued noncommuting monomials (yet)
    #Base.conj(x::Monomial{false}) = Monomial{false}(conj.(x.vars), x.z);

    function Base.real(x::PolyVar{C}) where {C}
        !iscomplex(x) && return x;
        # although not strictly necessary, we remove a possible conjugation
        return PolyVar{C}(-(x.id | IsNotConj) | IsRealPart, x.name);
    end;

    function Base.imag(x::PolyVar{C}) where {C}
        !iscomplex(x) && return x;
        # this is not type-stable, but we don't want to have another potential "variable", the imaginary part
        # of the conjugate
        if isconj(x)
            return Term{C,Int}(-1, PolyVar{C}(-(x.id | IsNotConj) | IsImagPart, x.name));
        else
            return PolyVar{C}(-(x.id | IsNotConj) | IsImagPart, x.name);
        end;
    end;

    # Finally, special handling of substitution. We redefine fillmap!, and in our case it is most
    # likely more efficient to use the noncommuting variant in every case.
    # Note that methods with type restriction is a bit problematic, as there are more generic methods that
    # we don't want to remove - so in order to be able to execute this cell multiple times, we choose a more
    # direct approach (which sadly relies on Julia internals...)
    map(Base.delete_method, m for m in methods(fillmap!)
                              if m.sig isa DataType &&
                                 (m.sig.parameters[2:end] == Core.svec(Any, Vector{DynamicPolynomials.PolyVar{false}}, DynamicPolynomials.MP.Substitution) ||
                                  m.sig.parameters[2:end] == Core.svec(Any, Vector{DynamicPolynomials.PolyVar{true}}, DynamicPolynomials.MP.Substitution)));
    function fillmap!(vals, vars::Vector{PolyVar{C}}, s::MP.Substitution) 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
        for j in eachindex(vars)
            # ordinary assignment
            if vars[j] == s.first
                vals[j] = s.second;
            # assignment to conjugate
            elseif conj(vars[j]) == s.first
                vals[j] = conj(s.second);
            # assignment to a complex variable, but only parts occur in the polynomial
            elseif MP.isrealpart(vars[j]) && vars[j] == real(s.first)
                vals[j] = real(s.second);
            elseif MP.isimagpart(vars[j])
                if vars[j] == imag(s.first)
                    vals[j] = imag(s.second);
                elseif vars[j] == imag(conj(s.first))
                    vals[j] = imag(conj(s.second));
                end;
            end;
            # assignment to parts of a complex variable, but the full appears in the polynomial
            # We don't support this for the moment (would require more complex logic, since two
            # substitutions may give the full value, or one part is left over - but can vals
            # hold such a leftover?).
        end;
    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 MP.$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;

    # define @polycvar macro. Note we have to be careful with the arrays, as complex variables have a negative id
    # and therefore their order must be reversed.
    function PolyVar{C}(name::AbstractString, ::Val{:complex}) where {C}
        id = parse(Int, string(gensym())[3:end]);
        return PolyVar{C}(-id, convert(String, name));
    end;

    function polyarraycvar(::Type{PV}, prefix, indices...) where {PV}
        return reverse(map(i -> PV("$(prefix)[$(join(i, ","))]", Val(:complex)), Iterators.reverse(Iterators.product(indices...))));
    end

    function buildpolycvar(::Type{PV}, var) where {PV}
        if isa(var, Symbol)
            return var, :($(esc(var)) = $PV($"$var", Val(:complex)));
        else
            isa(var, Expr) || error("Expected $var to be a variable name");
            Base.Meta.isexpr(var, :ref) || error("Expected $var to be of the form varname[idxset]");
            (2 ≤ length(var.args)) || error("Expected $var to have at least one index set");
            varname = var.args[1];
            prefix = string(varname);
            return varname, :($(esc(varname)) = polyarraycvar($PV, $prefix, $(esc.(var.args[2:end])...)));
        end
    end

    function buildpolycvars(::Type{PV}, args) where {PV}
        vars = Symbol[];
        exprs = [];
        for arg in args
            var, expr = buildpolycvar(PV, arg);
            push!(vars, var);
            push!(exprs, expr);
        end
        return vars, exprs;
    end;

    macro polycvar(args...)
        vars, exprs = buildpolycvars(PolyVar{true}, reverse(args));
        return :($(foldl((x, y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...)));
    end;
end;
	using MultivariatePolynomials
	import DynamicPolynomials

	# Here, we extend MultivariatePolynomials/DynamicPolynomials so that elementary complex-numbers related operations also work on polynomials.
	# For this, we introduce special variables (though their type does not change, only their behavior under printing, substitution and some complex operations).
	# Note that this requires a patch in DynamicPolynomials: a second internal constructor for PolyVar must be provided that allows to specify a custom id.
	# We also define a new macro for DynamicPolynomials, @polycvar, which allows to directly define complex-valued variables.
	# (The same could be done for the noncommuting case, which is left out here.)

	# Prototyping: generic functions on all variables. All generic variables are real by default.
	@eval MultivariatePolynomials begin
	iscomplex(x::AbstractVariable) = false;
	isrealpart(x::AbstractVariable) = false;
	isimagpart(x::AbstractVariable) = false;
	isconj(x::AbstractVariable) = false;
	ordvar(x::AbstractVariable) = x;
	Base.conj(x::AbstractVariable) = x;
	Base.real(x::AbstractVariable) = x;
	Base.imag(x::AbstractVariable) = 0;

	# 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)
	Base.conj(x::M) where {M<:AbstractMonomial} = convert(M, reduce(*, conj(var)^exp for (var, exp) in zip(variables(x), exponents(x)); init=constantmonomial(x)));
	Base.conj(x::V) where {V<:AbstractVector{<:AbstractMonomial}} = monovec(conj.(x));
	Base.conj(x::T) where {T<:AbstractTerm} = convert(T, conj(coefficient(x)) * conj(monomial(x)));
	Base.conj(x::P) where {P<:AbstractPolynomial} = convert(P, sum(conj(t) for t in x));
	LinearAlgebra.adjoint(x::AbstractVariable) = conj(x);
	LinearAlgebra.adjoint(x::AbstractMonomial) = conj(x);
	LinearAlgebra.adjoint(x::AbstractTerm) = _term(LinearAlgebra.adjoint(coefficient(x)), LinearAlgebra.adjoint(monomial(x)));

	realCoefficientTypePolynomialLike(::AbstractPolynomialLike{R}) where {R<:Real} = R;
	realCoefficientTypePolynomialLike(::AbstractPolynomialLike{Complex{R}}) where {R<:Real} = R;
	# 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::Union{AbstractMonomial,AbstractTerm})
	# We replace every complex variable by its decomposition into real and imaginary part
	substs = [var => real(var) + 1im*imag(var) for var in variables(x) if iscomplex(var)];
	# subs throws an error if it doesn't receive at least one substitution
	fullVersion = length(substs) > 0 ? subs(x, substs...) : polynomial(x);
	# Now everything that is imaginary can be recognized by its coefficient
	return convert(polynomialtype(fullVersion, realCoefficientTypePolynomialLike(fullVersion)),
	mapcoefficients!($fun, fullVersion));
	end;
	Base.$fun(x::AbstractVector{<:AbstractMonomial}) = map(Base.$fun, x);
	Base.$fun(x::AbstractPolynomial{T}) where {T} = sum($fun(t) for t in x);
	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.
	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::AbstractMonomial, v::AbstractVariable)
	deg = 0;
	degC = 0;
	cV = conj(v);
	for (var, exp) in powers(m)
	@assert(!isrealpart(var) && !isimagpart(var));
	if var == v
	deg += exp;
	elseif var == cV
	degC += exp;
	end;
	end;
	return max(deg, degC);
	end;

	mindegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
	isempty(X) ? 0 : minimum(t -> degree_complex(t, args...), X);
	mindegree_complex(p::AbstractPolynomialLike, args...) = mindegree_complex(terms(p), args...);

	maxdegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
	mapreduce(t -> degree_complex(t, args...), max, X, init=0);
	maxdegree_complex(p::AbstractPolynomialLike, args...) = maxdegree_complex(terms(p), args...);

	extdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, args...) =
	(mindegree_complex(p, args...), maxdegree_complex(p, args...));

	function ordvar(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 = ordvar(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;

	function _show(io::IO, mime::MIME, var::AbstractVariable)
	base, indices = name_base_indices(var);
	if isconj(var)
	for c in base
	print(io, c, '\u0305');
	end;
	else
	print(io, base);
	end;
	if !isempty(indices)
	print_subscript(io, mime, indices);
	end;
	isrealpart(var) && print(io, "ᵣ");
	isimagpart(var) && print(io, "ᵢ");
	end;

	function _show(io::IO, mime::MIME"text/print", var::AbstractVariable)
	if isconj(var)
	for c in name(var)
	print(io, c, '\u0305');
	end;
	else
	print(io, name(var));
	end;
	isrealpart(var) && print(io, "ᵣ");
	isimagpart(var) && print(io, "ᵢ");
	end;
	end;

	using MultivariatePolynomials: iscomplex, isrealpart, isimagpart, isconj, ordvar, degree_complex,
	mindegree_complex, maxdegree_complex, extdegree_complex

	# In the special case for DynamicPolynomials, we define complex-valued variables.
	# We hijack the ordinary PolyVar for this and use the three most significant bits to encode
	# which behavior our variable represents. Indeed, this will run us into trouble if we create
	# more than 2^61 variables...
	@eval DynamicPolynomials begin
	# Make sure the following inner constructor for PolyVar is inserted into the source:
	# PolyVar{C}(id::Int, name::AbstractString) where {C} = new(id, name);

	# We choose the following custom encoding
	# - look at the id of the PolyVar
	# - ordinary real variable: id ≥ 0
	# - complex variable: id < 0
	# - real part of complex variable: id & 0x4000000000000000 is set
	# - imaginary part of complex variable: id & 0x2000000000000000 is set
	# - conjugate variable: id < 0, but id & 0x6000000000000000 is unset
	const IsRealPart = Int(0x4000000000000000);
	const IsImagPart = Int(0x2000000000000000);
	const IsNotConj = IsRealPart \| IsImagPart;

	using MultivariatePolynomials: iscomplex, isrealpart, isimagpart, isconj, ordvar, degree_complex,
	mindegree_complex, maxdegree_complex, extdegree_complex

	MP.iscomplex(x::PolyVar{C}) where {C} = x.id < 0;
	MP.isrealpart(x::PolyVar{C}) where {C} = x.id ≥ 0 && (x.id & IsRealPart) != 0;
	MP.isimagpart(x::PolyVar{C}) where {C} = x.id ≥ 0 && (x.id & IsImagPart) != 0;
	MP.isconj(x::PolyVar{C}) where {C} = x.id < 0 && (x.id & IsNotConj) == 0;
	function MP.ordvar(x::PolyVar{C}) where {C}
	if x.id ≥ 0
	if x.id & IsNotConj == 0
	return x; # this never was a complex variable
	else
	return PolyVar{C}(-(x.id & ~IsNotConj), x.name);
	end;
	else
	return PolyVar{C}(x.id \| IsNotConj, x.name);
	end;
	end;

	function MP.similarvariable(p::PolyVar{C}, ::Type{Val{V}}) where {C,V}
	var = PolyVar{C}(string(V), Val(:complex));
	if isconj(p)
	return conj(var);
	elseif isrealpart(p)
	return real(var);
	elseif isimagpart(p)
	return imag(var);
	else
	return var;
	end;
	end;
	MP.similarvariable(p::PolyVar{C}, V::Symbol) where {C} = similarvariable(p, Val{V});

	function Base.conj(x::PolyVar{C}) where {C}
	!iscomplex(x) && return x;
	if isconj(x)
	return PolyVar{C}(x.id \| IsNotConj, x.name);
	else
	return PolyVar{C}(x.id & ~IsNotConj, x.name);
	end;
	end;
	# for efficiency reasons
	function Base.conj(x::Monomial{true})
	cv = conj.(x.vars);
	perm = sortperm(cv, rev=true);
	return Monomial{true}(cv[perm], x.z[perm]);
	end;
	# we don't define complex-valued noncommuting monomials (yet)
	#Base.conj(x::Monomial{false}) = Monomial{false}(conj.(x.vars), x.z);

	function Base.real(x::PolyVar{C}) where {C}
	!iscomplex(x) && return x;
	# although not strictly necessary, we remove a possible conjugation
	return PolyVar{C}(-(x.id \| IsNotConj) \| IsRealPart, x.name);
	end;

	function Base.imag(x::PolyVar{C}) where {C}
	!iscomplex(x) && return x;
	# this is not type-stable, but we don't want to have another potential "variable", the imaginary part
	# of the conjugate
	if isconj(x)
	return Term{C,Int}(-1, PolyVar{C}(-(x.id \| IsNotConj) \| IsImagPart, x.name));
	else
	return PolyVar{C}(-(x.id \| IsNotConj) \| IsImagPart, x.name);
	end;
	end;

	# Finally, special handling of substitution. We redefine fillmap!, and in our case it is most
	# likely more efficient to use the noncommuting variant in every case.
	# Note that methods with type restriction is a bit problematic, as there are more generic methods that
	# we don't want to remove - so in order to be able to execute this cell multiple times, we choose a more
	# direct approach (which sadly relies on Julia internals...)
	map(Base.delete_method, m for m in methods(fillmap!)
	if m.sig isa DataType &&
	(m.sig.parameters[2:end] == Core.svec(Any, Vector{DynamicPolynomials.PolyVar{false}}, DynamicPolynomials.MP.Substitution) \|\|
	m.sig.parameters[2:end] == Core.svec(Any, Vector{DynamicPolynomials.PolyVar{true}}, DynamicPolynomials.MP.Substitution)));
	function fillmap!(vals, vars::Vector{PolyVar{C}}, s::MP.Substitution) 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
	for j in eachindex(vars)
	# ordinary assignment
	if vars[j] == s.first
	vals[j] = s.second;
	# assignment to conjugate
	elseif conj(vars[j]) == s.first
	vals[j] = conj(s.second);
	# assignment to a complex variable, but only parts occur in the polynomial
	elseif MP.isrealpart(vars[j]) && vars[j] == real(s.first)
	vals[j] = real(s.second);
	elseif MP.isimagpart(vars[j])
	if vars[j] == imag(s.first)
	vals[j] = imag(s.second);
	elseif vars[j] == imag(conj(s.first))
	vals[j] = imag(conj(s.second));
	end;
	end;
	# assignment to parts of a complex variable, but the full appears in the polynomial
	# We don't support this for the moment (would require more complex logic, since two
	# substitutions may give the full value, or one part is left over - but can vals
	# hold such a leftover?).
	end;
	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 MP.$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;

	# define @polycvar macro. Note we have to be careful with the arrays, as complex variables have a negative id
	# and therefore their order must be reversed.
	function PolyVar{C}(name::AbstractString, ::Val{:complex}) where {C}
	id = parse(Int, string(gensym())[3:end]);
	return PolyVar{C}(-id, convert(String, name));
	end;

	function polyarraycvar(::Type{PV}, prefix, indices...) where {PV}
	return reverse(map(i -> PV("$(prefix)[$(join(i, ","))]", Val(:complex)), Iterators.reverse(Iterators.product(indices...))));
	end

	function buildpolycvar(::Type{PV}, var) where {PV}
	if isa(var, Symbol)
	return var, :($(esc(var)) = $PV($"$var", Val(:complex)));
	else
	isa(var, Expr) \|\| error("Expected $var to be a variable name");
	Base.Meta.isexpr(var, :ref) \|\| error("Expected $var to be of the form varname[idxset]");
	(2 ≤ length(var.args)) \|\| error("Expected $var to have at least one index set");
	varname = var.args[1];
	prefix = string(varname);
	return varname, :($(esc(varname)) = polyarraycvar($PV, $prefix, $(esc.(var.args[2:end])...)));
	end
	end

	function buildpolycvars(::Type{PV}, args) where {PV}
	vars = Symbol[];
	exprs = [];
	for arg in args
	var, expr = buildpolycvar(PV, arg);
	push!(vars, var);
	push!(exprs, expr);
	end
	return vars, exprs;
	end;

	macro polycvar(args...)
	vars, exprs = buildpolycvars(PolyVar{true}, reverse(args));
	return :($(foldl((x, y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...)));
	end;
	end;