Miller–Rabin primality test: Difference between revisions
Content added Content deleted
(Optimisation of existing solution) |
|||
Line 1,417: | Line 1,417: | ||
to permit use of integers of arbitrary precision. |
to permit use of integers of arbitrary precision. |
||
<lang erlang> |
<lang erlang> |
||
-module(miller_rabin). |
|||
-export([is_prime/1, power/2]). |
-export([is_prime/1, power/2]). |
||
Line 1,507: | Line 1,508: | ||
power(B, E, Acc) -> |
power(B, E, Acc) -> |
||
power(B, E - 1, B * Acc).</lang> |
power(B, E - 1, B * Acc).</lang> |
||
The above code optimised as follows: |
|||
- more efficient exponentiation operation |
|||
- parallel execution of tests |
|||
The parallel executions reduced the time to run the test from |
|||
53s to 11s on a quad-core 17 with 16 GB ram. The performance |
|||
gain from the improved exponentiation was not evaluated. |
|||
<lang erlang> |
|||
%%% @author Tony Wallace <tony@resurrection> |
|||
%%% @copyright (C) 2021, Tony Wallace |
|||
%%% @doc |
|||
%%% For details of the algorithms used see |
|||
%%% https://en.wikipedia.org/wiki/Modular_exponentiation |
|||
%%% @end |
|||
%%% Created : 21 Jul 2021 by Tony Wallace <tony@resurrection> |
|||
-module mod. |
|||
-export [mod_mult/3,mod_exp/3,binary_exp/2,test/0]. |
|||
mod_mult(I1,I2,Mod) when |
|||
I1 > Mod, |
|||
is_integer(I1), is_integer(I2), is_integer(Mod) -> |
|||
mod_mult(I1 rem Mod,I2,Mod); |
|||
mod_mult(I1,I2,Mod) when |
|||
I2 > Mod, |
|||
is_integer(I1), is_integer(I2), is_integer(Mod) -> |
|||
mod_mult(I1,I2 rem Mod,Mod); |
|||
mod_mult(I1,I2,Mod) when |
|||
is_integer(I1), is_integer(I2), is_integer(Mod) -> |
|||
(I1 * I2) rem Mod. |
|||
mod_exp(Base,Exp,Mod) when |
|||
is_integer(Base), |
|||
is_integer(Exp), |
|||
is_integer(Mod), |
|||
Base > 0, |
|||
Exp > 0, |
|||
Mod > 0 -> |
|||
binary_exp_mod(Base,Exp,Mod); |
|||
mod_exp(_,0,_) -> 1. |
|||
binary_exp(Base,Exponent) when |
|||
is_integer(Base), |
|||
is_integer(Exponent), |
|||
Base > 0, |
|||
Exponent > 0 -> |
|||
binary_exp(Base,Exponent,1); |
|||
binary_exp(_,0) -> |
|||
1. |
|||
binary_exp(_,0,Result) -> |
|||
Result; |
|||
binary_exp(Base,Exponent,Acc) -> |
|||
binary_exp(Base*Base,Exponent bsr 1,Acc * exp_factor(Base,Exponent)). |
|||
binary_exp_mod(Base,Exponent,Mod) -> |
|||
binary_exp_mod(Base rem Mod,Exponent,Mod,1). |
|||
binary_exp_mod(_,0,_,Result) -> |
|||
Result; |
|||
binary_exp_mod(Base,Exponent,Mod,Acc) -> |
|||
binary_exp_mod((Base*Base) rem Mod, |
|||
Exponent bsr 1,Mod,(Acc * exp_factor(Base,Exponent))rem Mod). |
|||
exp_factor(_,0) -> |
|||
1; |
|||
exp_factor(Base,1) -> |
|||
Base; |
|||
exp_factor(Base,Exponent) -> |
|||
exp_factor(Base,Exponent band 1). |
|||
test() -> |
|||
445 = mod_exp(4,13,497), |
|||
%% Rosetta code example: |
|||
R = 1527229998585248450016808958343740453059 = |
|||
mod_exp(2988348162058574136915891421498819466320163312926952423791023078876139, |
|||
2351399303373464486466122544523690094744975233415544072992656881240319, |
|||
binary_exp(10,40)), |
|||
R. |
|||
% mod module ends here |
|||
%% Modified version of rosetta code entry |
|||
%% Modification was more efficient exponentiation |
|||
%% Modification - use of rpc:pmap to utilise multithreaded CPUs |
|||
-module(miller_rabin). |
|||
-export([is_prime/1,mr_series_test/4,mersennes/1,test/0]). |
|||
is_prime(1) -> false; |
|||
is_prime(2) -> true; |
|||
is_prime(3) -> true; |
|||
is_prime(N) when N > 3, ((N rem 2) == 0) -> false; |
|||
is_prime(N) when ((N rem 2) ==1), N < 341550071728321 -> |
|||
is_mr_prime(N, proving_bases(N)); |
|||
is_prime(N) when ((N rem 2) == 1) -> |
|||
is_mr_prime(N, random_bases(N, 100)). |
|||
proving_bases(N) when N < 1373653 -> |
|||
[2, 3]; |
|||
proving_bases(N) when N < 9080191 -> |
|||
[31, 73]; |
|||
proving_bases(N) when N < 25326001 -> |
|||
[2, 3, 5]; |
|||
proving_bases(N) when N < 3215031751 -> |
|||
[2, 3, 5, 7]; |
|||
proving_bases(N) when N < 4759123141 -> |
|||
[2, 7, 61]; |
|||
proving_bases(N) when N < 1122004669633 -> |
|||
[2, 13, 23, 1662803]; |
|||
proving_bases(N) when N < 2152302898747 -> |
|||
[2, 3, 5, 7, 11]; |
|||
proving_bases(N) when N < 3474749660383 -> |
|||
[2, 3, 5, 7, 11, 13]; |
|||
proving_bases(N) when N < 341550071728321 -> |
|||
[2, 3, 5, 7, 11, 13, 17]. |
|||
is_mr_prime(N, As) when N>2, N rem 2 == 1 -> |
|||
% TStart = erlang:monotonic_time(), |
|||
{D, S} = find_ds(N), |
|||
% elapsed(TStart,"find_ds took ~p.~p seconds~n"), |
|||
%% this is a test for compositeness; the two case patterns disprove |
|||
%% compositeness. |
|||
TestResults = |
|||
rpc:pmap({miller_rabin,mr_series_test},[N,D,S],As), |
|||
R= not lists:any(fun(X) -> X end,TestResults), |
|||
% elapsed(TStart,"is_mr_prime took ~p.~p seconds~n"), |
|||
R. |
|||
mr_series_test(A,N,D,S) -> |
|||
% TMrS = erlang:monotonic_time(), |
|||
R = case mr_series(N, A, D, S) of |
|||
[1|_] -> false; % first elem of list = 1 |
|||
L -> not lists:member(N-1, L) % some elem of list = N-1 |
|||
end, |
|||
% elapsed(TMrS,"mr_series took ~p.~p seconds~n"), |
|||
R. |
|||
%elapsed(TStart,Msg) -> |
|||
% TElapsed_ms = erlang:convert_time_unit(erlang:monotonic_time()-TStart,native,1000), |
|||
% TSec = TElapsed_ms div 1000, |
|||
% Tms = TElapsed_ms rem 1000, |
|||
% io:format(Msg, [TSec,Tms]). |
|||
find_ds(N) -> |
|||
find_ds(N-1, 0). |
|||
find_ds(D, S) -> |
|||
case D rem 2 == 0 of |
|||
true -> |
|||
find_ds(D div 2, S+1); |
|||
false -> |
|||
{D, S} |
|||
end. |
|||
mr_series(N, A, D, S) when N rem 2 == 1 -> |
|||
Js = lists:seq(0, S), |
|||
lists:map(fun(J) -> mod:mod_exp(A, mod:binary_exp(2, J)*D, N) end, Js). |
|||
random_bases(N, K) -> |
|||
[basis(N) || _ <- lists:seq(1, K)]. |
|||
basis(N) when N>2 -> |
|||
% random:uniform returns a single random number in range 1 -> N-3, |
|||
% to which is added 1, shifting the range to 2 -> N-2 |
|||
1 + rand:uniform(N-3). |
|||
mersennes(N) when N>0, is_integer(N) -> |
|||
1 bsl N - 1. |
|||
test() -> |
|||
TStart = erlang:monotonic_time(), |
|||
true = is_prime(7), |
|||
true = is_prime(41), |
|||
false = is_prime(42), |
|||
true = is_prime(mersennes(31)), |
|||
true = is_prime(mersennes(127)), % M(127) checks okay if 64 bit word size exceeded, |
|||
true = is_prime(mersennes(3217)), % about the size of an rsa key, |
|||
TFinish = erlang:monotonic_time(), |
|||
ElapsedSeconds = erlang:convert_time_unit(TFinish - TStart,native,1), |
|||
io:format("Time seconds = ~p~n",[ElapsedSeconds]), |
|||
ok |
|||
. |
|||
</lang> |
|||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |