Should large exponentials not yield Inf?

General Usage
, ,
1

I am wondering why some large exponentials yield wrong values, specifically always 1 or 0

julia> 123^(456^789)
1

julia> (0.1)^(456^789)
1.0

julia> 456^789
0

This results in wrong values inside functions

julia> atan((0.1)^(456^789))
0.7853981633974483

julia> atan(1)
0.7853981633974483

I’d expect Inf, like

julia> exp(4567891231242141)
Inf

also note that

julia> 456^789.1
Inf

so the result comes from Integer exponentiation I think.

2
5 Likes
3

To summarize that very long discussion:

No, this is expected behavior. Julia integers are machine integers, ^ on integers does integer exponentiation. Checking for overflow & throwing an error/giving a warning would prevent a lot of optimizations like vectorization to take place. If you want to do floating point computations (with possibly Inf as a result), use floating point literals, which come at the cost of precision depending on the size of the number.

Note that this doesn’t happen “in the atan function” - ^ is computed first, and the result is then passed into atan. atan does not know where the number came from. It’s equivalent to writing

a = 456^789
b = 0.1^a
c = atan(b)
2 Likes
4

So this is where the difference lies:

julia> 456^789
0

julia> 456.0^789
Inf
2 Likes
5

And Inf is not defined for integers anyways… Only for floating point numbers…

2 Likes
6

This is right(!) in some technical (modular) sense. Modular arithmetic is the default for machine integers, and they can’t be Inf. So Int^Int, CAN’T result in Inf, if, given that it results in Int. I would argue it should result in Float64, similar to division, for the same reasons (at least for signed integers), when the power is (possibly) negative. Right now, power returns Int for purity reasons, that arguably do not apply to it, only the other operators.

I’m not sure a lot of code relies on integer exponentiation, but it it does then plan b could be to check for overflow and throw an error, and I’m not sure it will be slow or prevent optimization. Because I think there’s no such assembly instruction.

7

There has been some serious discussion of doing an overflow check for integer^integer non-literal powers, though I can’t find it at the moment (is there an issue)? Rationale:

  • integer^integer powers rarely seem to be performance sensitive in real applications, and are expensive enough that an overflow check might not add too much overhead anyway. (We could omit the overflow case for small literal powers like ^2 and ^3, which are inlined via literal_pow.)
  • We already check for overflow in factorial(integer), which is somewhat analogous.
  • integer^integer overflow is an extremely common source of user confusion.
7 Likes
8

Yes, I was searching for it in the very long discussion, but gave up. I believe Karpinski or some such was open to it, but I just opened an issue on rather the more general Float64 result.

9

Sorry for reviving this issue, I should have checked more thoroughly. I think I now get where the problem’s coming from, just have to be more careful with the modular nature of integer exponentials on a finite memory computer.

2 Likes
10

I guess an error message would do it for starters because I could only imagine hunting for such a problem in a larger codebase, that would certainly not be pleasant. Admittedly, I don’t know how easy such an implementation would be in a general case or possible at all, where values aren’t given beforehand. Dynamic checks might do the trick but those are supposedly not the best solution performance-wise. The pow() from C uses Floats and doubles, so that could serve as a potential solution, yet from a mathematical perspective the integers are indeed closed under exponentiaton so the result being an integer again also seems reasonable.This was discussed previously so there’s no point arguing about that again. I understand that each solution has its own ups and downsides, what I just find weird and confusing (especially for starters coming from python) is the fact that integer exponentiation can overflow to 0, while the floating variant does result in Inf.

11

That’s one of the possible design choices for integer exponentiation and frequently chosen when speed gets the highest priority. (Strictly speaking it doesn’t overflow specifically to 0 but the overflowing result is computed modulo 2^N, which for many numbers end up at 0.)

And you don’t have to leave Python to find examples of this design choice. Try

import numpy
numpy.int64(456) ** 789

or

from numba import jit

@jit
def f():
    return 456 ** 789
4 Likes
12

You imply modulo 2^64 (or 2^32) but that’s only correct for UInt. For Int it’s modulo 2^63 or even worse modulo 2^31, because of the sign bit.

My question, does anyone in practice rely on/want modulo 2^63 for Int (as opposed to for UInt)? All the number-theoretic code I could imagine would rely on modulo 2^64 or higher. Also if you rely on even 2^63 big enough, and not overflowing, then your code is wrong, and much more likely to overflow on 32-bit, where typemax(Int) is 2147483647. On of the reasons I opened an issue arguing for the correct Float64 result (currently has 3 downvotes, no upvotes… so less likely to happen).

13

A precise formulation for signed integers is that the overflowed result y relates to the non-overflowed result z such that mod(y - z, 2^N) = 0, where N is the number of bits in y. Another formulation can be found in the docstring for x % T where T is a type.

Example:

julia> 119^3
1685159

julia> Int8(119)^3
-89

julia> (119^3) % Int8
-89

julia> mod(-89 - 1685159, 2^8)
0

Nothing is computed modulo 2^63 or 2^31. Signed integers overflow earlier than their unsigned counterparts and positive numbers can overflow into negative numbers but you don’t lose a bit in the range.

14

Note that unsigned 64 bit integers have [0, 2^64 - 1] as their domain, while signed 64 bit integers have [-2^63, 2^63 - 1]. The 2^31 would be for 32 bit integers - you only “lose” a single bit, so a factor of 2, not half of your bits.

15

My point with mentioning 2^31 is that you do have half the bits on 32-bit systems, where Int is Int32. Yes, they’re a dying breed, but if you write portable code the not going to Inf (or overflow) is worse there.

Yes, technically it’s 64 or 32 bits, but rather useless fact to most, because they do not wnat (modulo math), negative numbers, thus in practice 63 correct or only 31. But you have no easy way to know if those bits are correct, unless the sign bit is set…

16

The operation is exactly the same for signed or unsigned, the only difference is the interpretation of the values. Both UIntN and IntN are rings modulo 2^N and they use the exact same operations even, we just choose a different set of “standard representatives” for the modular equivalence classes: for UIntN we choose 0 through 2^N-1 whereas for IntN we choose -2^(N-1) through 2^(N-1)-1.

3 Likes