/* * if n < 1,373,653, it is enough to test a = 2 and 3; * if n < 9,080,191, it is enough to test a = 31 and 73; * if n < 4,759,123,141, it is enough to test a = 2, 7, and 61; * if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803; * if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11; * if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13; * if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17. */

// There is no built-in method for generating random BigInteger values.// Instead, random BigIntegers are constructed from randomly generated// byte arrays of the same length as the source. RandomNumberGenerator rng = RandomNumberGenerator.Create();byte[] bytes =newbyte[source.ToByteArray().LongLength]; BigInteger a;

for(int i =0; i < certainty; i++){do{// This may raise an exception in Mono 2.10.8 and earlier.// http://bugzilla.xamarin.com/show_bug.cgi?id=2761 rng.GetBytes(bytes); a =new BigInteger(bytes);}while(a <2|| a >= source -2);

is_mr_prime(N,As)whenN>2,Nrem2==1->{D,S}=find_ds(N),%% this is a test for compositeness; the two case patterns disprove%% compositeness.notlists:any(fun(A)->casemr_series(N,A,D,S)of[1|_]->false;% first elem of list = 1L->notlists:member(N-1,L)% some elem of list = N-1endend,As).

! this is limited since we're not using a bignum libcall do_test((/(i, i=1, 29)/))

contains

subroutine do_test(a)integer(huge), dimension(:), intent(in)::a

integer::i

do i =1, size(a,1) print *, a(i), miller_rabin_test(a(i), prec)enddo

endsubroutine do_test

endprogram TestMiller

Possible improvements: create bindings to the GMP library, change integer(huge) into something like type(huge_integer), write a lots of interfaces to allow to use big nums naturally (so that the code will be unchanged, except for the changes said above)

Integer overflow is a severe risk, and even 64-bit integers won't get you far when the formulae are translated as MOD(A**D,N) - what is needed is a method for raising to a power that incorporates the modulus along the way. There is no library routine for that, so...

In this run, 32-bit integers falter for 19 in calculating 179, and 64-bit integers falter for 31 with 1915 by showing a negative number. Other 64-bit overflows however do not show a negative (as with 2315) because there is about an even chance that the high-order bit will be on or off. The compiler option for checking integer overflow does not report such faults with 64-bit integers, at least with the Compaq V6.6 F90/95 compiler. In this context, one misses the IF OVERFLOW ... that was part of Fortran II but which has been omitted from later versions.

Thus, there is no avoiding a special MODEXP function, even for small test numbers.

Go has it in math/big in standard library as ProbablyPrime. The argument n to ProbablyPrime is the input k of the pseudocode in the task description.

Deterministic

Below is a deterministic test for 32 bit unsigned integers. Intermediate results in the code below include a 64 bit result from multiplying two 32 bit numbers. Since 64 bits is the largest fixed integer type in Go, a 32 bit number is the largest that is convenient to test.

The main difference between this algorithm and the pseudocode in the task description is that k numbers are not chosen randomly, but instead are the three numbers 2, 7, and 61. These numbers provide a deterministic primality test up to 2^32.

/* Miller-Rabin primality test. For n < 341550071728321, the test is deterministic; for larger n, the number of bases tested is given by the option variable primep_number_of_tests, which is used by Maxima in primep. The bound for deterministic test is also the same as in primep. */

This naive implementation of a Miller-Rabin method is adapted from
the Prolog version on this page. The use of the form integer(N) to
permit use of integers of arbitrary precision as done here is not
efficient. It is better to define a tabled version of each known
integer and to use the tabled versions. For example, suppose you want
integer(2), then do

and use n2 as the integer in your code. Performance will be greatly
improved. Also Mercury has a package using Tom's Math for integers of
arbitrary precision and another package to some of the functions of the
GMP library for much faster operation with long integers. These can be
found with instructions for use in Github.

A basic deterministic test can be obtained by an appeal to the ERH (as proposed by Gary Miller) and a result of Eric Bach (improving on Joseph Oesterlé). Calculations of Jan Feitsma can be used to speed calculations below 264 (by a factor of about 250).

While normally one would use is_prob_prime, is_prime, or is_provable_prime, which will do a BPSW test and possibly more, we can use just the Miller-Rabin test if desired. For large values we can use an object (e.g. bigint, Math::GMP, Math::Pari, etc.) or just a numeric string.

use ntheory qw/is_strong_pseudoprime miller_rabin_random/;sub is_prime_mr {my$n=shift;# If 32-bit, we can do this with 3 bases.return is_strong_pseudoprime($n,2,7,61)if($n>>32)==0;# If 64-bit, 7 is all we need.return is_strong_pseudoprime($n,2,325,9375,28178,450775,9780504,1795265022)if($n>>64)==0;# Otherwise, perform a number of random base tests, and the result is a probable prime test.return miller_rabin_random($n,20);}

Math::Primality also has this functionality, though its function takes only one base and requires the input number to be less than the base.

Math::Pari can be used in a fashion similar to the Pari/GP custom function. The builtin accessed using a second argument to ispseudoprime was added to a later version of Pari (the Perl module uses version 2.1.7) so is not accessible directly from Perl.

This versions will give answers with a very small probability of being false. That probability being dependent on _mrpt_num_trials and the random numbers used for name a passed to function try_composite.

importrandom

_mrpt_num_trials =5# number of bases to test

def is_probable_prime(n):""" Miller-Rabin primality test.

A return value of False means n is certainly not prime. A return value of True means n is very likely a prime.

# test the base a to see whether it is a witness for the compositeness of ndef try_composite(a):ifpow(a, d, n)==1:returnFalsefor i inrange(s):ifpow(a,2**i * d, n)== n-1:returnFalsereturnTrue# n is definitely composite

for i inrange(_mrpt_num_trials): a =random.randrange(2, n)if try_composite(a):returnFalse

This versions will give correct answers for n less than 341550071728321 and then reverting to the probabilistic form of the first solution. By selecting predetermined values for the a values to use instead of random values, the results can be shown to be deterministically correct below certain thresholds.
For 341550071728321 and beyond, I have followed the pattern in choosing a from the set of prime numbers.While this uses the best sets known in 1993, there are better sets known, and at most 7 are needed for 64-bit numbers.