This code requires the Multi-Precision Integer Arithmetic package.
% rho(Integer, Factor) is true if Factor is a factor of the Integer,
% where Integer and Factor are expressed as strings. If Factor = `1`,
% no factor has been found.
% e.g. rho(`11337773311`, `151451`).
% rho(`999846002329`, X).
% rho(`1824032775457`, X).
% rho(`576460752303423487`, X).
rho(N0, G0):-
string_bigint(N0, N),
rho_1(N, G0).
rho_1(N, `2`):-
divide(N, [2], _, Mod),
Mod = [0], !.
rho_1(N, G0):-
string_bigint( `3`, A), % Or any other integer less than N
string_bigint(`17`, X), % Or any other integer less than N
rho_2(N, A, X, X, G),
bigint_string(G, G0).
rho_2(N, A, X, Y, G):-
% X1 is (X * X + A) mod N
multiply(X, X, X_2),
add(X_2, A, X_2_A),
divide(X_2_A, N, _, X1),
% T is (Y * Y + A) mod N
multiply(Y, Y, Y_2),
add(Y_2, A, Y_2_A),
divide(Y_2_A, N, _, T),
% Y1 is (T * T + A) mod N
multiply(T, T, T_2),
add(T_2, A, T_2_A),
divide(T_2_A, N, _, Y1),
% G0 is gcd(abs(X1 - Y1), N)
absolute_difference(X1, Y1, X_Y),
X_Y \= [0], !,
gcd(X_Y, N, G0),
(less_than([1], G0) ->
G = G0
;
rho_2(N, A, X1, Y1, G)
).
rho_2(_, _, _, _, [1]).