Factoring (Fermat's Method)

This code requires the Multi-Precision Integer Arithmetic package.

/* fermat(N, P, Q) is true if P is the largest factor of N less than or      */
/*   equal to sqrt(N), and Q = N / P, where N, P and Q are strings           */
/*   containing integer values, and the value of N is odd and greater than 1.*/
/*   This is coded using assert/retract in order to reduce heap usage.       */
/* e.g. fermat(`8616460799`, P, Q).                                          */
/*      fermat(`999846002329`, P, Q).                                        */
fermat(N0, P0, Q0):-
  string_bigint(N0, N),
  less_than([1], N),           % n > 1
  divide_short(N, 2, _, 1),    % n is odd
  square_root(N, S),           % s=floor(sqrt(n))
  add(S, S, TwoS),
  add(TwoS, [1], X),           % x=2*s+1
  multiply(S, S, S2),
  subtract(N, S2, R),          % r=n-s^2
  retractall(f(_,_,_,_,_)),
  assert(f(R,X,[1],[1],[1])),
  fermat_1,
  retract(f(_,_,_,P,Q)), !,
  bigint_string(P, P0),
  bigint_string(Q, Q0).

fermat_1:-
  f([0],_,_,_,_),
  retract(f(_,X,Y,_,_)), !,    % Factor found
  subtract(X, Y, Diff),
  divide_short(Diff, 2, P, _), % p=(x-y)/2
  add(X, Y, Sum),
  subtract(Sum, [2], Sum1),
  divide_short(Sum1, 2, Q, _), % q=(x+y-2)/2
  assert(f([0],X,Y,P,Q)).
fermat_1:-                     % Factor not yet found
  retract(f(R,X,Y,P,Q)), !,
  add(X, [2], X1),             
  subtract(X, R, X_R),         
  fermat_2(X_R, Y, Y1, R1),
  assert(f(R1,X1,Y1,P,Q)),
  fermat_1.

fermat_2(X_R0, Y0, Y, R):-
  less_than_or_equal(X_R0, Y0), !,
  add(Y0, [2], Y),          
  subtract(Y0, X_R0, R).    
fermat_2(X_R0, Y0, Y, R):-
  subtract(X_R0, Y0, X_R1),  
  add(Y0, [2], Y1),         
  fermat_2(X_R1, Y1, Y, R).

LPA Index     Home Page