Montgomery Modular Multiplication and Exponentiation in Julia

mme.jl

const N = 32

"""
Montgomery Modular Multiplication.
"""
function MMM(a::UInt32, b::UInt32, n::UInt32)::UInt32
    s::UInt64 = 0
    q::UInt32 = 0
    for i::UInt32 in 0:N-1
        # i-th bit of a shifted to the least significant position
        a_i = (a & (1 << i)) >> i
        q = (s + a_i * b) % 2
        s = (s + (q * n) + (a_i * b)) >> 1
    end
    return if s >= n
        s - n
    else
        s
    end
end

@assert MMM(0x01234567, 0x89abcdef, 0x93849ca7) == 0x93817349
@assert MMM(0x00000001, 0x00000001, 0x00000001) == 0x00000000

"""
Montgomery Modular Multiplication using only bitshifting operators
"""
function MMM_shift(a::UInt32, b::UInt32, n::UInt32)::UInt32
    s::UInt64 = 0
    q::UInt32 = 0
    a_i_b::UInt32 = 0
    for _ in 0:N-1
        # i-th bit of a shifted to the least significant position
        a_i = (a << 31) >> 31
        a = a >> 1
        a_i_b = a_i * b
        q = ((s + a_i_b) % 2)
        s = (s + (q * n) + a_i_b) >> 1
    end
    return if s >= n
        s - n
    else
        s
    end
end

MMM_shift(0x01234567, 0x89abcdef, 0x93849ca7)


function calc_k(n, N)
    convert(UInt32, convert(UInt128, 2) ^ convert(UInt128, 2 * N) % n)
end

const E = 64

"""
Montgomery Modular Exponentiation
"""
function MME(m::UInt32, e::UInt32, n::UInt32)
    k::UInt32 = calc_k(n, N)
    m = MMM(m, k, n)
    r::UInt32 = MMM(0x00000001, k, n)
    for _ in 1:E-1
        e_i = (e << 31) >> 31
        e = e >> 1
        if e_i == 0x00000001
            r = MMM(r, m, n)
        end
        m = MMM(m, m, n)
    end
    return MMM(r, 0x00000001, n)
end

MME(0x01234567, 0x89abcdef, 0x93849ca7) # 0x31a9137c
MMM(0x00000001, 0x8c8d9129, 0x93849ca7)