Williams’ P+1 Factorization Algorithm

June 4, 2010

Hugh Williams invented the p+1 integer factorization method in 1982 based on John Pollard’s p-1 integer factorization method; the p+1 method finds a factor p when p is smooth with respect to a bound B. We noted in a previous exercise the similar structure of Pollard’s p-1 method and Lenstra’s elliptic curve method, and Williams’ p+1 method shares the same two-stage structure; Pollard’s p-1 method performs multiplications modulo p, Lenstra’s elliptic curve method performs multiplications over an elliptic curve, and Williams’ p+1 method performs multiplications over a quadratic field using Lucas sequences (similar to the Lucas chain in the Baillie-Wagstaff primality-testing exercise). A good description of Williams’ p+1 method is given at http://www.mersennewiki.org/index.php/P_Plus_1_Factorization_Method.

Williams’ p+1 method uses the Lucas sequence V0 = 2, V1 = A, Vj = A · Vj−1Vj−2, with all operations modulo N, where A is an integer greater than 2. Multiplication by successive prime powers is done as in Pollard’s p-1 algorithm; when all the products to the bound B are accumulated, the greatest common divisor of N and the result VB − 2 may be a divisor of N if all the factors of p+1 are less than B. However, if A2 − 4 is a quadratic non-residue of p, the calculation will fail, so several attempts must be made with different values of A before concluding that p+1 is not B-smooth.

Given the addition formula Vm+n = Vm VnVmn and the doubling formula V2m = Vm × Vm − 2, multiplication is done by means of a Lucas sequence that computes two values at a time. Starting with the pair Vkn and V(k+1)n, and looking at the bits in the binary representation of the multiplier M, excluding the most significant bit (which is always 1) and working from most significant to least significant, calculate the pair V2kn, V(2k+1)n when the bit is zero and the pair V(2k+1)n, V2(k+1)n when the bit is one.

The second stage finds factors that are B1-smooth except for a single prime factor in the range B1 to B2. It can be done in a similar manner to the second stage of Pollard’s p-1 method, but a little bit of algebra gives us a better second stage. Assuming that Vn is the point that survives the first stage, multiply the products VkVn for each k divisible by 6 between B1 and B2, then take the greatest common divisor of the product with n to reveal the factor.

Your task is to write a program that factors integers using Williams’ 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.


Pages: 1 2

4 Responses to “Williams’ P+1 Factorization Algorithm”

  1. Lucas A. Brown said
    def primegen():
        Generates primes lazily via the sieve of Eratosthenes
        Input: none
            Sequence of integers
        >>> list(takewhile(lambda x: x < 100, primegen()))
        [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
        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
                    next, step = q + 2*p, 2*p
                    while next in sieve: next += step
                    sieve[next] = step
                    p = ps.next()
                    q = p**2
                step = sieve.pop(n)
                next = n + step
                while next in sieve: next += step
                sieve[next] = step
            n += 2
    def mlucas(v, a, n):
        """ Helper function for williams_pp1().  Multiplies along a Lucas sequence modulo n. """
        v1, v2 = v, (v**2 - 2) % n
        for bit in bin(a)[3:]: v1, v2 = ((v1**2 - 2) % n, (v1*v2 - v) % n) if bit == "0" else ((v1*v2 - v) % n, (v2**2 - 2) % n)
        return v1
    def williams_pp1(n):
        Williams' p+1 integer factoring algorithm
            n -- integer to factor
            Integer.  A nontrivial factor of n.
        >>> williams_pp1(315951348188966255352482641444979927)
        for v in count(1):
            for p in primegen():
                e = ilog(sqrt(n), p)
                if e == 0: break
                for _ in xrange(e): v = mlucas(v, p, n)
                g = gcd(v - 2, n)
                if 1 < g < n: return g
                if g == n: break
  2. Ruslan Khamidullin said

    Hello, what is the ilog function?

  3. programmingpraxis said

    @Ruslan: The integer logarithm (ilog) of n to the base b is the largest integer e such that ben. You can see an implementation by clicking to run the code.

  4. Pierre said

    Hi, when I run this code Python says” name ‘count’ is not defined”, where is the mistake?

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 )

Connecting to %s

%d bloggers like this: