Consider changing backend for PolyForm arithmetic type promotion

[bowenszhu]

Short description

The current implementation of PolyForm division relies on DynamicPolynomials.Polynomial, which in turn utilizes MutableArithmetics.jl for type promotion. MutableArithmetics.jl, designed for numerical software such as JuMP, promotes division of two integers to a Float64. This behavior might not be ideal for symbolic computations where preserving integer types is often desired.

Detailed description

If the input expression contains a / operation, simplify and simplify_fractions call simplify_div.

SymbolicUtils.jl/src/simplify.jl

Lines 16 to 43 in 6713fa0

	function simplify(x;
	expand=false,
	polynorm=nothing,
	threaded=false,
	simplify_fractions=true,
	thread_subtree_cutoff=100,
	rewriter=nothing)
	if polynorm !== nothing
	Base.depwarn("simplify(..; polynorm=$polynorm) is deprecated, use simplify(..; expand=$polynorm) instead",
	:simplify)
	end
	f = if rewriter === nothing
	if threaded
	threaded_simplifier(thread_subtree_cutoff)
	elseif expand
	serial_expand_simplifier
	else
	serial_simplifier
	end
	else
	Fixpoint(rewriter)
	end
	x = PassThrough(f)(x)
	simplify_fractions && has_operation(x, /) ?
	SymbolicUtils.simplify_fractions(x) : x

SymbolicUtils.jl/src/polyform.jl

Lines 327 to 337 in 6713fa0

	function simplify_fractions(x; polyform=false)
	x = Postwalk(quick_cancel)(x)
	!needs_div_rules(x) && return x
	sdiv(a) = isdiv(a) ? simplify_div(a) : a
	expr = Postwalk(sdiv ∘ quick_cancel,
	similarterm=frac_similarterm)(Postwalk(add_with_div,
	similarterm=frac_similarterm)(x))

simplify_div transforms both numerator and denominator to PolyForms, and removes their greatest common divisor.

SymbolicUtils.jl/src/polyform.jl

Lines 277 to 280 in 6713fa0

	function simplify_div(d)
	d.simplified && return d
	ns, ds = polyform_factors(d, get_pvar2sym(), get_sym2term())
	ns, ds = rm_gcds(ns, ds)

rm_gcds involves the division of PolyForms.

SymbolicUtils.jl/src/polyform.jl

Lines 534 to 543 in 6713fa0

	function rm_gcds(ns, ds)
	ns = flatten_pows(ns)
	ds = flatten_pows(ds)
	for i = 1:length(ns)
	for j = 1:length(ds)
	g = _gcd(ns[i], ds[j])
	if !_isone(g)
	ns[i] = div(ns[i], g)
	ds[j] = div(ds[j], g)

SymbolicUtils.jl/src/polyform.jl

Line 86 in 6713fa0

Base.div(x::PolyForm, y::PolyForm) = $PF(div(x.p, y.p), mix_dicts(x, y)...)

As shown in the above code, our current implementation of PolyForm division relies on DynamicPolynomials.Polynomial division via div(x.p, y.p).

JuliaAlgebra/MultivariatePolynomials.jl uses jump-dev/MutableArithmetics.jl for arithmetic type promotion.
https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/blob/9eb6bf4973e27fdbf6f53d0d2043b476ad2a38c5/src/division.jl#L148-L150
https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/blob/9eb6bf4973e27fdbf6f53d0d2043b476ad2a38c5/src/division.jl#L388-L390
https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/blob/9eb6bf4973e27fdbf6f53d0d2043b476ad2a38c5/src/division.jl#L379-L384
https://github.com/jump-dev/MutableArithmetics.jl/blob/7344c9468072bd663feec7e4cdc1a0e98e34ea19/src/interface.jl#L112-L114
https://github.com/jump-dev/MutableArithmetics.jl/blob/7344c9468072bd663feec7e4cdc1a0e98e34ea19/src/interface.jl#L38-L44

In the last piece of code above, for inferring the resulting type of two Ints with the / operation, it obtains Float64 instead of Int or Rational by essentially executing typeof(Int(0) / Int(1)).
This is reasonable for numerical software, but we might want it to behave differently in our symbolic-focused software.

[hersle]

The simple change in JuliaAlgebra/MultivariatePolynomials.jl#294 would fix all examples in JuliaSymbolics/Symbolics.jl#1083.
My knowledge of that library is limited, but I think it might break other things there, and not be ready to merge (?)

[blegat]

I didn't want to force div of polynomials with integer coefficients to be floating point, it was just a case I didn't give much thought on. Indeed, div and rem for polynomials only make sense when the coefficient type is a field. For instance divrem(5x, 3x) should be 5/3 with zero remainder, not 1 with remainder 2x because the leading monomial of the remainder cannot be divisible by the leading monomial of the divisor. There is also pseudo_rem in case you have non-field coefficient like Int. For that reason, I assumed the user that needs to divide polynomials with integer coefficients would promote it to either rational or float in case there was a strong opinion on what the result should be. So there is a quick fix on your end:

julia> function to_field(p)
           T = coefficient_type(p)
           if T <: Integer
               return map_coefficients(c -> Rational{T}(c), p)
           else
               return p
           end
       end
to_field (generic function with 1 method)

julia> to_field(2x + 1)
1//1 + 2//1x

Now, defaulting to Rational could make sense, I'm having an attempt in JuliaAlgebra/MultivariatePolynomials.jl#296