named functions (du_dx, du_dy, dv_dx, dv_dy) to make it completely transparent that the cauchy-riemann equations are being used, and what the actual check is.

holomorphic.jl

# f(z = x + iy) = u(x, y) + i * v(x, y)
function finite_difference{T <: Complex}(f::Function,
                                         z::T,
                                         alongx::Bool = true,
                                         alongu::Bool = true)
    epsilon = sqrt(eps(max(abs(one(T)), abs(z))))
    if alongx
        zplusdz = z + epsilon
    else
        zplusdz = z + epsilon * im
    end
    if alongu
        return real(f(zplusdz) - f(z)) / epsilon
    else
        return imag(f(zplusdz) - f(z)) / epsilon
    end
end

function du_dx(f::Function, z::Complex)
    return finite_difference(f, z, true, true)
end

function du_dy(f::Function, z::Complex)
    return finite_difference(f, z, false, true)
end

function dv_dx(f::Function, z::Complex)
    return finite_difference(f, z, true, false)
end

function dv_dy(f::Function, z::Complex)
    return finite_difference(f, z, false, false)
end

function isholomorphic(f::Function, z::Complex)
    cond1 = isapprox(du_dx(f, z), dv_dy(f, z))
    cond2 = isapprox(du_dy(f, z),-dv_dx(f, z))
    return cond1 && cond2
end

function isholomorphic(f::Function, N::Integer = 1000)
    res = 0
    for i in 1:N
        z = randn() + randn() * im
        res += int(isholomorphic(f, z))
    end
    return res / N
end

isholomorphic(sin) # => 0.979
isholomorphic(cos) # => 0.992
isholomorphic(abs) # => 0.0