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