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−1 − Vj−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 Vn − Vm−n 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 Vk − Vn 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.
Programming Praxis 
def primegen(): """ Generates primes lazily via the sieve of Eratosthenes Input: none Output: Sequence of integers Examples: >>> 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 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 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 Input: n -- integer to factor Output: Integer. A nontrivial factor of n. Example: >>> williams_pp1(315951348188966255352482641444979927) 12403590655726899403L """ 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: breakHello, what is the ilog function?
@Ruslan: The integer logarithm (ilog) of n to the base b is the largest integer e such that be ≤ n. You can see an implementation by clicking to run the code.
Hi, when I run this code Python says” name ‘count’ is not defined”, where is the mistake?