This is not necessarily a bug but it is a difficulty when using symengine with complex numbers, and I'm mostly asking about what a suitable course of action could be. I'll use the Python bindings in my examples for convenience, though the "problem" comes from the C++ engine.

If I have a Mul expression that involves a symbol multiplied by a complex float, for example:

import symengine
expr = 1j * symengine.Symbol("x")
# expr = (0.0 + 1.0*I) * x

If I now substitute x for some concrete value, I would ideally like to get bit-for-bit exact floating-point results, but I get an extra error term:

expr.subs({"x": 1})
# 6.12323399573677e-17 + 1.0*I

This error term comes in during SubsVisitor::bvisit(const Mul &) (inheriting from XReplaceVisitor) - the multiplication is evaluated involving a (0.0 + 1.0*I) ** 1.0 step in Mul::dict_add_term_new(), which eventually calls (essentially) std::pow(1i, 1.0) which is where the error term validly arises.

I would like to avoid introducing the error term here, since it's a symptom of how symengine performs the evaluation rather than one that's expected from how the expression is written. I'm guessing that general architecture wants the Mul evaluation tree have the coefficient written as a general exponential, but would it be possible/appropriate to insert a skip if the exponent is 1.0 in the appropriate power methods? It's not a great general solution, but for complex-valued expressions that don't naturally involve powers, it would be a big convenience.