Pollard Rho and Brent-Pollard Rho factoring algorithm

brent_pollard_rho.py

# brent_pollard_rho.py

import sys
import gmpy2
from random import randint
from math import gcd

def miller_rabin(n, k=20):
    s, d = 0, n-1

    while d % 2 == 0:
        s += 1
        d //= 2

    for i in range(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1:
            continue
        for r in range(s):
            if x == n-1:
                break
            x = (x*x) % n
        else:
            return False

    return True

# https://xn--2-umb.com/09/12/brent-pollard-rho-factorisation/
# https://qiita.com/Kiri8128/items/eca965fe86ea5f4cbb98
def brent_pollard_rho(n):
    m, is_exact = gmpy2.iroot(n, 8)
    y, r, q, d = 2, 1, 1, 1
    while d == 1:
        x = y
        for i in range(r):
            y = (y*y + 1) % n
        for k in range(0, r, m):
            ys = y
            for i in range(min(m, r - k)):
                y = (y*y + 1) % n
                q = (q * abs(x - y)) % n
            d = gcd(q, n)
            if d > 1:
                break
        r <<= 1
    if d == n:
        while True:
            ys = (ys*ys + 1) % n
            d = gcd(abs(x - ys), n)
            if d > 1:
                break
    return d

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

gistfile1.txt

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

real    0m16.533s
user    0m16.520s
sys     0m0.008s

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

real    0m10.135s
user    0m10.119s
sys     0m0.008s

pollard_rho.py

# pollard_rho.py

import sys
from random import randint
from math import gcd

def miller_rabin(n, k=20):
    s, d = 0, n-1

    while d % 2 == 0:
        s += 1
        d //= 2

    for i in range(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1:
            continue
        for r in range(s):
            if x == n-1:
                break
            x = (x*x) % n
        else:
            return False

    return True

def pollard_rho(n):
    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

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