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.

About these ads

Pages: 1 2

5 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. Nice Blog… Going to go through one post a time and read this whole website like it was designed for a magazine.

  4. [...] 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 [...]

  5. 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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 627 other followers

%d bloggers like this: