diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 558d55c118b34..ea35da0910736 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -218,7 +218,7 @@ gcdx(a::T, b::T) where T<:Real = throw(MethodError(gcdx, (a,b))) # multiplicative inverse of n mod m, error if none """ - invmod(n, m) + invmod(n::Integer, m::Integer) Take the inverse of `n` modulo `m`: `y` such that ``n y = 1 \\pmod m``, and ``div(y,m) = 0``. This will throw an error if ``m = 0``, or if @@ -257,6 +257,43 @@ function invmod(n::Integer, m::Integer) return mod(x, m) end +""" + invmod(n::Integer, T) where {T <: Base.BitInteger} + invmod(n::T) where {T <: Base.BitInteger} + +Compute the modular inverse of `n` in the integer ring of type `T`, i.e. modulo +`2^N` where `N = 8*sizeof(T)` (e.g. `N = 32` for `Int32`). In other words these +methods satisfy the following identities: +``` +n * invmod(n) == 1 +(n * invmod(n, T)) % T == 1 +(n % T) * invmod(n, T) == 1 +``` +Note that `*` here is modular multiplication in the integer ring, `T`. + +Specifying the modulus implied by an integer type as an explicit value is often +inconvenient since the modulus is by definition too big to be represented by the +type. + +The modular inverse is computed much more efficiently than the general case +using the algorithm described in https://arxiv.org/pdf/2204.04342.pdf. + +!!! compat "Julia 1.11" + The `invmod(n)` and `invmod(n, T)` methods require Julia 1.11 or later. +""" +invmod(n::Integer, ::Type{T}) where {T<:BitInteger} = invmod(n % T) + +function invmod(n::T) where {T<:BitInteger} + isodd(n) || throw(DomainError(n, "Argument must be odd.")) + x = (3*n ⊻ 2) % T + y = (1 - n*x) % T + for _ = 1:trailing_zeros(2*sizeof(T)) + x *= y + true + y *= y + end + return x +end + # ^ for any x supporting * to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x) @noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p, LazyString( diff --git a/test/intfuncs.jl b/test/intfuncs.jl index ceaac235a3da9..7fa457e63cd4a 100644 --- a/test/intfuncs.jl +++ b/test/intfuncs.jl @@ -221,7 +221,7 @@ end @test_throws MethodError gcdx(MyOtherRational(2//3), MyOtherRational(3//4)) end -@testset "invmod" begin +@testset "invmod(n, m)" begin @test invmod(6, 31) === 26 @test invmod(-1, 3) === 2 @test invmod(1, -3) === -2 @@ -256,6 +256,37 @@ end end end +@testset "invmod(n)" begin + for T in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128) + if sizeof(T) ≤ 2 + # test full domain for small types + for a = typemin(T)+true:T(2):typemax(T) + b = invmod(a) + @test a*b == 1 + end + else + # test random sample for large types + for _ = 1:2^12 + a = rand(T) | true + b = invmod(a) + @test a*b == 1 + end + end + end +end + +@testset "invmod(n, T)" begin + for S in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128), + T in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128) + for _ = 1:2^8 + a = rand(S) | true + b = invmod(a, T) + @test (a * b) % T == 1 + @test (a % T) * b == 1 + end + end +end + @testset "powermod" begin @test powermod(2, 3, 5) == 3 @test powermod(2, 3, -5) == -2