Modular multiplication without overflow

General Usage Performance
4

Translating the (completely opaque and unverified) wikipedia code to an equally dodgy julia version:

Listing: `mul_mod(a,b,c)` 
function mul_mod(a::UInt128, b::UInt128, m::Integer)

    magic1 = (UInt128(0xFFFFFFFF) << 32)
    magic2 = UInt128(0x8000000000000000)

    if iszero(((a | b) & magic1))
        return (a * b) % m
    end

    d = zero(UInt128)
    mp2 = m >> 1
    
    if a >= m; a %= m; end
    if b >= m; b %= m; end

    for _ in 1:64
        (d > mp2) ? d = ((d << 1) - m) : d = (d << 1)
        if !iszero(a & magic2)
            d += b
        end
        if d >= m
            d -= m 
        end
        a <<= 1
    end

    return d
end

function mul_mod(a::Integer, b::Integer, m::Integer)
    sa, sb = UInt128.(unsigned.((a,b)))
    return mul_mod(sa, sb, m)
end

“Testing”:

julia> a, b, c = (4856388669903036, 6493747681980967, 3860134847225580)
(4856388669903036, 6493747681980967, 3860134847225580)

julia> (big(a) * big(b)) % big(c)  ## Correct, but using BigInt
3059575088694852

julia> (a * b) % c ## Incorrect
2569083244915888

julia> ((a % c) * (b % c)) % c  ## Also incorrect due to overflow of multiplication
226995574897440

julia> mul_mod(a, b, c) |> Int128  ## Correct, without BigInt conversion
3059575088694852
2 Likes