Miller–Rabin primality test

Miller–Rabin primality test
You are encouraged to solve this task according to the task description, using any language you may know.

This page uses content from Wikipedia. The original article was at Miller–Rabin primality test. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not. The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm.

The pseudocode, from Wikipedia is:

Input: n > 2, an odd integer to be tested for primality;
       k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
LOOP: repeat k times:
   pick a randomly in the range [2, n − 2]
   x ← ad mod n
   if x = 1 or x = n − 1 then do next LOOP
   for r = 1 .. s − 1
      x ← x2 mod n
      if x = 1 then return composite
      if x = n − 1 then do next LOOP
   return composite
return probably prime

• The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
• Deterministic variants of the test exists and can be implemented as extra (not mandatory to complete the task)

Fortran

For the module PrimeDecompose, see Prime decomposition.

<lang fortran>module Miller_Rabin

 use PrimeDecompose
 !use iso_fortran_env
 implicit none

 integer, parameter :: max_decompose = 100

 private :: int_rrand, max_decompose

contains

 function int_rrand(from, to)
   integer(huge) :: int_rrand
   integer(huge), intent(in) :: from, to

   real :: o
   call random_number(o)
   int_rrand = floor(from + o * real(max(from,to) - min(from, to)))
 end function int_rrand

 function miller_rabin_test(n, k) result(res)
   logical :: res
   integer(huge), intent(in) :: n
   integer, intent(in) :: k
   
   integer(huge), dimension(max_decompose) :: f
   integer(huge)                     :: s, d, i, a, x, r

   res = .true.
   f = 0

   if ( (n <= 2) .and. (n > 0) ) return
   if ( mod(n, 2) == 0 ) then
      res = .false.
      return
   end if

   call find_factors(n-1, f)
   s = count(f == 2)
   d = (n-1) / (2 ** s)
   loop:  do i = 1, k
      a = int_rrand(2_huge, n-2)
      x = mod(a ** d, n)
      
      if ( x == 1 ) cycle
      do r = 0, s-1
         if ( x == ( n - 1 ) ) cycle loop
         x = mod(x*x, n)
      end do
      if ( x == (n-1) ) cycle
      res = .false.
      return
   end do loop
   res = .true.
 end function miller_rabin_test

end module Miller_Rabin</lang>

Testing

<lang fortran>program TestMiller

 use Miller_Rabin
 implicit none

 integer, parameter :: prec = 30
 integer(huge) :: i

 ! this is limited since we're not using a bignum lib
 call 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)
   end do

 end subroutine do_test
 

end program TestMiller</lang>

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)