Spurious exponents in sympy expressions

Question

I build sp expressions (import sympy as sp) by multiplying / dividing others. Sometimes I find results like (expr1)

1.0*K1**1.11e-16*K2**1.11e-16/K3**1.11e-16*K4**0.8

while the correct expression is (expr2)

1.0*K4**0.8

My first question is what are possible ways of converting expr1 into expr2?
I can set an arbitrary threshold for the exponents of, say, 1e-10 to consider them spurious. I can concoct a recursive checking of the expression, but perhaps there is something more concise/pythonic. For the time being I can assume all powers are only at the first mul/div level, there are no additions/subtractions involved.

My second question is what are possible reasons for these spurious powers to appear?
There are no sp.N involved, e.g. The way I build my expressions is quite intricate to post here, so I welcome possible causes for this. I guess this is well within reach of knowledgeable fellows. In the meantime, I will try conceiving reasonable and representative MCVEs.

Related (some, marginally):

This happens because you are using floats rather than exact rational numbers. If you are not explicitly using float expressions like `0.2` then look instead for division like `1/5` and rewrite as `sympy.S(1)/5` or `sympy.Rational(1, 5)`. If you start with exact numbers then sympy will preserve that. The floats in the output can only be there because you gave sympy floats as input. — Oscar Benjamin, Jan 21 '23 at 23:49
You could try `nsimplify()` on your (input) expressions, as in [How do I replace floats with rationals in a sympy expression?](https://stackoverflow.com/questions/21005132/how-do-i-replace-floats-with-rationals-in-a-sympy-expression). See also [the sympy gotchas](https://docs.sympy.org/latest/tutorials/intro-tutorial/gotchas.html#two-final-notes-and) — JohanC, Jan 22 '23 at 00:14
@OscarBenjamin - Replacing with `sympy.Rational(...` in all expressions was tedious, but it worked. You could post it as an answer. — sancho.s ReinstateMonicaCellio, Jan 24 '23 at 02:39
@JohanC - I didn't try your options since `sympy.Rational` worked right away. — sancho.s ReinstateMonicaCellio, Jan 24 '23 at 02:40
`nsimplify()` works on a complete expression, while `Rational` goes one value at a time. E.g.`nsimplify(2.51 *x**0.5+x*0.33333333)` gives `251*sqrt(x)/100 + x/3`. `Rational(0.333333333)`gives `3002399748577931/9007199254740992`. Sometimes you might prefer one, sometimes the other. — JohanC, Jan 24 '23 at 06:27

score 1 · Accepted Answer · answered Jan 24 '23 at 15:38

As noted in the comments this happens because you are using floats rather than exact rational numbers:

In [16]: from sympy import *

In [17]: x = symbols('x')

In [18]: x**0.11/x**0.1/x**0.01
Out[18]: 
 -5.20417042793042e-18
x

The problem here is that in (binary) floating point 0.11 for example is not exactly equal to 0.1 + 0.01 because in fact none of these numbers floats is equal to the intended number that cannot be represented exactly.

The best approach in general is that in cases where you did not mean for a number to be approximate you should use exact rational numbers. There are different ways of creating exact numbers rational numbers:

In [19]: x**Rational('0.11')/x**Rational('0.1')/x**Rational('0.01')
Out[19]: 1

In [20]: x**Rational(11,100)/x**Rational(10,100)/x**Rational(1,100)
Out[20]: 1

In [22]: x**(S(11)/100)/x**(S(10)/100)/x**(S(1)/100)
Out[22]: 1

Another approach also suggested in the comments is to use nsimplify:

In [23]: x**0.11
Out[23]: 
 0.11
x    

In [24]: nsimplify(x**0.11)
Out[24]: 
  11
 ───
 100
x   

In [25]: x**Rational(0.11)
Out[25]: 
  7926335344172073
 ─────────────────
 72057594037927936
x

Observe that the exact value of the float 0.11 is not in fact equal to the mathematical number 0.11. What nsimplify does is to try to guess what number you really intended so it in fact performs an inexact conversion from float to Rational. This is useful as a convenience but cannot be expected to be always reliable so it is better to use rational numbers in the first place and have exact calculations throughout.

Another reason you might have these floats is because of calling evalf:

In [35]: e = sqrt(x)*pi

In [36]: e
Out[36]: π⋅√x

In [37]: e.evalf()
Out[37]: 
                  0.5
3.14159265358979⋅x   

In [38]: nfloat(e)
Out[38]: 3.14159265358979⋅√x

Here the nfloat function can be used to avoid calling evalf on the exponent. This is useful because having floats in exponents of symbols is particularly problematic.

Another approach is using functions for making roots e.g.:

In [39]: x**0.5
Out[39]: 
 0.5
x   

In [40]: x**(1/3)
Out[40]: 
 0.333333333333333
x                 

In [41]: sqrt(x)
Out[41]: √x

In [42]: cbrt(x)
Out[42]: 
3 ___In [39]: x**0.5
Out[39]: 
 0.5
x   

In [40]: x**(1/3)
Out[40]: 
 0.333333333333333
x                 

In [41]: sqrt(x)
Out[41]: √x

In [42]: cbrt(x)
Out[42]: 
3 ___
╲╱ x 

In [43]: root(x, 4)
Out[43]: 
4 ___
╲╱ x 
╲╱ x 

In [43]: root(x, 4)
Out[43]: 
4 ___
╲╱ x

Finally if you have an expression like the one that you showed it is not possible to use nsimplify at that stage because it will not guess that small numbers are supposed to be zero. Instead you can manually replace small floats though:

In [49]: K1, K2, K3, K4 = symbols('K1:5')

In [50]: e = 1.0*K1**1.11e-16*K2**1.11e-16/K3**1.11e-16*K4**0.8

In [51]: e
Out[51]: 
      1.11e-16   1.11e-16   -1.11e-16   0.8
1.0⋅K₁        ⋅K₂        ⋅K₃         ⋅K₄   

In [52]: nsimplify(e)
Out[52]: 
          111                   111              
  ───────────────────   ───────────────────      
  1000000000000000000   1000000000000000000   4/5
K₁                   ⋅K₂                   ⋅K₄   
─────────────────────────────────────────────────
                        111                      
                ───────────────────              
                1000000000000000000              
              K₃                                 

In [53]: e.replace(lambda t: isinstance(t, Float), lambda f: f if abs(f) > 1e-10 else 0)
Out[53]: 
      0.8
1.0⋅K₄

Spurious exponents in sympy expressions

1 Answers1