Factoring (Williams' p+1 Method)

We want to find a prime factor p of N when p+1 is smooth.

Choose some integer A>2 and define the sequences:

  U[0]=0, U[1]=1, U[j]=AU[j-1]-U[j-2]
  V[0]=2, V[1]=A, V[j]=AV[j-1]-V[j-2]

where all operations are performed modulo N.

Then any odd prime p divides both gcd(N,U[M]) and gcd(N,V[M]-2) whenever M is a multiple of p-(D/p), where D=A^2-4 and (D/p) is the Jacobi symbol. We require that (D/p) be -1, that is, D should be a quadratic non-residue modulo p. But as we don't know p beforehand, more than one value of A may be required before finding a solution. If (D/p)=1, this algorithm degenerates into a slow version of Pollard's p-1 algorithm. The values of M used are successive factorials.

The elements of the Lucas sequence (A,B) are defined by:

  U[0]=0, U[1]=1, U[j]=AU[j-1]-BU[j-2]
  V[0]=2, V[1]=A, V[j]=AV[j-1]-BV[j-2]

and it can be shown that:

  U[2j]=U[j]V[j]
  U[2j+1]=U[j+1]V[j]-B^j
  V[2j]=V[j]^2-2B^j
  V[2j+1]=V[j+1]V[j]-AB^j

This implies that U[i] and V[i] can be calculated in O(log i) steps.

For the P+1 algorithm, B=1 and we don't use the elements U[j], and all operations are performed modulo N. So:

  V[2j]=V[j]^2-2
  V[2j+1]=V[j+1]V[j]-A

and the calculation of V[M] becomes:

  x=A             /* V[j]=V[1]   */
  y=(A^2-2)%N     /* V[j+1]=V[2] */
  for each bit of M to the right of the most significant bit
    if the bit is 1
      x=(x*y-a)%N /* V[2j+1] */
      y=(y^2-2)%N /* V[2j+2] */
    else
      y=(x*y-a)%N /* V[2j+1] */
      x=(x^2-2)%N /* V[2j]   */
  V[M]=x

This is similar to "left-to-right" exponentiation:

  V^(2j)=(V^j)^2
  V^(2j+1)=V(V^j)^2

where the calculation of V^M is:

  result=V
  for each bit of M to the right of the most significant bit
    result=result*result
    if the bit is 1
      result=result*V

The following code requires the Multi-Precision Integer Arithmetic package.

/* pplus1(N, F) is true if F is a factor of N, where N and F are both        */
/*   strings. If no factor is found after 10000 iterations, `1` is returned. */
/*   The initial value of A is 4.  F is not guaranteed to be prime -- for    */
/*   example, pplus1(`207`, F) gives F=`9`.                                  */
/* e.g. pplus1(`27198662590716548097867889`, F).                             */
pplus1(N, F):-pplus1(N, 4, 10000, F).

/* pplus1(N, A, MaxI, F) is true if F is a factor of N, where N and F are    */
/*   both strings, and A and MaxI are integers. If no factor is found after  */
/*   MaxI iterations, `1` is returned. It may be necessary to try different  */
/*   values of A, as the algorithm requires that A^2-4 be a quadratic        */
/*   non-residue modulo F.  F is not guaranteed to be prime -- for example,  */
/*   pplus1(`207`, 4, 100, F) gives F=`9`.                                   */
/*   This is coded using assert/retract in order to reduce heap usage.       */
/* e.g. pplus1(`25443025601020650513668093`, 3, 1000, F).                    */
pplus1(`0`, _, _, `1`):-!.
pplus1(`1`, _, _, `1`):-!.
pplus1(N0, A, MaxI, F0):-
  string_bigint(N0, N),
  retractall(q(_,_,_,_,_)),
  int_bigint(A, BigA),
  assert(q(1,MaxI,N,BigA,[1])),
  pplus1_1,
  retract(q(_,_,_,_,F)), !,
  bigint_string(F, F0).

pplus1_1:-                         % No factor found after MaxM iterations
  q(M,MaxM,_,_,_),
  M=MaxM, !.
pplus1_1:-
  q(M,MaxM,N,V,_),
  (M - 1) mod 10 =:= 0,            % Calculate the gcd every 10 iterations
  subtract(V, [2], V_2),
  gcd(V_2, N, F),
  less_than([1], F),
  less_than(F, N),                 % Factor found
  retract(q(_,_,_,_,_)), !,
  assert(q(M,MaxM,N,V,F)).
pplus1_1:-                         % Factor not yet found
  retract(q(M,MaxM,N,V,F)), !,
  bitcount(M, K),
  multmod(V, V, N, V2),
  subtmod(V2, [2], N, Y),
  K1 is K - 1,
  pplus1_2(K1, M, N, V, V, Y, V1), % x=a and y=(a^2-2)%n
  M1 is M + 1,
  assert(q(M1,MaxM,N,V1,F)),
  pplus1_1.

% For each bit of the integer M starting with the msb...
pplus1_2(0, _, _, _, X, _, X):-!.
pplus1_2(K, M, N, A, X, Y, V):-
  is_one(M, K), !,
  multmod(X, Y, N, XY),
  subtmod(XY, A, N, X1),           % x=(x*y-a)%n
  multmod(Y, Y, N, YY),  
  subtmod(YY, [2], N, Y1),         % y=(y^2-2)%n
  K1 is K - 1,
  pplus1_2(K1, M, N, A, X1, Y1, V).
pplus1_2(K, M, N, A, X, Y, V):-
  multmod(X, Y, N, XY),
  subtmod(XY, A, N, Y1),           % y=(x*y-a)%n
  multmod(X, X, N, XX),
  subtmod(XX, [2], N, X1),         % x=(x^2-2)%n
  K1 is K - 1,
  pplus1_2(K1, M, N, A, X1, Y1, V).
  
/* bitcount(M, Bits) is true if Bits is the number of significant binary     */
/*   digits in the positive integer M.                                       */
/* e.g. bitcount(19, 5).                                                     */
bitcount(M, Bits):-
  M >= 0,
  bitcount_1(M, 0, Bits).

bitcount_1(0, Bits, Bits):-!.
bitcount_1(M, Bits0, Bits):-
  Bits1 is Bits0 + 1,
  M1 is M // 2,
  bitcount_1(M1, Bits1, Bits).
  
/* is_one(M, K) is true if the K-th bit of the positive integer M is 1,      */
/*   where the least significant bit is bit 1.                               */
is_one(M, 1):-!,
  1 =:= M mod 2.
is_one(M, K):-
  M1 is M // 2,
  K1 is K - 1,
  is_one(M1, K1).

LPA Index     Home Page

Valid HTML 4.0 Strict