Issues with Rationals in Symbolics.jl

Lately I have been playing with the idea of writing a simple symbolic integrator based on the book by Bronstein “Symbolic Integration 1”, and have been more or less writing it from the ground up. So I have been writing things like a polynomial division algorithm.

This was all going well until I started getting rounding issues with the function that I wrote to divide polynomials.

At this point the solution that I decided to try was to convert the values of the coefficients of all of the equations into rationals. This seems to have worked for the issue that I was having however this brings up two new issues, both of witch seem to be part of symbolics.

The first issue, and so far this seems to be just a minor issue, is that if I do

using  Symbolics
@variables x
A = 1//1 * x
B = Symbolics.coeff(A, x)

I would expect that this would return a rational{Int64} for the typeof(B) like it will if you use
A = 5//1 *x
however it seems to insist on returning a Int64. This seems to happen in all situations so it is not just a consequence of how I did it here. So far this has not seemed to give me any issues however it’s not the result that I expected.

The real issue that I am having is that if I do

using Symbolics
@variables t
B = (((t^2)^3)*(((800//1)*(t^3) / ((800//1)*(t^3) 
    - (14//1)*t))^2)) / ((( (60//1)*(t^2))^2)
    * ( (14//1)*t)^3)

B = Symbolics.simplify(B)

The result is

ERROR: MethodError: no method matching div(::PolyForm{Real}, ::Rational{Int64}, ::RoundingMode{:ToZero})

Closest candidates are:
  div(::MultivariatePolynomials.AbstractPolynomialLike, ::Number, ::Any...)
   @ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/ckbfK/src/division.jl:8
  div(::P, ::Real, ::RoundingMode) where P<:Dates.Period
   @ Dates /usr/share/julia/stdlib/v1.9/Dates/src/periods.jl:81
  div(::Integer, ::Rational, ::RoundingMode)
   @ Base rational.jl:449

  [1] div(a::PolyForm{Real}, b::Rational{Int64})
    @ Base ./div.jl:47
  [2] rm_gcds(ns::Vector{Any}, ds::Vector{Any})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:538
  [3] simplify_div(d::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:275
  [4] (::SymbolicUtils.var"#sdiv#126")(a::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:328
  [5] call_composed
    @ ./operators.jl:1034 [inlined]
  [6] (::ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)})(x::SymbolicUtils.BasicSymbolic{Real}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./operators.jl:1031
  [7] (::ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)})(x::SymbolicUtils.BasicSymbolic{Real})
    @ Base ./operators.jl:1031
  [8] (::SymbolicUtils.Rewriters.Walk{:post, ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)}, typeof(SymbolicUtils.frac_similarterm), false})(x::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils.Rewriters ~/.julia/packages/SymbolicUtils/Oyu8Z/src/rewriters.jl:200
  [9] simplify_fractions(x::SymbolicUtils.BasicSymbolic{Real}; polyform::Bool)
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:330
 [10] simplify_fractions(x::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:322
 [11] simplify(x::SymbolicUtils.BasicSymbolic{Real}; expand::Bool, polynorm::Nothing, threaded::Bool, simplify_fractions::Bool, thread_subtree_cutoff::Int64, rewriter::Nothing)
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/simplify.jl:42
 [12] simplify(x::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/simplify.jl:16
 [13] simplify(n::Num; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:146
 [14] simplify(n::Num)
    @ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:146
 [15] top-level scope
    @ Untitled-1:17

The Julia VERSION that is running is v"1.9.3" and the Symbolics version is Symbolics v5.5.1
This is not the equation that originally gave me this error as I simplified it as much as I could without losing the actual error message. However I don’t really know what I need to do to stop this from coming up or if I need to go back and reconsider how I am doing things.

1 Like

Seems worth an issue

Looks like this was fixed? If I try the MWE on Julia 1.9.3 and Symbolics 5.10.0 it doesn’t give the MethodError anymore:

using Symbolics
@variables t
B = (((t^2)^3)*(((800//1)*(t^3) / ((800//1)*(t^3) 
    - (14//1)*t))^2)) / ((( (60//1)*(t^2))^2)
    * ( (14//1)*t)^3)

B = Symbolics.simplify(B)

with the result

(50(t^3)) / (3087((-7.0 + 400.0(t^2))^2))

What’s a bit weird is that using simplify turns some of the fractions into Float64 which is probably not what we want in the end (can lead to rounding errors again later).