roundingsympyexponentiation

Spurious exponents in sympy expressions


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):

  1. Evaluating coefficients but not exponents in sympy
  2. How to isolate the exponent in Sympy (how to `match` large formulas in Sympy)
  3. Sympy - Find all expressions of a certain level

Solution

  • 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₄