Weird result in power operator

New to Julia
1

hi
there is a non stable result in power operator. this is the example. when I call the power with the variable the power not works. but when calling directly it
julia> a = -1.9345e6
-1.9345e6

julia> a^1.5
ERROR: DomainError with -1.9345e6:
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(::Float64) at .\math.jl:35
[2] ^(::Float64, ::Float64) at .\math.jl:782
[3] top-level scope at none:0

calling directly:
julia> -1.9345e6^1.5
-2.6906245536352706e9

matlab result:
-1.9345e+06^1.5

ans =

-2.6906e+09

calling with complex

julia> complex(a)^1.5
0.0 - 2.6906245536352706e9im

2

The - operator has lower precedence than the ^ operator, so the expression -1.9345e6^1.5 means -(1.9345e6^1.5) rather than (-1.9345e6)^1.5. If you try the a = -1.9345e6; a^1.5 version in Matlab you’ll get the same imaginary value as complex(a)^1.5 gives in Julia.

Also, please quote code in posts:

5 Likes
3

While preparing an introductory talk about Julia on my University, I am looking and Chris Rackauckas’s notes why Julia. I am puzzled with his example with powers, as 5^5 should return an Int (and it does), whereas 5^-5 should fail. Nevertheless 5^-5 executes correctly in repl, but when I define function

julia> f(x, y) = x^y - 5
julia> f(5, -5)
ERROR: DomainError with -5:

it fails as it should. I noticed the thread make x^-n equivalent to inv(x)^n for literal n by stevengj · Pull Request #24240 · JuliaLang/julia · GitHub, which says that 5^-5 should be lowered to inv(5)^5, which more or less explains the behavior in repl, but does not explain the discrepancy. Moreover, when I ask, which version is executed

julia> @which 5^-5
^(x::T, p::T) where T<:Integer in Base at intfuncs.jl:220

it points me to an inline function to ^(x::T, p::T) where {T<:Integer} = power_by_squaring(x,p), which makes perfect sense and indeed,

julia> Base.power_by_squaring(5, -5)
ERROR: DomainError with -5:

works as expected.

To conclude, there is something that I do not understand, and I would very much like to understand.

4

Literals are treated differently:

julia> Base.Meta.@lower 5^-5
:($(Expr(:thunk, CodeInfo(
1 ─ %1 = (Core.apply_type)(Base.Val, -5)
│   %2 = (%1)()
│   %3 = (Base.literal_pow)(^, 5, %2)
└──      return %3
))))

which is equivalent to

julia> Base.literal_pow(^, 5, Val(-5))
0.0003200000000000001

julia> f(x, y) = Base.literal_pow(^,x,y)
f (generic function with 1 method)

julia> f(5, Val(-5))
0.0003200000000000001

I don’t know what the reason behind this design is. Presumably to avoid having to type literal floats.

2 Likes
5

Because

julia> Meta.lower(Main, :(5^-5))
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope'
1 ─ %1 = Core.apply_type(Base.Val, -5)
│   %2 = (%1)()
│   %3 = Base.literal_pow(^, 5, %2)
└──      return %3
))))

note the literal_pow (which indeed comes from that PR), while in the function it is not constant, so is not parsed/lowered like that, but as an application of ^:

julia> @code_warntype f(5,-5)
Variables
  #self#::Core.Compiler.Const(f, false)
  x::Int64
  y::Int64

Body::Int64
1 ─ %1 = (x ^ y)::Int64
│   %2 = (%1 - 5)::Int64
└──      return %2

See the tools above to help you in general.

Is there something specific you want to achieve, eg have your code work in general, or is this just about understanding what is going on?

6

I have figure that it is hitting literal_pow, but I was quite confused that @which has not pointed me to the correct method that was used.

7

Thanks @mauro3 and @Tamas_Papp. I do not want to achieve anything special, I am just interested in internals, how does it work. Your explanation is great, yet it adds another thing I do not understand well.

8

Precisely, see #3024 (from the dawn of time).

Personally, I would prefer that users type inv(a)^n if that’s what they want, but that adds yet another item to the list of unexpected things that surprises people who started using Julia without bothering to read the manual.

In exchange for resolving that, the confusion is amplified and backshifted in time to the point where users become interested in finer points of the language, and encounter the magic that was used in the implementation.

1 Like
9

I confess that the last time I have read the manual was for Julia v0.3. I think I have a homework.

4 Likes
10

That remark was not targeted at you (but of course, re-reading the manual from time to time is great, a lot of effort goes into it and it and it is really useful), it was more about the reasons that I assume lead to the introduction of this feature.

Namely, the pivotal role of type stability in Julia is very easy to miss if one just skims through a short tutorial or similar. Yet it explains the design choices that answer questions like “why doesn’t √-2 work?”.

1 Like
11

I think your remark was very sound. Thanks a lot for opening me eyes.

12

See also `@which` and `@code_warntype` don't show `^Val{p}` lowering · Issue #21014 · JuliaLang/julia · GitHub for the @which bug.

13

I am unsure whether this should be a new thread, the title of this thread fits well…

Typing 16^16 in the REPL gives a result of 0 (Julia v1.4).

I guess this is related to the size of Int64 (and 2^64 logically gives the same answer, so does 16^17). But the fact that no error is raised here may be a problem, as one may take long to find out where the error comes from (in a larger calculation).

14

Integer overflow is a different issue which comes up regularly. The most recent discussion about it is here, including reasons for this behaviour, ways around it and some suggestions to make it easier on newcomers.

4 Likes