Miller-Rabin Primality Test

This code requires the Multi-Precision Integer Arithmetic package.

% miller_rabin(N, T) succeeds if N (a bigint) is "almost surely prime" and  
%   fails if N is "definitely not prime".  This is because if the test says 
%   T (an integer) times in a row that N is "probably prime", we can say    
%   that N is "almost surely prime" with a probability of less than (1/4)^T 
%   of being wrong.                                                         
miller_rabin([2], _):-!.
miller_rabin(N, T):-
  less_than([2], N),
  divide_short(N, 2, _, 1),
  T > 0,
  subtract(N, [1], N_1),
  miller_rabin_2(N_1, D, 0, S),  % N - 1 = 2^S * D
  miller_rabin_1(T, N, S, D).

% Repeat T times...
miller_rabin_1(0, _, _, _):-!.
miller_rabin_1(I, N, S, D):-
  repeat,
    random_bigint(N, B),
  less_than([1], B), !,
  powermod(B, D, N, X),          % X=B^D mod N
  miller_rabin_3(0, S, N, X),
  I1 is I - 1,
  miller_rabin_1(I1, N, S, D).

% miller_rabin_2(N, D, 0, S) is true if N = 2^S*D, and D is odd.            
miller_rabin_2(N, N, S, S):-
  divide_short(N, 2, _, 1), !.   % N mod 2 = 1
miller_rabin_2(N, D, S0, S):-    % N mod 2 = 0
  S1 is S0 + 1,
  divide_short(N, 2, N1, 0),     % N1 = N div 2
  miller_rabin_2(N1, D, S1, S).

% Succeeds if N is probably prime.
% For each R in the range [0, S-1]...
miller_rabin_3(0, _, _, [1]):-!. % B^D = 1(mod N)
miller_rabin_3(_, _, N, X):-  
  subtract(N, [1], X), !.        % (B^D)^(2*R) = -1(mod N)
miller_rabin_3(R, S, N, X):-
  R1 is R + 1,
  R1 < S,
  multmod(X, X, N, X1),
  miller_rabin_3(R1, S, N, X1).

LPA Index     Home Page

Valid HTML 4.0!