Small example with different values when using PyPI symengine and shared library

[iyanmv]

I recently noticed a numerical difference when using the symengine wheel package from PyPI and the one provided by Arch Linux. This leads to a qiskit unit test to fail.

Here is a small example :

import symengine

f = symengine.Symbol("f")
expr = f / 1000

expr1 = expr * 1e9
expr2 = expr * 10**9

value1 = float(expr1.subs("f", 3.14))
value2 = float(expr2.subs("f", 3.14))

print(value1)
# 3139999.9999999995 with Arch Linux python-symengine
# 3140000.0 with PyPI
print(value2)
# 3140000.0 in both cases

I can imagine it is not a good practice to do expr * 1e9 in the first place, but is this documented somewhere?
And why if I manually create the expression I get the correct result? I.e.,

float((1000000.0*f).subs("f", 3.14))
# 3140000.0 in both cases

[isuruf]

Can you do print(expr1, expr2) in both cases?

[iyanmv]

Can you do print(expr1, expr2) in both cases?

# PyPI package
print(expr1, expr2)
# 1000000.0*f 1000000*f

# Arch packages
print(expr1, expr2)
# 1000000.0*f 1000000*f

There is no visible difference here, only later when substituting the symbol.

[iyanmv]

I think the difference is already visible here:

float(expr1.args[0])
# 999999.9999999998

[rikardn]

Interesting! The difference between expr1 and expr2 is that the constant factors are different on the python side. 1e9 has the type float and 10**9 has the type int. It makes sense that expr2 is exact since it is working with integers until the final substitution. expr1 is affected by a small rounding which is strange to me since the operation 1000000000.0 / 1000.0 * 3.14 should have no problems with rounding errors (I checked it for double precision floats).

[iyanmv]

Yes, exactly, but if there is any rounding error it should be in the 14 or 15 digit, not in the 10th. Were you able to replicate the issue on your side?

There are the flags used by Arch Linux to compile symengine:

-DCMAKE_INSTALL_PREFIX=/usr \
-DBUILD_SHARED_LIBS=ON \
-DWITH_TCMALLOC=ON \
-DWITH_PTHREAD=ON \
-DWITH_SYMENGINE_THREAD_SAFE=ON \
-DWITH_ARB=ON \
-DWITH_ECM=ON \
-DWITH_LLVM=ON \
-DWITH_MPFR=ON \
-DWITH_MPC=ON \
-DWITH_PRIMESIEVE=ON \
-DWITH_BOOST=ON \
-DWITH_COTIRE=OFF \
-DWITH_SYSTEM_CEREAL=ON

I can start to disable libraries to check if some of them are causing the issue. Any other ideas?

[rikardn]

No, I couldn't reproduce it, but I only tried with my installation.

Could you try the following using the "bad" symengine version?

from symengine import *

RealDouble(1000000000.0) / Integer(1000) * RealDouble(3.14)

[iyanmv]

That works as expected:

RealDouble(1000000000.0) / Integer(1000) * RealDouble(3.14)
# 3140000.0

[iyanmv]

Oh, I think the problem is with Rational:

from symengine import *

float(Rational(1, 1000))
# 0.0009999999999999998

So the f / 1000 is actually a Rational and the problem starts there.

expr = f / 1000
float(expr.subs("f", 1))
# 0.0009999999999999998

This works even if the factor is a float:

expr = f * 1e-3
float(expr.subs("f", 1))
# 0.001