SymEngine: dropping terms in a polynomial with near-zero coefficients

General Usage
1

I’ve been using SymEngine to get a symbolic polynomial representation of a function and have been impressed by its speed and amazed that it even works seemlessly with BigFloat. However, at least partially because the function I care about is based on floating point parameters parsed from a file, which were serialized with limited precision, I end up with a lot of terms with coefficients like

2.632918509595140993295845084408890582769689363958590871803066577611530482365e-116

I’d like to drop those terms. Is there a simple way to do that?

2

This might help:

EDIT: No, zchop on an expression returned by SymEngine throws a StackOverflow. This is a bug in ZChop.

3

But, this doesn’t look good to me:

julia> using SymEngine;

julia> a, b = symbols("a b");

julia> expr = a * b
a*b

julia> typeof(expr)
Basic

julia> Base.isiterable(Basic)
true

julia> for x in expr
       global z = x
       println(x)
       end
a*b

julia> z  === expr
true

This is consistent with iteration over objects whose type is a subtype of Number. It follows from

julia> supertype(Basic)
Number

Even so, it’s counterintutive.

4

I fixed this bug on the master branch of ZChop. zchop(x::Number) now returns x if x is neither Real nor Complex.

Unfortunately this does not provide a solution to the question. And I was unable to find documents or examples related to accessing parts of expressions in SymEngine.

5

Here’s a prototype that you can improve for your needs. (I’d like to get rid of if/elseif/else after implementing Add a method to create a object from class id and args · Issue #1478 · symengine/symengine · GitHub)

function chop_expr(x::Basic)
  cls = SymEngine.get_symengine_class(x)
  args = [chop_expr(arg) for arg in get_args(x)]
  if cls == :Add
     return sum(args)
  elseif cls == :Mul
     return prod(args)
  elseif cls == :Pow
     return args[1]^args[2]
  elseif cls in [:RealDouble, :RealMPFR, :ComplexDouble, :ComplexMPFR]
     y = N(x)
     # Do chopping here
     if abs(y) <= 1e-100
        y = 0
     end
     return Basic(y)
  else
     return x
  end
end

Btw, if all you are using are polynomials, you’ll find that SymEngine is slow for polynomials and a specialized polynomial library will be much faster.

1 Like
6

Great, thank you very much. @jlapeyre, thanks as well.