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).