The issue is that the Python code uses arbitrary precision while Julia code uses 64 bits. Python may be slower because of that. One way to get rid of this discrepancy is to compile with typed Cython, or with Numba, which is what I did. The other way would be to benchmark Julia with BigInt. Either way would be fair IMHO.