Curious solution

curious.jl

#!/usr/bin/env julia

# Adapted from https://github.com/pablocelayes/rsa-wiener-attack

using ContinuedFractions # Note: this needs my fork of ContinuedFractions.jl

# fast perfect square test (from project euler tools)
# returns -1 if n is not a perfect square
function perfect_square(n::Integer)
    h = n & 0xF  # last four bits
    h > 9 && return -1
    if h == 9 || h <= 1 || h == 4
        t = isqrt(n)
        t*t == n && return t
        return -1
    end
    -1
end

function find_d{T<:Integer}(e::T, n::T)
    for (i, rat) in ContinuedFraction(e//n) |> convergents |> enumerate
        k, d = rat.num, rat.den
        if (k != 0) && ((e*d - 1) % k == 0)
            ϕ = div((e*d - 1), k)
            s = n - ϕ + 1
            # check if the equation x^2 - s*x + n = 0 has integer roots
            discr = s*s - 4n
            discr < 0 && continue
            t = perfect_square(discr)
            t < 0 && continue # not a perfect square
            iseven(s+t) && return d # :)
        end
    end
    0 # default
end

# This is ugly, I know
function recover_primes{T<:Integer}(d::T, e::T, N::T)
    k = d*e - 1
    isodd(k) && return (0, 0) # not found
    t = 0
    while iseven(k)
        t += 1
        k >>= 1
    end
    x, y = 0, 0
    for i in 1:100
        x_one, x_n = false, false
        g = rand(0:N-1)
        y = powermod(g, k, N)
        (y == 1 || y == N-1) && continue
        for j in 1:t-1
            x = y*y % N
            if x == 1
                x_one = true
                break
            end
            if x == N-1
                x_n = true
                break
            end
            y = x
        end
        x_one && break
        x_n && continue
        x = y*y % N
        x == 1 && break
    end
    p = gcd(y-1, N)
    q = div(N, p)
    p, q
end

function case(_)
    N = readline() |> BigInt
    e = readline() |> BigInt
    c = readline() |> BigInt
    d = find_d(e, N)
    if d == 0
        println("Could not find d.")
        return
    end
    p, q = recover_primes(d, e, N)
    println("d = $d")
    println("p = $p")
    println("q = $q")
    m = powermod(c, d, N) # decrypt message
    bytes = base(16, m) |> hex2bytes
    msg = map(Char, bytes) |> join
    println("Message: $msg")
end

function main()
    VERSION < v"0.4-" && error("Julia 0.3.x overflows the stack here :(")
    map(case, 1:101)
end
main()