Shanks' SQUFOF Factoring Algorithm

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

until Q[i] is a perfect square
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, [0]),         % 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([1], F),
  \+ equal(F, N),                  % F is a non-trivial factor of N
  divide(KN, F, _, [0]).           % 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).

LPA Index     Home Page

Valid HTML 4.01 Strict