Extending Pollard’s P-1 Factorization Algorithm
March 19, 2010
We begin by restating the one-stage algorithm, in a somewhat different form than it appeared in the previous exercise:
(define (pollard1 n b1)
(let stage1 ((a 2) (i 2))
(if (< i b1)
(stage1 (expm a i n) (+ i 1))
(let ((d (gcd (- a 1) n)))
(if (< 1 d n) d #f)))))
The two-stage algorithm starts the same, but the continuation loops over the integers from B1 to B2, making the same computation as the first stage:
(define (pollard2 n b1 b2)
(let stage1 ((a 2) (i 2))
(if (< i b1)
(stage1 (expm a i n) (+ i 1))
(let ((d (gcd (- a 1) n)))
(if (< 1 d n) (list 'stage1 d)
(let stage2 ((j b1))
(if (= j b2) #f
(let ((d (gcd (- (expm a j n) 1) n)))
(if (< 1 d n) (list 'stage2 d)
(stage2 (+ j 1)))))))))))
We can see how this works using the sample problem given above:
> (pollard1 15770708441 150)
#f
> (pollard1 15770708441 180)
135979
> (pollard2 15770708441 150 180)
(stage2 135979)
We use expm from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/VpaXvPEN.
Programming Praxis 
[…] Praxis – Extending Pollard’s P-1 Factorization Algorithm By Remco Niemeijer In today’s Programming Praxis exercise we need to write an improved version of a factorization algorithm. I […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/03/19/programming-praxis-extending-pollard%E2%80%99s-p-1-factorization-algorithm/ for a version with comments):
import Data.Bits import Data.List expm :: Integer -> Integer -> Integer -> Integer expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 . filter (flip testBit 0 . snd) . zip (iterate (flip mod m . (^ 2)) b) . takeWhile (> 0) $ iterate (`shiftR` 1) e pollard :: (Integer -> t) -> (Integer -> t) -> Integer -> Integer -> t pollard found notFound n b1 = f 2 2 where f a i | i < b1 = f (expm a i n) (i + 1) | 1 < d && d < n = found d | otherwise = notFound a where d = gcd (a - 1) n pollard1 :: Integer -> Integer -> Maybe Integer pollard1 = pollard Just (const Nothing) pollard2 :: Integer -> Integer -> Integer -> Maybe (String, Integer) pollard2 n b1 b2 = pollard (Just . ((,) "stage1")) (f b1) n b1 where f j a | j == b2 = Nothing | 1 < d && d < n = Just ("stage2", d) | otherwise = f (j + 1) a where d = gcd (expm a j n - 1) n[…] have studied John Pollard’s p−1 algorithm for integer factorization on two previous occasions, giving first the basic single-stage algorithm and later adding a second stage. In today’s […]
def primegen(): """ Generates primes lazily via the sieve of Eratosthenes. Code shamelessly stolen from someone else's work. No idea where I got it. Input: none Output: Sequence of integers """ yield 2; yield 3; yield 5; yield 7; yield 11; yield 13 ps = primegen() # yay recursion p = ps.next() and ps.next() q, sieve, n = p**2, {}, 13 while True: if n not in sieve: if n < q: yield n else: next, step = q + 2*p, 2*p while next in sieve: next += step sieve[next] = step p = ps.next() q = p**2 else: step = sieve.pop(n) next = n + step while next in sieve: next += step sieve[next] = step n += 2 def ispower(n): """ If n is a perfect power, return the largest integer (in terms of absolute value) that, when squared/cubed/etc, yields n. If n is not a perfect power, return 0. Input: n -- an integer Output: An integer Examples: >>> [ispower(n) for n in [64, 25, -729, 1729]] [8, 5, -9, 0] """ for p in primegen(): r = introot(n, p) if r is None: continue if r ** p == n: return r if r == 1: return 0 def pollard_pm1(n, B1=100, B2=1000): # TODO: What are the best default bounds and way to increment them? """ Integer factoring function. Uses Pollard's p-1 algorithm. Note that this is only efficient iff the number to be factored has a prime factor p such that p-1's largest prime factor is "small". In this implementation, that tends to mean less than 10,000,000 or so. Input: n -- number to factor B1 -- Natural number. Bound for phase 1. Default == 100. B2 -- Natural number. Bound for phase 2. Must be > B1. Default == 100. Output: A factor of n. Examples: >>> pollard_pm1(1275683450258216989546041841) 1224040923709997L """ if isprime(n): return n m = ispower(n) if m: return m while True: pg = primegen() q = 2 # TODO: what about other initial values of q? p = pg.next() while p <= B1: q, p = pow(q, p**ilog(B1, p), n), pg.next() g = gcd(q-1, n) if 1 < g < n: return g while p <= B2: q, p = pow(q, p, n), pg.next() g = gcd(q-1, n) if 1 < g < n: return g # These bounds failed. Increase and try again. B1 *= 10 B2 *= 10