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),
  add(Vp, AUp, VpAUp),
  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), !,
  vector_add(Es, Es2, Es1),
  % 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),
  vector_add(Es0, Es1, Es2),
  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).

vector_add([], [], []).
vector_add([E|Es], [F|Fs], [G|Gs]):-
  G is E + F,
  vector_add(Es, Fs, Gs).

/* 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.

LPA Index     Home Page

Valid HTML 4.01 Strict