#### Shanks' SQUFOF Factoring Algorithm

N is the integer to be factored, which must be neither a prime number nor a perfect square.

```P=floor(sqrt(N))
Q=1
Q=N-P^2
Repeat
b=floor(floor(sqrt(N)+P[i-1])/Q[i])
P[i]=bQ[i]-P[i-1]
Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
until Q[i] is a perfect square
Q=sqrt(Q[i])
b=floor((floor(sqrt[N])+P[i])/Q)
P=bQ-P[i]
Q=(N-P^2)/Q
Repeat
b=floor(floor(sqrt(N)+P[i-1])/Q[i])
P[i]=bQ[i]-P[i-1]
Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
until P[i+1]=P[i]
```

Then gcd(N,P[i]) is a non-trivial factor of N.

The following code requires the Multi-Precision Integer Arithmetic and Miller-Rabin Primality Test packages.

With LPA Win-Prolog, this will factor integers with up to 30 digits or so, often remarkably quickly. With other Prolog implementations, the limit may be 2^60-1.

```% prime_factors(N, Fs) is true if the elements of the list Fs are the prime
%   factors of N in ascending order of sequence, where N and the elements
%   of Fs are all strings.
% e.g. prime_factors(`119999999999999999999999911`, Fs)
%   gives Fs = [`451356383`,`265865299616245817`]
prime_factors(N0, Fs0):-
string_bigint(N0, N),
prime_factors_1(N, [], Fs),
bigints_strings(Fs, Fs0).

% Here N and the elements of Fs are bigints
prime_factors_1(N, Fs0, Fs):-
prime_factors_2(1, N, F),     % F is a factor of N
\+ equal(N, F), !,            % N is composite
divide(N, F, G, ),         % N = F * G
prime_factors_1(F, Fs0, Fs1), % Insert the factors of F
prime_factors_1(G, Fs1, Fs).  % Insert the factors of G
prime_factors_1(N, Fs0, Fs):-   % N is prime
insert(Fs0, N, Fs).           % Insert N

% Finds a multiplier K that enables N to be factored
prime_factors_2(K, N, F):-
squfof_1(N, K, F), !.         % F is a factor of N
prime_factors_2(K, N, F):-
K1 is K + 1,                  % Try the next value of K
prime_factors_2(K1, N, F).

% insert(Ys, X, Zs) is true if Zs is the result of inserting the bigint X
%   into the ordered list of bigints Ys.
insert([Y|Ys], X, [Y|Zs]):-less_than(Y, X), !, insert(Ys, X, Zs).
insert(Ys, X, [X|Ys]).

% squfof(N, F) is true if F is a factor, not necessarily prime, of N, where
%   N and F are both strings.
% e.g. squfof(`11111111111111111`, F) gives F=`2071723`
%      squfof(`11111118881111111`, F) gives F=`11111118881111111` (a prime)
%      squfof(`314159265358979323`, F) gives F=`317213509`
%      squfof(`4016178294651270164667643`, F) gives F=`596396179789`
%      squfof(`119999999999999999999999911`, F) gives F=`451356383`
squfof(N0, F0):-
string_bigint(N0, N),
for(1, 1000, K),                 % Find a suitable multiplier
squfof_1(N, K, F), !,
bigint_string(F, F0).

squfof_1(N, _, N):-
T=5,                             % Or perhaps a slightly bigger number
miller_rabin(N, T), !.           % N is almost surely prime
squfof_1(N, _, F):-
square_root(N, F),
multiply(F, F, N), !.            % N is a perfect square
squfof_1(N, K, F):-
multiply_short(N, K, KN),
square_root(KN, P0b),
bigint_int(P0b, P0),             % P0 = floor(sqrt(KN)); Q0 = 1
multiply(P0b, P0b, T1),
subtract(KN, T1, Q1b),
bigint_int(Q1b, Q1),             % Q1 = KN - P0 * P0
reduce_1(P0, P0, 1, Q1, Pi, S0), % S0 = sqrt(Qi)
R0 is ((P0 + Pi) // S0) * S0 - Pi,
int_bigint(R0, R0b),
multiply(R0b, R0b, T2),
subtract(KN, T2, T3),
divide_short(T3, S0, S1b, 0),
bigint_int(S1b, S1),             % S1 = (KN - R0 * R0) div S0
reduce_3(P0, R0, S0, S1, Pj),
int_bigint(Pj, Pjb),
gcd(N, Pjb, F),
\+ equal(, F),
\+ equal(F, N),                  % F is a non-trivial factor of N
divide(KN, F, _, ).           % Sanity check

reduce_1(X, P0, Q0, Q1, P, Q):-
Q1 > 0,                          % Avoids division by 0 if KN is a square
B1 is (X + P0) // Q1,
P1 is B1 * Q1 - P0,
P1 =\= P0,
Q2 is Q0 + B1 * (P0 - P1),
B2 is (X + P1) // Q2,
P2 is B2 * Q2 - P1,
Q3 is Q1 + B2 * (P1 - P2),
reduce_2(X, P1, P2, Q2, Q3, P, Q).

reduce_2(_, _, P, Q2, _, P, Q):-
Q is sqrt(Q2),
Q =:= int(Q), !.                 % Q2 = Q * Q
reduce_2(X, P1, P2, Q2, Q3, P, Q):-
P2 =\= P1,
reduce_1(X, P2, Q2, Q3, P, Q).

reduce_3(X, P0, Q0, Q1, P):-
B1 is (X + P0) // Q1,
P1 is B1 * Q1 - P0,
P1 =\= P0,
Q2 is Q0 + B1 * (P0 - P1), !,
reduce_3(X, P1, Q1, Q2, P).
reduce_3(_, P, _, _, P).

for(I, J, I):-I =< J.
for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K).
```