N is the integer to be factored, which must be neither a prime number nor a perfect square.
P[0]=floor(sqrt(N)) Q[0]=1 Q[1]=N-P[0]^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[0]=sqrt(Q[i]) b=floor((floor(sqrt[N])+P[i])/Q[0]) P[0]=bQ[0]-P[i] Q[1]=(N-P[0]^2)/Q[0] 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, [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).