Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Exponentiation by a negative power throws a DomainError #3024

Closed
bensadeghi opened this issue May 5, 2013 · 50 comments · Fixed by #24240
Closed

Exponentiation by a negative power throws a DomainError #3024

bensadeghi opened this issue May 5, 2013 · 50 comments · Fixed by #24240

Comments

@bensadeghi
Copy link
Member

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
Copy link
Member

pao commented May 6, 2013

Did you intend to close this?

@JeffBezanson
Copy link
Member

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

@bensadeghi
Copy link
Member Author

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

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

@diegozea
Copy link
Contributor

diegozea commented May 6, 2013

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
Copy link
Member

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
Copy link

bsxfan commented May 9, 2013

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
Copy link

bsxfan commented May 9, 2013

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

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

@StefanKarpinski
Copy link
Member

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
Copy link

bsxfan commented May 9, 2013

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
Copy link

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

@jiahao
Copy link
Member

jiahao commented Dec 15, 2013

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

@miguelbarao
Copy link

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
Copy link
Member

@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
Copy link
Member

timholy commented Oct 20, 2014

@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
Copy link
Member

ivarne commented Oct 20, 2014

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
Copy link
Member

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

What is this, JavaScript? :-P

@StefanKarpinski
Copy link
Member

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

@pao
Copy link
Member

pao commented Oct 20, 2014

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

@StefanKarpinski
Copy link
Member

That it certainly does.

@ivarne
Copy link
Member

ivarne commented Oct 20, 2014

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
Copy link
Member

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

@andreasnoack
Copy link
Member

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
Copy link
Member

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
Copy link

bsxfan commented Oct 21, 2014

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
Copy link
Member

ivarne commented Oct 21, 2014

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

@miguelbarao
Copy link

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
Copy link
Member

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
Copy link
Contributor

toivoh commented Oct 22, 2014

@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
Copy link
Member

timholy commented Oct 22, 2014

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

@miguelbarao
Copy link

@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
Copy link
Member

timholy commented Oct 22, 2014

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

@toivoh
Copy link
Contributor

toivoh commented Oct 22, 2014

@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
Copy link
Member

@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
Copy link

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
Copy link
Member

@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
Copy link
Contributor

tkelman commented Oct 22, 2014

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
Copy link

@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
Copy link
Member

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
Copy link
Contributor

@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
Copy link
Member

timholy commented Oct 23, 2014

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

@StefanKarpinski
Copy link
Member

👏 @timholy

@timholy
Copy link
Member

timholy commented Oct 24, 2014

@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.

StefanKarpinski added a commit that referenced this issue Oct 24, 2014
Raise informative error for exponentiation to negative power (see #3024)
@miguelbarao
Copy link

@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
Copy link
Member

timholy commented Oct 26, 2014

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
Copy link
Member

timholy commented Oct 26, 2014

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

timholy added a commit that referenced this issue Nov 7, 2014
(cherry picked from commit 79c9dfc)
make #8784 work on Mac
(cherry picked from commit d2252d0)
@nalimilan
Copy link
Member

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
Copy link
Member

stevengj commented Oct 20, 2017

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

@EvergreenTree
Copy link

@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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet