Exponentiation by a negative power throws a DomainError

[bensadeghi]

julia> 2 ^ -2
ERROR: DomainError()
 in power_by_squaring at intfuncs.jl:84
 in ^ at intfuncs.jl:109

This fix works for me:

function power_by_squaring(x, p::Integer)
    x = to_power_type(x)
    if p == 1
        return x
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
##        throw(DomainError())
        return 1.0 / power_by_squaring(x, -p)
    end

[pao]

Did you intend to close this?

[JeffBezanson]

Either way, I don't think we're going to change this.

[bensadeghi]

Yes, I realized that it was a silly misunderstanding. Negative float exponents work fine:

2 ^ -2.0
3 .^ -[0:6.0]

[diegozea]

I don't understand why this error on Julia.

On Julia:

julia> 2^(-2)
ERROR: DomainError()
 in power_by_squaring at intfuncs.jl:84
 in ^ at intfuncs.jl:109

On R:

> 2^(-2)
[1] 0.25

On Python:

>>> 2**(-2)
0.25

People can't do this:

julia> 3 * 10 ^ (-6)
ERROR: DomainError()
 in power_by_squaring at intfuncs.jl:84
 in ^ at intfuncs.jl:109

Maybe a better error message?

[StefanKarpinski]

Type-stability, which other dynamic language completely ignore. There's a long discussion in an old issue about this. This was the best compromise we could come up with, although it is admittedly annoying.

[bsxfan]

So the integers are not closed w.r.t. division and inversion (raising to negative integer powers). Julia addresses this by:

• automatically converting inv(2) to float.
• automatically converting 4/2 and 2/4 etc. to float.
• but throwing DomainError for 2^(-1). Can't this just be converted to float automatically like the other cases? But then again, I see that this solution would break the rule that argument values should not determine output type.

P.S. In contrast, rationals CAN of course be raised to negative integer powers:

(2//3)^(-1) == 3//2

But be careful: 2//3^(-1) does throw DomainError, because ^ binds tighter than //. If 3^(-1) would instead automatically produce a float, we would get the same behaviour as for 2//3^(-1.0), which---quite reasonably---throws: ERROR: no method //(Int32,Float64).

[bsxfan]

Hmm, would it perhaps be possible to define a new operator: "^-", defined so that:

^-(x::Integer,n::Integer) = ^(float(x)^(-n)) 

[StefanKarpinski]

Oh, that's dastardly. I think I might like it. Because most of the time this is a purely syntactic negation that is actually quite safe. The case of a^b where a and b are both integers and b only sometimes ends up negative is unusual and often would actually constitute a bug.

[bsxfan]

This would make the relationship between ^ and ^- similar to the relationship between * and / for integer arguments:

• The * and ^always give integer results, while / and ^- always give floats.
• For both / and ^- there are special argument values for which the result could be accurately represented by integers, e.g. 4/2 and 2^-(-1), but which are returned as floats regardless.

[panlanfeng]

Perhaps you guys can add some more message for this error. It is really surprising to see it the first time.

[jiahao]

We should improve the error message as part of #4744.

[miguelbarao]

Negative powers are just too common.
It would be really helpful and EXPECTED if 2^-1 worked the same as 1/2.
There could be another operator like 2^^3 for type stability, if that's needed.

Some time ago I used a lot from __future__ import division in python and wouldn't like to follow the same path again here. I am really sorry this issue is closed.

(just a user)

[nalimilan]

@miguelbarao Replacing 2^-1 with 2.0^-1 is quite easy when porting code, and type stability is essential there if you want your code to be fast. It's more annoying when working at the REPL.

[timholy]

@miguelbarao, the problem is that currently "fast = type stable." You can have a slow ^, but then no one will use it, and everyone will face the same issue with ^^.

[ivarne]

I think the problem here is really that 2 and -1 parse as Int, instead of Float64. If we treat all numbers as Float64, this problem would go away, but I don't think losing convenient syntax for Integers is worth the convenience. Currently you get a hard DomainError, and adding the . is a pretty easy fix.

[StefanKarpinski]

If we treat all numbers as Float64, this problem would go away

What is this, JavaScript? :-P

[StefanKarpinski]

I'm still not entirely convinced that having a ^- operator is a terrible idea, but it's certainly is tricksy.

[pao]

I think ^- opens you up to "why doesn't a = -1; 2^a work when 2^-1 does?" sorts of questions.

[StefanKarpinski]

That it certainly does.

[ivarne]

Or even worse; Why did my huge formatting clanup create bugs that were only discovered in production, when I changed 2^-1 to 2 ^ -1 for clarity?

[nalimilan]

If the error message was explicit enough, I guess raising an error wouldn't be a problem.

[andreasnoack]

Isn't the obvious analogue to / to always return Float64 from ^ and have pow for integers similarly to div? I'm not proposing this. Just wondered why it wasn't mentioned here.

[StefanKarpinski]

The trouble with that is that one wants x^2 to be equivalent to x*x. If anything, I would do it the other way around and have pow be the version that gives a Float64 even for integer arguments.

[bsxfan]

As the proposer of the new operator ^- , I admit that it would lead to more confusion and more problems than it would solve, so let's forget about that one.

However, I do think the error-message for negative exponents of integer arguments could be more helpful. Let me propose that we do not look at integer exponentiation in islolation. There are other functions with similar behaviour, for example log() and sqrt(), which in other languages would automatically map negative real arguments to complex outputs, but in Julia give DomainError.

I get:

julia> log(-1)
ERROR: DomainError
log will only return a complex result if called with a complex argument.
try log(complex(x))
in log at math.jl:125

julia> sqrt(-1)
ERROR: DomainError
sqrt will only return a complex result if called with a complex argument.
try sqrt(complex(x))
in sqrt at math.jl:132

julia> 2^(-1)
ERROR: DomainError
in power_by_squaring at intfuncs.jl:60
in ^ at intfuncs.jl:84

It would be nice if this one could give a similar message to the ones given by log() and sqrt().

[ivarne]

+1 for a helpful message. I'll fix that.

[miguelbarao]

I believe it all boils down to the intended audience of the Julia language. I think most people coming from physical sciences and engineering disciplines (like me) will expect 2^-1 to produce a float just because floats are the most frequent used types, the most used scientific languages have that behavior (matlab/octave/python2/python3/R) and out of mathematical consistency.

I would expect the behavior to be:

• 2^3 multiplies 2 three times
• 2^-3 multiplies 2 three times (like above), then take the inverse: 1/2^3 = 1/8 = 0.125.
• 2.0^-3 multiplies 2.0 three times, then take the inverse: 1/2.0^3 = 1/8.0 = 0.125.
• 2^1.2 and 2.0^1.2 slow version of exponentiation.

Making the parallel with python, they replaced the integer division 1/2=0 with the new operator 1//2=0 in version 3, so that 1/2=0.5 is always a float. Inspired on this I had proposed that 2^-3 would always produce a float, while a new operator (e.g. 2^^-3) would always produce an integer.

Not sure about this, but another route could be to ensure type stability by distinguishing between natural and integer numbers (unsigned/signed)? A natural number would be promoted to integer just like an integer can be promoted to a float.

• 2^3 has natural arguments, therefore the result is a natural number.

• (-2)^3 has an integer base, therefore the result is integer.

• 1.2^3 has a float base, therefore the result is a float.

• 2^-3 has an integer exponent, so the result is a float...

  This way the results depend on types and not values of the arguments and it seems consistent with 2^-3=1/2^3. Another way of seeing it is 2^(x-y) = 2^x / 2^y with natural x and y. Presently this is not true for some combinations of x and y.

There are other places where the current way of doing exponentials can lead to problems:

• (-1)^2.0 gives 1.0
• (-1)^2.00000001 produces DomainError in ^ at promotion.jl:162

When the exponent is the result of a computation with error accumulation, this does not work. Example:

• (-1)^(0.1+0.2+0.3+0.4) equals -1.0
• (-1)^(0.4+0.3+0.2+0.1) gives DomainError in ^ at promotion.jl:162

I understand that it is not easy to draw a line. Sometimes we want a domain error, other times we want a complex number, but I feel there is much inconsistency here.

I'm currently considering doing things like
^(x::Int64,n::Int64) = ^(float(x),n)
to run my code, depending on the kind of problems considered.
Are there unforeseen side effects of redefining this operator?

@ivarne the error message in (-1)^0.5 also needs a fix.

[johnmyleswhite]

FWIW, I don't really think there's any inconsistency here. I think there's a tension caused by Julia using integers, when most of the other languages (matlab/octave/R) use floating numbers for everything and invest a lot of energy into pretending that floating point numbers are integers. There's clearly a reason that scientific languages have done that, but I think Julia is intentionally encouraging users to work with integers more often than traditional scientific languages have.

[toivoh]

@miguelbarao: The methods of ^ are shared globally (it's actually ^ from the Base module), so adding that method definition can poison other code in other modules that you run.

There are ways to make a ^ function that is local to your module, but it's a bit roundabout because operators are imported by default.

[timholy]

Although "try it and see" might be relevant advice here.

[miguelbarao]

@toivoh Thank you. Redefining ^ propagates to 2.^-[1,2,3] which is exactly what I intended.
@johnmyleswhite How is 1/0 ok but not 1^-1? This is the kind of consistency that I'm talking about. It just doesn't make sense.

[timholy]

I do think it's worth considering whether ^(x::Integer,n::Signed) should promote x to a Float64.

[toivoh]

@miguelbarao: Yes, the downside to that is that you might break other people's code that relies on the current behavior for performance or correctness.

[johnmyleswhite]

@miguelbarao The "consistent" behavior you want is easily attainable: use floating point numbers for all operations. The result type of an operation is always something you have to learn about for each operation individually, so it's not at all surprising that division and exponentiation should produce different results. In modular arithmetic, for example, integer exponentiation of the sort Julia provides is trivially definable for moduluses, whereas the behavior of division is something that depends strongly on the modulus and requires some subtlety to define meaningfully. This is also true of matrices: exponentiation to a positive integer is trivial to define, but division is much harder. So I don't find contrasting 1/2 with 1^2 convincing: if anything, I find it to be strong evidence in favor of the behavior you dislike.

@timholy's point is a good one, but it involves a change that's potentially undesirable: it makes exponentiation and multiplication more dissimilar because you can't replace x^2 with x*x. I personally think that's much worse than throwing a domain error when given a non-positive integer exponent.

[miguelbarao]

Lets put it in another way (user point of view, like me):

• Domain: natural^natural=natural, integer^integer=rational. Raising to a negative integer power is defined as b^-n = 1/b^n in mathematics. Negative powers are in the domain of integer exponentiation, contrary to the error message in Julia.
• Performance: What is exactly the difference in performance?

julia> x = rand(1:62, 10000000)   # ten million random integers between 1 and 62 (to avoid overflow)
julia> @elapsed 2.^x              # integer base, integer exponent
0.300111225        <- average of several runs, discarding the first and ocasional outliers with higher times

Now with floats:

julia> @elapsed float(2).^x              # float base, integer exponent
0.270533889        <- average...

So floats are actually faster. It could be argued that integer powers are usually small numbers. Let's see: raising only to very low powers, say 2 to 5:

julia> x = rand(2:5, 10000000)
julia> @elapsed 2.^x
0.142693083

So it's nearly double the speed for small exponents. We could call it a draw performance-wise. Am I doing something wrong?
(Core I7 2.4GHz MacBook Pro, OSX)

• Types: 2^(a-b) has always a different type than 2^a / 2^b, and it may not even work in the former case.
• Mathematical consistency: 2^(a+b) is different from 2^a * 2^b if one of them is negative but their sum is still positive (e.g. a=2, b=-1). We may get an integer in the first case and a domain error in the second.

If all this doesn't make a point (and there are still other issues in the exponentiation), I guess Julia is just not for my mindset and should look elsewhere. Really a pity because Julia has a lot of potential and I really admire your work.
@johnmyleswhite All the other operators are fine +, -, *, /. The issue is just with ^.

[johnmyleswhite]

@miguelbarao: Why not just use floating point numbers all the time? That's clearly what you want and is already available in Julia without any changes being made. Instead of writing 2, you can write 2.0 and get exactly the behavior you're asking for.

[tkelman]

That's actually a very interesting point. What if integer^integer was made to always return a rational? Correct me if I'm wrong but that should be type-stable, the question is how much overhead do rationals introduce right now.

[miguelbarao]

@johnmyleswhite Because I'm promoting the use of Julia with my students and whenever I have a problem like "Consider a bag containing N balls, what is the probability of getting the same ball M times", I get these answers:
(1/N)^M
1/N^M
N^-M
and 1/3 of them banging their heads on why it doesn't work.
Since N and M are naturally integers, nobody remembers to do a float(N)^-M or a N^float(-M) which is difficult to explain!
@tkelman Agree. Why not ^ gives float like / does, ^^ gives rational like // does, and pow() gives an integer like div() does. Easy to remember.

[StefanKarpinski]

If all this doesn't make a point (and there are still other issues in the exponentiation), I guess Julia is just not for my mindset and should look elsewhere. Really a pity because Julia has a lot of potential and I really admire your work.

I'm always perplexed that a single nit with one of the many tough decisions and tradeoffs that are necessary to create a good general-purpose language can put some people off of something which they otherwise thought was great. Usually it's overflow of 10^100 or 1-based indexing. I have to wonder what magical language has managed to pass this gauntlet, making no decisions that didn't please.

We should definitely fix the error message and then be done with this. Even students (who are notoriously unable to learn new things) should be able to handle the current behavior if they're given a nice clear error message when it happens.

[grayclhn]

@miguelbarao, it might be worth putting together a "Finesse.jl" package for your students that throws specific warnings for errors that students frequently make, and either auto-converts or throws an error (depending on the level of the class, your goals, etc). The documentation is pretty good at directions for setting up packages, and the users mailing list can help as well.

[timholy]

There's more breath being wasted on the conversation than in just fixing the problem. See #8784.

[StefanKarpinski]

👏 @timholy

[timholy]

@miguelbarao, my adding the error doesn't necessarily mean that no further changes are possible. But you're going to have to do the work (this is open source, after all) and address the concerns. If you go with a floating-point base, you'll need to address the one raised by @johnmyleswhite, how one explains for x = 10^23, why x^2 != x*x. If you go with Rational, then you need to figure out the performance implications.

[miguelbarao]

@timholy 'x = 10^23' overflows as Int64 and x*x overflows again. I presume @johnmyleswhite was referring to types x*x and x^2 being different, not values. Thanks anyway.

[timholy]

OK, try this:

julia> (float(2^27+1))^2 == (2^54+2^28+1)
false

julia> (2^27+1)^2 == (2^54+2^28+1)
true

[timholy]

In other words, it's a problem with values, not types. 2 == 2.0 is true.

[nalimilan]

FWIW, when introducing integers Lua decided that ^ would always return a float so that 2^-2 works. Cf. http://www.inf.puc-rio.br/~roberto/talks/ws2014.pdf (p. 19). Their choice may have to do with backward compatibility.

[stevengj]

Note that 2^-2 and similar expressions (with literal integer exponents) will work fine when/if now that #24240 is merged.

[EvergreenTree]

@stevengj After reading this long discussion, I found that ^ is now type-unstable (now that #24240 is merged), just as @StefanKarpinski said in this post in #20240. So, if I understand it correctly, we gave up the principle of type-stability, only for the sake of coding convenience? (in this case, 2^-2.0 is necessary for efficient compilation, 2^-2 will work but leads to inefficiency)

[StefanKarpinski]

Better to ask this kind of question on discourse than here. In short (please continue on discourse), referential transparency that was sacrificed, not type stability: ^ to a literal integer power is different than raising to a variable with the same integer value.