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]).