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)
println(B)
typeof(B)

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
  ...

Stacktrace:
  [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).

Recently retried the OP’s example.

# test_simplify_rational.jl
using Symbolics

begin
    @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)
 
    @show B = Symbolics.simplify(B)
end

Now the output result is

B = Symbolics.simplify(B) =
(50(t^3)) / (3087((-(7//1) + (400//1)*(t^2))^2))

The coefficients are now returned as integer or rational numbers rather than being converted to floats.

Pleased to see this result as this was the sticking point that kept me from switching from Python/Sympy to Julia/Symbolics for symbolic computation and it now appears to have been addressed.

1 Like