## Extending Pollard’s P-1 Factorization Algorithm

### March 19, 2010

We studied John Pollard’s p-1 factorization algorithm in a previous exercise. You may recall that the algorithm finds factors of a number n by calculating the least common multiple of the integers up to some bound B, call it k, then calculates the greatest common divisor of 2k-1 and n; if the greatest common divisor is between 1 and n, it is a factor of n.

What is happening mathematically is that we are trying to find a factor p|n (that’s “p divides n“, meaning that p is a factor of n, for those not familiar with the mathematical notation), for which we know the factorization of p-1. Consider the number 15770708441 = 135979 × 115979. If we apply Pollard’s p-1 algorithm with a bound of 150, no factors are found, but if we apply Pollard’s p-1 algorithm with a bound of 180 the 135979 factor is found, because 135979 – 1 = 2 × 3 × 131 × 173; increasing the bound to include the factor 173 makes Pollard’s p-1 algorithm work. The 135979 factor is found first because 115979 – 1 = 2 × 103 × 563, and 563 is out-of-bounds.

An alternative to increasing the bound is to call a second stage that looks for a p-1 for which all factors are less than the first-stage bound except one factor that is between the first-stage bound and the second-stage bound. That is, instead of calculating lcm(1..180), we calculate lcm(1..150) × j, where j ranges from 151 to 180. For small numbers like 150 and 180, the difference doesn’t matter, but for larger numbers like B1 = 106 and B2 = 109, the difference in the computational cost is noticeable.

Your task is to write a two-stage version of Pollard’s p-1 algorithm. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisements

Pages: 1 2

### 4 Responses to “Extending Pollard’s P-1 Factorization Algorithm”

1. […] 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 […]

2. Remco Niemeijer said

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
```
3. […] 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 […]

4. Lucas A. Brown said
```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
```