#### 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).
```