pythonnumpysympysymengine

How to apply evalf() to an expression from symengine.py?


We couldn't find any evalf() in the package symengine for python. Does there exist any?

More detailed explanation of our problem.

We are trying to switch to symengine instead of sympy because we percept essential speedup. We use symbolic algebra tools to solve problem of optimisation which has applications in combinatorial random structure generation.

After some testing on sympy, we discovered that solving some "innocent-looking" optimisation problems requires intermediate computations in range 1e+24 and beyond (but in the end they finally provide the correct answer). Then I discovered that I misused sympy's function subs() because it is less precise than evalf() with substitution dictionary.

This issue is descried in the manual of the function subs:

If the substitution will be followed by numerical
evaluation, it is better to pass the substitution to
evalf as

>>> (1/x).evalf(subs={x: 3.0}, n=21)
0.333333333333333333333

rather than

>>> (1/x).subs({x: 3.0}).evalf(21)
0.333333333333333314830

as the former will ensure that the desired level of precision is obtained. In fact, substituting e+24 using standard sympy.subs() function, or symengine.subs() throws an infinity, though e+24 is still in the range of numpy.float64, which I suspect to be a standard type that my machine casts into, when I call built-in float(...) function.

We are aware of "lambdify", which is a recommended step if you do repetitive substitutions and want high numerical precision, but this will be the next step. We couldn't use the manuals from python interpreter because the code is precompiled, and on the internet it is also difficult to find any. Looking inside the source also didn't help (but maybe I missed something).


Solution

  • There's no evalf on symengine yet. Here's how the above can be achieved,

    In [14]: (1/x).subs({x: 3}).n(73, real=True)
    Out[14]: 0.333333333333333333333
    

    Note that if you do like below, it will do the substitution in precision 53 (15 decimal digits)

    In [15]: (1/x).subs({x: 3.0}).n(73, real=True)
    Out[15]: 0.333333333333333314830
    

    You can use subs with Float(3.0, dps=21) instead of just 3.0 to do the calculations in higher precision. Below example will work with both symengine and sympy.

    In [16]: (1/x).subs({x: Float(3.0, dps=21)})
    Out[16]: 0.333333333333333333333
    

    Note that .n() in symengine prefers precision in binary digits instead of decimal digits and real=True needs to be given for real domain, otherwise it'll assume complex and give you a complex number with imaginary part 0. and will be a little slower than real=True.