How does sympy handle exponents to the 0.5 power?

Question

((gamma-(gamma**2-omega**2)**0.5)*(gamma+(gamma**2-omega**2)**0.5)).simplify()

The output is: gamma^2 - (gamma^2 -omega^2)^{1.0} $

However, I expected the result to be omega^2. I know in the sympy docs, it warns about being careful with floating point numbers, but I was under the impression that integers and also fractional powers of 2 (which can be represented exactly) were fine.

The following code correctly reproduces omega^2:

((gamma-(gamma**2-omega**2)**sym.Rational(1,2))*(gamma+(gamma**2-omega**2)**sym.Rational(1,2))).simplify()

Why does the first code not produce the expected result?

Answer 1

SymPy considers that there is a distinction between exact and inexact numbers. In this context floats like 0.5 and 1.0 are considered to be inexact and therefore it is not clear that x**1.0 is really equal to x or equal to something slightly different like say x**1.00000000000000000000001. That is because floats usually arise from floating point calculations which can have rounding errors. In your example the result is:

In [5]: from sympy import *

In [6]: gamma, omega = symbols('gamma, omega')

In [7]: e = ((gamma-(gamma**2-omega**2)**0.5)*(gamma+(gamma**2-omega**2)**0.5)).simplify()

In [8]: e
Out[8]: 
              1.0
 2   ⎛ 2    2⎞   
γ  - ⎝γ  - ω ⎠ 

If you want to tell SymPy that the 1.0 should be treated as an exact 1 then you can use SymPy's nsimplify function:

In [9]: nsimplify(e)
Out[9]: 
 2
ω