WIP: Introduce Overflow Checking for ^(::Integer, ::Integer)

[Keno]

I haven't updated the tests yet, mostly because we need to get a consensus on what to do.
The thinking behind this is that by far the most common case of people running into overflows is trying to do things like 10^x, when they meant 10.0^x. Since this functionality is very well isolated, we can add checking without generally killing performance. There's a couple of details to work out though:

• What to do about small powers ^2, ^3, etc. @StefanKarpinski expressed a preference for those to keep behaving like x*x, x*x*x.
• What to do about literals. Do we keep referential transparency (e.g. do the literal ^2, ^3 no check for overflow, but the non-literal one does).
• Are we ok with this? Initial performance numbers look like about a 10-20% performance penalty for this change.

cc @JeffBezanson @StefanKarpinski @stevengj @timholy

[stevengj]

I don't like the idea of slowing down a really primitive operation like ^, particular when the overflow checking is only done sporadically.

[Keno]

Is integer pow really used that often though? The range of usable inputs is actually fairly small (particularly for say Int32).

[TotalVerb]

It seems a little inconsistent to me to only check overflow in this function. Wouldn't it be better to have a flag/environment/module where all functions are overflowed check?

[Keno]

It would be, but we can't implement that in an efficient manner. This is supposed to catch a case where this happens very frequently that's implementable without terrible performance regressions.

[timholy]

If it's only a 10-20% slowdown, and if it seems unlikely that future changes in CPUs will significantly widen the gap, then I'm certainly willing to consider it. I'd think we'd want a consistent policy with regards to ^, however, even for small integer powers. (Meaning that, in contrast to @StefanKarpinski's suggestion above, I'd say if we have it for any powers of integers we should have it for all.)

[tkelman]

probably want an instance, not the type

[Keno]

yes, already fixed locally. I just figured I'd start the policy discussion with some code and numbers ;).

[Keno]

I'm happy to make it always checked.

[timholy]

I wonder if we could use @fastmath to skip the check? I'm mostly thinking of @simd, where any branch in the inner loop prevents vectorization.

[PallHaraldsson]

"where any branch in the inner loop prevents vectorization."

Do you need that? I'm not sure with LLVM, but in assembly can't you get a flag; and can you then OR them together, and check out of the loop?

[Or is there a fast way to rule out overflow, something like (64-clz(a))b(64-clz(b)) :

https://en.wikipedia.org/wiki/Find_first_set#Hardware_support ]

[Keno]

Does anybody have an example of real code where this is performance critical? It would be good to look at some real world examples, because I suspect in a lot of cases we might be able to hoist the check anyway.

[oscardssmith]

Adding a triage label to see if we want to do this in light of discussion on discourse https://discourse.julialang.org/t/discussion-about-integer-overflow/69627

[stevengj]

Probably the details of the implementation need to be revisited, but I've grown to be in favor of the overall idea.
I'm inclined to think that @Keno is right and that Int^Int is rarely performance critical in real code, so that overflow checking is worth the overhead, at least in ^(::IntXX, ::IntXX).

(To be conservative, we could keep small literal powers x^2 == x*x and x^3 == x*x*x as-is for now: those are the most likely to be performance critical, and the least likely to incur inadvertent overflow.)

[oscardssmith]

Yeah, I have an approach that I think will be easier to implement, and might be faster. The basic idea is to get a lower bound for x^y by using x^y >= (2^floor(log2(x)))^y = 2^(floor(log2(x) * y). Therefore, we can check
(63-leading_zeros(x))*y < 63 && 1<<(63-leading_zeros(x)) < x^y (replace 63 with 31 for Int32). This is an overflow check that is exact, and will 7 instructions if your processor has a fast ctlz_int which is true on Armv8, x86 with BMI1 (Haswell/Jaguar or newer), and Risc-V with "B" extension.

[JeffBezanson]

I'm not sure this is worth it. It just seems kind of complex/random/undisciplined to check for overflow only in this case. As it is, those who want overflow to be checked at least acknowledge that we made a decision and stuck with it 🤷

[oscardssmith]

The theory here is that powers are the one integer operation that is commonly used and commonly overflows with small-ish inputs. Also, with this implementation, the amount of overhead is approximately 0.

[stevengj]

We also throw OverflowError in factorial, which is the integer operation most similar to x^y.

[stevengj]

@oscardssmith, as I understand it, your suggestion would look something like:

function ^(x::IntXX, y::IntXX)
    xʸ = power_by_squaring(x, y)
    nbits = sizeof(x) * 8 - 1
    if !((nbits-leading_zeros(x))*y < nbits && 1<<(nbits-leading_zeros(x)) < xʸ)
        throw(OverflowError("some informative message"))
    end
    return xʸ
end

It doesn't seem quite right, because it throws for 10^200 but not 10^20…

[perrutquist]

Maybe this is stating the obvious, but...

For literal exponents, the overflow check only requires two integer comparisons at runtime, since the bounds on x can be pre-computed. For example:

function literal_pow_overflows(x::Int64, ::Val{p}) where {p}
    p <= 1 && return false
    p >= 64 && return x < -1 || x > 1
    (x < (-3037000499, -2097152, -55108, -6208, -1448, -512, -234, -128, -78, -52, -38, -28, 
        -22, -18, -15, -13, -11, -9, -8, -8, -7, -6, -6, -5, -5, -5, -4, -4, -4, -4, -3, -3, -3, -3, 
        -3, -3, -3, -3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, 
        -2, -2, -2, -2, -2)[p-1] ||
    x > (3037000499, 2097151, 55108, 6208, 1448, 511, 234, 127, 78, 52, 38, 28, 
         22, 18, 15, 13, 11, 9, 8, 7, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 
         3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
         2, 2, 2, 1)[p-1])
end

(The same method could obviously also be applied to non-literal exponents, although the table lookups would probably be too expensive.)

[PallHaraldsson]

I really think we should do something about overflow (since it's a very common criticism of Julia), even just for power, in some form, and I prefer A:

A.
Simply have a^b return Float64 (similar to, and logically the same rationale as, for division), since in general returning an integer is wrong (and power is the only type-unstable primitive operation), consider the previous design a bug (if b is Unsigned, that's an exception, and can keep status quo). It seems to be at worst 3.7% slower, but since a floating-point instruction (as opposed to larger integer code), I think we win it back in real-world-sized non-benchmark code.

Then my proposal here (made for physics code, speed-of-light is redundant):
tonyhffong/Lint.jl#271

We could try this (trivial PR, I'm willing to make it), in 1.9-DEV and PkgEval it.

[For small literal powers, if people want integer results, then a workaround is possible, by changing your code, or if wanted keep same literal power logic and return integers there.]

B.
The code, as is, has many checks, and it could start with if b > 3 and only overflow check for that. It wouldn't be any slower for b <= 3. We could do that for the b is Unsigned case if we bother do "fix" it at all.

[stevengj]

@PallHaraldsson, my sense is people are broadly receptive to this, but we need someone to step up and write an optimized implementation.

[perrutquist]

If this is implemented, then there also needs to be an easy way to opt out from checking when the overflow is intentional. My suggestion would be to let Base.powermod accept a type as the modulus, so that powermod(x,p,typeof(x)) can be used for non-checking x^p.

[LilithHafner]

I think that the idea to catch common overflows is great, but this is technically breaking because 17^29 is currently a valid way to compute 17 to the power of 29 mod Int64. Given the lack of decisive support for this, I think it is not worth slipping this minor breakage into a 1.x release.

Perhaps this could become part of a larger discussion on overflow behavior come 2.0.

[oscardssmith]

triage thinks this is probably too breaking for 1.x

[vtjnash]

Just saw this was implemented and merged in a later PR

[Keno]

I don't think so:

julia> 10^200
0

[oscardssmith]

It's Base.checked_pow in #52849.

[Keno]

That's not what this PR does though. This PR turns on overflow checking by default in ^.

[vtjnash]

This PR attempted to do both things at once it seems, which stalled it. We can change the default later, but the majority of the PR looks like it was implemented now

[PallHaraldsson]

this is technically breaking because 17^29 is currently a valid way to compute 17 to the power of 29 mod Int64.

No it isn't, this is mod Int32 on 32-bit computers. So getting two different answers is problematic in itself for me, unless for in case of calculations where there's actually no truncation (and for portability you need to think of 32-bit). And this is the most dangerous operation. [Technically my argument for widening goes for *, just even better there, then no change of overflow, and <<, + and -, but you want one simple machine instruction for those, and * much less dangerous, and ^ hardly ever speed-critical.]

How about widening always to Int128, to at least have a fast operation (rather than to BigInt...), and then you can mod all you want later (or cast to Float64 what you likely should have done)...?

This is probably a bug (or intended for literal?):

julia> 2^UInt128(2)
4

julia> typeof(ans)
Int64