JuliaPy / SymPy.jl

Julia interface to SymPy via PyCall
http://juliapy.github.io/SymPy.jl/
MIT License
266 stars 61 forks source link

Incorrect `lambdify`? (exponentiation of negative `real`) #510

Open LebedevRI opened 1 year ago

LebedevRI commented 1 year ago

simplify() moves exponentiation into the modulo operation, but that apparently changes semantics as per externally-observable side-effects:

using SymPy
@syms x::real
@show f_sym = 0.17*x/(abs(x))^(0.6)
@show f_sym(-1)
f = lambdify(f_sym)
@show f(-1)
@show f_sym_simplified = simplify(f_sym)
@show f_sym_simplified(-1)
f_simplified = lambdify(f_sym_simplified)
@show f_simplified(-1)
f_sym = (0.17x) / abs(x) ^ 0.6 = 0.17*x/Abs(x)^0.6
f_sym(-1) = -0.170000000000000
f(-1) = -0.17
f_sym_simplified = simplify(f_sym) = 0.17*x/Abs(x^0.6)
f_sym_simplified(-1) = -0.170000000000000
DomainError with -1.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

Stacktrace:
 [1] throw_exp_domainerror(x::Float64)
   @ Base.Math ./math.jl:37
 [2] ^(x::Float64, y::Float64)
   @ Base.Math ./math.jl:1123
 [3] ^
   @ ./promotion.jl:444 [inlined]
 [4] __POW__
   @ ~/.julia/packages/SymPy/mpN0u/src/lambdify.jl:45 [inlined]
 [5] var"##294"(x::Int64)
   @ SymPy ./none:0
 [6] #invokelatest#2
   @ ./essentials.jl:816 [inlined]
 [7] invokelatest
   @ ./essentials.jl:813 [inlined]
 [8] (::SymPy.var"#120#121"{SymPy.var"###294"})(args::Int64)
   @ SymPy ~/.julia/packages/SymPy/mpN0u/src/lambdify.jl:248
 [9] top-level scope
   @ show.jl:1128

This is Julia Version 1.9.1, SymPy v1.1.9, on debian amd64 sid.

jverzani commented 1 year ago

Thanks for pointing this out. I have to think how best to address it. The issue that errors is the expression when lambdified includes an x^(0.6) which for a negative, real x will error. You can avoid that with x=-1.0 + 0im.

But that is unsatisfying. It seems SymPy shouldn't simplify this the way it does.

julia> @syms x::real;

julia> simplify(abs(x)^(0.6))
   0.6
│x│   

julia> simplify(1/abs(x)^(0.6))
│ -0.6│
│x    │