square_root_ceil(n) = ceil(Int, sqrt(n))

function factor_fermat0(n)
    a = square_root_ceil(n)
    while true
        d = a*a-n
        b = square_root_ceil(d)
        if b^2 == d
            return a-b, a+b
        end
        a += 1
    end
end

 function factor_fermat1(n)
    a = square_root_ceil(n)
    b = square_root_ceil(a^2-n)

    while a^2-b^2 != n
        while a^2-b^2 < n; a += 1; end
        while a^2-b^2 > n; b += 1; end
    end

    return a-b, a+b
end

function factor_fermat2(n)
    a = square_root_ceil(n)
    a2 = a^2
    b = square_root_ceil(a^2-n)
    b2 = b^2
    while a2-b2 != n
        while a2-b2 < n; a2, a = a2+a+a+1, a+1; end
        while a2-b2 > n; b2, b = b2+b+b+1, b+1; end
    end

    return a-b, a+b
end

function factor_fermat3(n)
    a = square_root_ceil(n)
    a2 = a^2
    c = 2*a+1
    b = square_root_ceil(a^2-n)
    b2 = b^2
    d = 2*b+1
    while a2-b2 != n
        while a2-b2 < n; a2, c = a2+c, c+2 end
        while a2-b2 > n; b2, d = b2+d, d+2 end
    end

    return (c-d)÷2, (c+d)÷2-1
end

function factor_fermat4(n)
    a = square_root_ceil(n)
    a2 = a^2
    c = 2*a+1
    b = square_root_ceil(a^2-n)
    b2 = b^2
    d = 2*b+1
    while a2-b2 != n
        while a2-b2 < n
            a2 += c
            c += 2
        end
        while a2-b2 > n
            b2 += d
            d += 2
        end
    end

    return (c-d)÷2, (c+d)÷2-1
end

# import Pkg; Pkg.add("BenchmarkTools")
using BenchmarkTools
for f in (factor_fermat0, factor_fermat1, factor_fermat2,
         factor_fermat3, factor_fermat4)
    @btime $f(4736906887651)
end