Python Pollard P-1 factorization

Question

I'm trying to implement Pollard's P-1 factorization in Python. Note that the Rho method has some answers but this p-1 is different and the best I can give you here about p-1 is the wiki and the Wolfram:

http://en.wikipedia.org/wiki/Pollard's_p_%E2%88%92_1_algorithm

http://mathworld.wolfram.com/Pollardp-1FactorizationMethod.html

This is factoring something from n but consistently does not find p. The np and sp are from numpy and scipy respectively. So the built-in functions for sp.uint64 is an unsigned long 64 int (because of the size of the expected integers) and np.prod(p) is the cumulative product pi of the list p:

def pm1_attack(n,b):
  p = [2,3,5,7,11,13,17]; i=19; a=2
  while i<b:
    if is_prime(i,10): p.append(i)
    i+=2;
  k = sp.uint64(np.prod(p)); q = power2(a,k,n)
  g = euc_al_i((q-1),n)
  print "product pi: ",k
  print "q: ",q
  print "g: ",g
  #return a

print "pollard_pm1_attack(n,b): ",pollard_pm1_attack(n,2000)

Output does not find p:

Python 2.7 (r27:82525, Jul  4 2010, 09:01:59) [MSC v.1500 32 bit (Intel)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>> 
p = 1300199 
q = 2063507 

euler_totient = 2682966374188 

common n = 2682969737893 


public key e = 1588051820871 

secret key d = 2410616084843 

cleartext message = test 

encoded message = 1489995542681 

decoded message = test 

check_rsa = Successful encryption and decrytion. Message is authenticated. 

pollard_pm1_attack(n,b):  product pi:  18446744073460481730
q:  2391570546599
g:  1
None
>>> 

I'm learning Python so it may be some simple mistake. The power2() function uses exponentiation by squaring and is basically a super-charged pow() for very large integers. The euc_al_i() is just the gcd. You can use whatever gcd() you like, but since I'm learning I wanted to make these myself.

I'm trying to find out what went so horribly wrong here that it doesn't find p from even relatively small n (as small as 20 bit-length).

Answer 1

I don't know what np.prod and sp.uint64 do, but I can tell you about the p - 1 algorithm, which was invented by John Pollard in 1974.

Pollard's algorithm is based on Fermat's Little Theorem a ^ p == a (mod p), which when a != 0 can be stated a ^ (p - 1) == 1 (mod p) by dividing a through the expression. As a consequence, if p - 1 divides m, then p divides gcd(2^m - 1, n). Pollard's p - 1 algorithm computes m as the least common multiple of the integers less than a bound b, so that if all the factors of p - 1 are less than b, then gcd(2 ^ lcm(1..b) - 1, n) is a factor of n. The least common multiple is computed by multiplying primes less than b by their multiplicities less than b:

function pminus1(n, b)
    c := 2
    for p in primes(b)
        pp := p
        while pp < b
            c := powerMod(c, p, n)
            pp := pp * p
    g := gcd(c-1, n)
    if 1 < g < n return g
    error "factorization failed"

An optional second stage searches for a "large prime" between b1 and b2 that combines with the least common multiple of the first stage to find a factor. The second stage requires only a modular multiplication for each prime rather than a modular exponentiation, making it quite fast, and the second-stage gcds can be computed in batches. The cache is small but important to the efficiency of the function.

function pminus1(n, b1, b2)
    c := 2
    for p in primes(b1)
        pp := p
        while pp < b
            c := powerMod(c, p, n)
            pp := pp * p
    g := gcd(c-1, n)
    if 1 < g < n return g
    k := 0
    for q in primes(b1, b2)
        d := q - p
        if d is not in cache
            x := powerMod(c, d, n)
            store d, x in cache
        c := (c * x(d)) % n
        p := q
        k := k + 1
        if k % 100 == 0
            g := gcd(c-1, n)
            if 1 < g < n return g
    g := gcd(c-1, n)
    if 1 < g < n return g
    error "factorization failed"

It is possible that Pollard's p - 1 method may fail to find a factor of n; it depends on the factorization of n - 1 and the bounds you have chosen. The way to check is to factor n - 1 yourself, then call Pollard's method with a b that is larger than the largest factor of n - 1. For instance, if you want to factor n = 87463 = 149 * 587, note that n - 1 = 87462 = 2 * 3 * 3 * 43 * 113, so call the one-stage algorithm with b = 120 or the two-stage algorithm with b1 = 50 and b2 = 120 and see if you find a factor.

I discuss Pollard's p - 1 factorization algorithm, along with several other factorization algorithms, at my blog. I also give there implementations of the powerMod and gcd functions, in case you are having trouble with your implementations of those functions. And I wrote a simple implementation of the single-stage algorithm, in Python, at http://ideone.com/wdyjxK.

Answer 2

Here is a two stage version implemented in python.

from math import gcd, log

def pminus1(n, B1, B2):
    log_B1 = log(B1)
    M = 1
    primes = prime_sieve()
    for p in primes:
        if p > B1:
            break
        M *= p**int(log_B1/log(p))
    M = pow(2, M, n)
    g = gcd(M-1, n)
    if 1 < g < n:
        return True
    if g == n:
        return False

    # Start part 2.
    cache = {0:M}
    S = (M*M) % n
    for d in range(2, int(log(B2)**2), 2):
        cache[d] = cache[d-2] * S) % n

    HQ = M
    for k, q in enumerate(primes):
        if q > B2:
            break
        d = q - p
        HQ = (HQ * cache[d]) % n
        M = (M * (HQ-1)) % n
        p = q
        if k % 200 == 0:
            if 1 < gcd(M, n) < n:
                return True
    return 1 < gcd(M, n) < n