improve performance of powermod (#59930)

based off of and closes #54866

differences to #54866 are:

* checks `iszero(r)` in `mod(a, ::SignedMultiplicativeInverse)` for
cases like `mod(0, multiplicativeinverse(-5))`
* catches some cases of mismatched signed-ness of `x, m` that were
failing
* added a very very rough heuristic `(p > 2sizeof(mm))` so we don't
bother computing the mi when the power is very small
* switched to LSB multiplication loop from MSB

benchmark:
```
using BenchmarkTools

function setup()
    x,p,m = rand(NTuple{3, Int})
    while iszero(m) || ((p < 0) && !isone(gcd(x, m)))
        x,p,m = rand(NTuple{3, Int})
    end
    return (x, p, m)
end

@benchmark powermod(x, p, m) setup=((x,p,m) = setup())


#master
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
 Range (min … max):  1.304 μs …   8.938 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.779 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.803 μs ± 270.848 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

               ▂▄▅████▇▆▄▁▁                                    
  ▁▁▁▁▁▁▁▂▂▃▅▅▇█████████████▆▅▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  1.3 μs          Histogram: frequency by time        2.76 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.



#PR
julia> @benchmark powermod(x, p, m) setup=((x,p,m) = setup())
BenchmarkTools.Trial: 10000 samples with 135 evaluations per sample.
 Range (min … max):  667.896 ns …  1.267 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     854.015 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   886.407 ns ± 74.610 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                        ▂▃▆██▅▂                                 
  ▂▂▁▁▂▂▂▂▂▂▂▂▂▂▃▃▃▄▄▅▆▇████████▆▅▄▃▃▃▃▃▄▄▅▆▆▇▇████▇▆▆▅▄▃▃▃▃▂▂ ▄
  668 ns          Histogram: frequency by time         1.06 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.
```

Co-authored-by: Chris Elrod <[email protected]>

Author: adienes

base/intfuncs.jl

@@ -535,20 +535,47 @@ function powermod(x::Integer, p::Integer, m::T) where T<:Integer
         end
     end
     (m == 1 || m == -1) && return zero(m)
-    b = oftype(m,mod(x,m))  # this also checks for divide by zero
-
-    t = prevpow(2, p)
-    r = 1
-    while true
-        if p >= t
-            r = mod(widemul(r,b),m)
-            p -= t
+
+    mm = uabs(m)
+    rr = one(mm)
+    bb = oftype(mm, mod(x, mm))
+
+    # legal && profitable
+    if _powermod_mi_legal(mm) && (p > 2sizeof(mm))
+        if bb == 0
+            rr = zero(mm)
+        else
+            mis = MultiplicativeInverses.multiplicativeinverse(mm)
+            Base.@assume_effects :terminates_locally while true
+                if (p & 1) != 0
+                    rr = mod(rr * bb, mis)
+                end
+                p >>= 1
+                p == 0 && break
+                bb = mod(bb * bb, mis)
+            end
+        end
+    else
+        if bb == 0
+            rr = zero(mm)
+        else
+            Base.@assume_effects :terminates_locally while true
+                if (p & 1) != 0
+                    rr = oftype(mm, mod(widemul(rr, bb), mm))
+                end
+                p >>= 1
+                p == 0 && break
+                bb = oftype(mm, mod(widemul(bb, bb), mm))
+            end
         end
-        t >>>= 1
-        t <= 0 && break
-        r = mod(widemul(r,r),m)
     end
-    return r
+    r = oftype(m, rr)
+    return (iszero(r) || (m > 0)) ? r : r + m
+end
+
+_powermod_mi_legal(::Integer) = false
+function _powermod_mi_legal(mm::T) where {T<:Unsigned}
+    return Base.hastypemax(T) && (mm <= (typemax(T) >> (sizeof(T) << 2)))
 end
 
 # optimization: promote the modulus m to BigInt only once (cf. widemul in generic powermod above)

base/multinverses.jl

@@ -2,7 +2,7 @@
 
 module MultiplicativeInverses
 
-import Base: div, divrem, mul_hi, rem, unsigned
+import Base: div, divrem, mul_hi, rem, unsigned, mod
 using  Base: IndexLinear, IndexCartesian, tail
 export multiplicativeinverse
 
@@ -153,6 +153,13 @@ function divrem(a::T, b::MultiplicativeInverse{T}) where T
     (d, a - d*b.divisor)
 end
 
+mod(a::T, b::UnsignedMultiplicativeInverse{T}) where {T} = rem(a, b)
+
+function mod(a::T, b::SignedMultiplicativeInverse{T}) where {T}
+    r = rem(a, b)
+    return (iszero(r) || signbit(r) == signbit(b.divisor)) ? r : r + b.divisor
+end
+
 multiplicativeinverse(x::Signed) = SignedMultiplicativeInverse(x)
 multiplicativeinverse(x::Unsigned) = UnsignedMultiplicativeInverse(x)
 

test/intfuncs.jl

@@ -385,6 +385,10 @@ end
 
     @test powermod(-3, 0x80, 7) === 2
     @test powermod(0x03, 0x80, 0x07) === 0x02
+
+    @test powermod(511, 1, 0x00000021) === 0x00000010
+    @test powermod(Int8(-1), 0xff, Int8(33)) === Int8(32)
+    @test powermod(0, 10, -5) === 0
 end
 
 @testset "nextpow/prevpow" begin