#### Factoring Using Continued Fractions

This code requires the Multi-Precision Integer Arithmetic package.

```/* factor(N, K, B, F) is true if F is a factor of N (which is not a prime    */
/*   and not a perfect square), where N and F are expressed as strings, the  */
/*   small integer K is the multiplier, and the integer B is the number of   */
/*   primes to be used in the factor base.                                   */
/* Examples:                                                                 */
/*   factor(`38347921`,1,4,X)                   % 16381 * 2341               */
/*   factor(`999846002329`,1,3,X)               % 999983 * 999863            */
/*   factor(`15000000000001`,1,8,X)             % 1375951 * 10901551         */
/*   factor(`11111111111111111`,1,12,X)         % 5363222357 * 2071723       */
/*   factor(`314159265358979323`,1,11,X)        % 990371647 * 317213509      */
/*   factor(`3225913367241254663`,2,19,X)       % 1556963413 * 2071926251    */
/*   factor(`147573952589676412927`,1,18,X)     % 761838257287 * 193707721   */
/*   factor(`13191270754108226049301`,1,24,X)   % 91813 * 143675413657196977 */
factor(N0, K, B, F0):-
string_bigint(N0, N),
N \= [0],
K > 0,
multiply_short(N, K, KN),
square_root(KN, R),  % R=trunc(sqrt(KN))
multiply(R, R, R2),
subtract(KN, R2, V), % V=KN-R^2
V \= [0],
factor_base(3, B, KN, Ps0),
Ps=[2|Ps0],
reverse(Ps, [MaxP|_]),
add(R, R, Rp),       % R'=2R
multiply([5], [MaxP], B2), % B2 < MaxP^2. Put [1] to have only full relations
retractall(arg(_,_,_,_,_,_,_,_,_,_,_,_,_,_,_)),
assert(arg(Rp,Rp,[1],V,R,[1],[0],1,N,Rp,Ps,B2,[],[],[0])),
relations,
retract(arg(_,_,_,_,_,_,_,_,_,_,_,_,_,_,F)), !,
divide(N, F, _, [0]), % Sanity check
bigint_string(F, F0).

% Find the primes belonging to the factor base
factor_base(_, 1, _, []):-!.
factor_base(P, _, _, []):-
P > 32766, !,
write(`Prime in factor base exceeds 32766`), nl,
fail.
factor_base(P, B, N, [P|Qs]):-
is_prime(P),
jacobi(N, [P], J),
J >= 0, !,
P1 is P + 2,
B1 is B - 1,
factor_base(P1, B1, N, Qs).
factor_base(P, B, N, Qs):-
P1 is P + 2,
factor_base(P1, B, N, Qs).

% Find a factor by investigating relations
% This is coded using assert/retract in order to reduce heap usage
relations:-
arg(_,_,_,_,_,_,_,_,_,_,_,_,_,_,F),
\+ equal(F, [0]), !.
relations:-
retract(arg(U,Up,V,Vp,P,Pp,A,S,N,Rp,FactorBase,B2,Gs,Part,_)),
multiply(A, Up, AUp),
multiply(A, U, AU),
subtract(VpAUp, AU, V1),  % V1=A(U'-U)+V'
cast_out_primes(FactorBase, V1, Es, T),
relations_1(p(P,[S|Es],T), B2, FactorBase, N, Part, Part1, Gs, Gs1, F),
\+ equal([1], V1), !,
divide(U, V1, A1, Temp1), % A1=U div V
subtract(Rp, Temp1, U1),  % U1=R'-(U mod V)
S1 is 1 - S,
multmod(A1, P, N, Temp2),
addmod(Temp2, Pp, N, P1), % P1=(AP+P') mod N
assert(arg(U1,U,V1,V,P1,P,A1,S1,N,Rp,FactorBase,B2,Gs1,Part1,F)),
relations.
relations:-
write(`Cycling has started`), nl, % e.g. factor(`4294967297`,1,6,F)
fail.

% See if a relation has been discovered
relations_1(p(X,Es,[1]), _, FactorBase, N, Part, Part, Gs0, Gs, F):-!,
% A full relation has been discovered
solve(f(X,Es), FactorBase, N, Gs0, Gs, F).
relations_1(p(X,Es,T), B2, FactorBase, N, Part, Part, Gs0, Gs, F):-
less_than_or_equal(T, B2),
% A partial relation has been discovered
member(p(X2,Es2,T), Part),
multmod(X, X2, N, Temp),
divmod(Temp, T, N, X1), !,
% The partial relation has been combined with a previous partial relation
%   to give a full relation
solve(f(X1,Es1), FactorBase, N, Gs0, Gs, F).
relations_1(p(X,Es,T), B2, _, _, Part, [p(X,Es,T)|Part], Gs, Gs, [0]):-
less_than_or_equal(T, B2),
% A partial relation has been discovered
!.
relations_1(_, _, _, _, Part, Part, Gs, Gs, [0]).

% See if the full relation gives a factor
solve(f(X0,Es0), FactorBase, N, Gs0, Gs, Factor):-
first_odd_element(Es0, 1, K),
member(g(K,X1,Es1), Gs0), !,
% Es0 contains its first odd element at index K, and so does Es1 in Gs
multmod(X0, X1, N, X2),
solve(f(X2,Es2), FactorBase, N, Gs0, Gs, Factor).
solve(f(X0,Es0), _, _, Gs, [g(K,X0,Es0)|Gs], [0]):-
first_odd_element(Es0, 1, K),
% Es0 contains its first odd element at index K, but no element of Gs does
!.
solve(f(X,Es0), FactorBase, N, Gs, Gs, Factor):-
% Es0 contains only even elements
Es0=[_|Es], % Ignore the (even) exponent of -1
evaluate(Es, FactorBase, N, [1], Y),
X \= Y,
subtract(N, Y, Temp1),
X \= Temp1, !,
absolute_difference(X, Y, X_Y),
gcd(X_Y, N, Factor).
solve(_, _, _, Gs, Gs, [0]).
% Es0 contains only even elements, but does not give a non-trivial factor

% Find the index of the first odd element of a list
first_odd_element([X|_], I, I):-
X mod 2 =:= 1, !.
first_odd_element([_|Xs], I0, I):-
I1 is I0 + 1,
first_odd_element(Xs, I1, I).

evaluate([], [], _, Y, Y).
evaluate([E|Es], [P|FactorBase], N, Y0, Y):-
E_2 is E // 2,
powermod([P], [E_2], N, P_E_2),
multmod(Y0, P_E_2, N, Y1),
evaluate(Es, FactorBase, N, Y1, Y).

G is E + F,

/* e.g. cast_out_primes([2,3,5,7,11], 125, [0,0,3,0,0], 1) */
/*      cast_out_primes(+Ps, +T0, -Es, -T). */
cast_out_primes([], T, [], T).
cast_out_primes([P|FactorBase], T0, [E|Es], T):-
cast_out_prime(T0, P, T1, 0, E),
cast_out_primes(FactorBase, T1, Es, T).

/* e.g. cast_out_prime(125, 5, 1, 0, 3). */
/*      cast_out_prime(+T0, +P, -T, 0, -E). */
cast_out_prime(T0, P, T, E0, E):-
divide(T0, [P], T1, [0]), !,
E1 is E0 + 1,
cast_out_prime(T1, P, T, E1, E).
cast_out_prime(T, _, T, E, E).

/* is_prime(P) is true if P is a prime number.                               */
is_prime(2):-!.
is_prime(P):-
P > 2,
P mod 2 =:= 1,
SqrtP is sqrt(P),
is_prime_1(3, P, SqrtP).

is_prime_1(I, _, SqrtP):-
I > SqrtP, !.
is_prime_1(I, P, SqrtP):-
P mod I > 0,
I1 is I + 2,
is_prime_1(I1, P, SqrtP).

/* jacobi(A, N, J) is true if J is the Jacobi symbol for A and N.            */
jacobi(A, N, J):-
less_than([1], N),
jacobi_1(A, N, J).

jacobi_1([0], _, J):-!,
J=0.
jacobi_1([1], _, J):-!,
J=1.
jacobi_1(A, N, J):-
divide_short(A, 2, _, 1),
divide_short(A, 4, _, 3),
divide_short(N, 4, _, 3), !,
divide(N, A, _, N_A),
jacobi(N_A, A, J0),
J is -J0.
jacobi_1(A, N, J):-
divide_short(A, 2, _, 1), !,
divide(N, A, _, N_A),
jacobi(N_A, A, J).
jacobi_1(A, N, J):-
divide_short(N, 8, _, 1), !,
divide_short(A, 2, A_2, _),
jacobi(A_2, N, J).
jacobi_1(A, N, J):-
divide_short(N, 8, _, 7), !,
divide_short(A, 2, A_2, _),
jacobi(A_2, N, J).
jacobi_1(A, N, J):-
divide_short(A, 2, A_2, _),
jacobi(A_2, N, J0),
J is -J0.
```