Pollard-rho algorithm in Julia and GMP C

gistfile1.txt

$ lsb_release -a
No LSB modules are available.
Distributor ID:	Ubuntu
Description:	Ubuntu 22.04.3 LTS
Release:	22.04
Codename:	jammy

$ grep "processor\|model name" /proc/cpuinfo
processor	: 0
model name	: Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz
processor	: 1
model name	: Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz
processor	: 2
model name	: Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz
processor	: 3
model name	: Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz

$ python3 --version
Python 3.10.12

$ ./pypy3.10-v7.3.12-linux64/bin/pypy3 --version
Python 3.10.12 (af44d0b8114cb82c40a07bb9ee9c1ca8a1b3688c, Jun 15 2023, 12:39:27)
[PyPy 7.3.12 with GCC 10.2.1 20210130 (Red Hat 10.2.1-11)]

$ ./julia-1.9.2/bin/julia --version
julia version 1.9.2

$ gcc --version
gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0
Copyright (C) 2021 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ gcc -O2 -o gmp_c pollard_rho.c -lgmp

$ time python3 pollard_rho.py 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m14.730s
user	0m14.717s
sys	0m0.008s

$ time python3 pollard_rho_gmpy2.py 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m7.452s
user	0m7.439s
sys	0m0.009s

$ time ./pypy3.10-v7.3.12-linux64/bin/pypy3 pollard_rho.py 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m14.228s
user	0m14.109s
sys	0m0.068s

$ time ./julia-1.9.2/bin/julia pollard_rho.jl 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m20.532s
user	0m20.452s
sys	0m0.241s

$ time ./julia-1.9.2/bin/julia pollard_rho_opt.jl 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m5.199s
user	0m5.077s
sys	0m0.288s

$ time ./gmp_c 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real	0m3.705s
user	0m3.699s
sys	0m0.005s

pollard_rho.c

#include <stdio.h>
#include <gmp.h>

int pollard_rho(mpz_t *ptr_n, mpz_t *ptr_p)
{
    mpz_t x, y, d;
    int ret;

    mpz_init_set_ui(x, 2);
    mpz_init_set_ui(y, 2);
    mpz_init_set_ui(d, 1);

    while (mpz_cmp_ui(d, 1) == 0) {
        mpz_mul(x, x, x); mpz_add_ui(x, x, 1); mpz_mod(x, x, *ptr_n);
        mpz_mul(y, y, y); mpz_add_ui(y, y, 1); mpz_mod(y, y, *ptr_n);
        mpz_mul(y, y, y); mpz_add_ui(y, y, 1); mpz_mod(y, y, *ptr_n);
        mpz_sub(d, x, y); mpz_abs(d, d); mpz_gcd(d, d, *ptr_n);
    }

    if (mpz_cmp(d, *ptr_n) != 0) {
        mpz_set(*ptr_p, d);
        ret = 1;
    } else {
        ret = 0;
    }

    mpz_clears(x, y, d, NULL);
    return ret;
}

int main(int argc, char *argv[])
{
    mpz_t n, p, q;

    mpz_inits(n, p, q, NULL);
    mpz_set_str(n, argv[1], 0);

    int ret = pollard_rho(&n, &p);
    if (ret) {
        mpz_fdiv_q(q, n, p);
        gmp_printf("%Zd = %Zd * %Zd\n", n, p, q);
    } else {
        gmp_printf("%Zd is prime\n", n);
    }

    mpz_clears(n, p, q, NULL);
    return 0;
}

pollard_rho.jl

using Printf

# single-line comments

#=
multi-line comments
=#

function pollard_rho(n::BigInt)::Union{Some{BigInt}, Nothing}
    x, y, d = 2, 2, 1

    while d == 1
        x = (x*x + 1) % n
        y = (y*y + 1) % n
        y = (y*y + 1) % n
        d = gcd(abs(x-y), n)
    end

    if d != n
        return Some(d)
    else
        return nothing
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    n = parse(BigInt, ARGS[1])
    result = pollard_rho(n)
    try
        p = something(result)
        @printf "%d = %d * %d\n" n p (n÷p)
    catch e
        @printf "%d is prime\n" n
    end
end

pollard_rho.py

import sys
from math import gcd

# single-line comments

"""
multi-line comments
"""

def pollard_rho(n: int) -> int | None:
    x, y, d = 2, 2, 1

    while d == 1:
        x = (x*x + 1) % n
        y = (y*y + 1) % n
        y = (y*y + 1) % n
        d = gcd(abs(x-y), n)

    if d != n:
        return d
    else:
        return None

if __name__ == '__main__':
    n = int(sys.argv[1])
    p = pollard_rho(n)
    if p:
        print('{} = {} * {}'.format(n, p, n//p))
    else:
        print('{} is prime'.format(n))

pollard_rho_gmpy2.py

import sys
import gmpy2

# single-line comments

"""
multi-line comments
"""

def pollard_rho(n: int) -> int | None:
    n = gmpy2.mpz(n)
    x, y, d = gmpy2.mpz(2), gmpy2.mpz(2), gmpy2.mpz(1)

    while d == 1:
        x = (x*x + 1) % n
        y = (y*y + 1) % n
        y = (y*y + 1) % n
        d = gmpy2.gcd(abs(x-y), n)

    if d != n:
        return int(d)
    else:
        return None

if __name__ == '__main__':
    n = int(sys.argv[1])
    p = pollard_rho(n)
    if p:
        print('{} = {} * {}'.format(n, p, n//p))
    else:
        print('{} is prime'.format(n))

pollard_rho_opt.jl

using Printf
import Base.GMP.MPZ: add!, sub!, mul!, fdiv_r!, gcd!

# single-line comments

#=
multi-line comments
=#

# https://stackoverflow.com/a/69000702
function pollard_rho(n::BigInt)::Union{Some{BigInt}, Nothing}
    x, y, d = BigInt(2), BigInt(2), BigInt(1)

    while d == big"1"
        mul!(x, x); add!(x, big"1"); fdiv_r!(x, n)
        mul!(y, y); add!(y, big"1"); fdiv_r!(y, n)
        mul!(y, y); add!(y, big"1"); fdiv_r!(y, n)
        sub!(d, x, y); gcd!(d, abs(d), n)
    end

    if d != n
        return Some(d)
    else
        return nothing
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    n = parse(BigInt, ARGS[1])
    result = pollard_rho(n)
    try
        p = something(result)
        @printf "%d = %d * %d\n" n p (n÷p)
    catch e
        @printf "%d is prime\n" n
    end
end