Trying to implement the Miller-Rabin test in Fortran but encountered a bug

Help
1

Hi there, I am trying to implement the Miller-Rabin test for a prime number in Fortran, however I ran into a weird instance (bug? maybe a mistake on my part ?)
Where the number 29 turns up as a composite when it should be showing up as a prime.
Compiler: ifort
optimization: -O3
Here’s the code:

PROGRAM pr3
  IMPLICIT NONE

  INTEGER(KIND=8),PARAMETER :: given_number=29
  INTEGER(KIND=8) :: number_of_primes
  INTEGER(KIND=8),ALLOCATABLE :: prime_array(:)


  IF (is_prime(given_number) .EQ. .TRUE.) THEN
     PRINT *, "The given number ",given_number," was prime "
  ELSE
     PRINT *, "The number was composite"
  END IF




CONTAINS

  LOGICAL FUNCTION is_prime(number) RESULT(test_result)
    INTEGER(KIND=8),INTENT(IN) :: number
    INTEGER(KIND=8) :: test_number

    test_result = .TRUE.
    test_number=2

    DO WHILE (test_number .LT. INT(CEILING(SQRT(REAL(number,8))),8))
       IF (MOD(test_number**number,number) .EQ. test_number) THEN
          test_result = .TRUE.
       ELSE
          test_result = .FALSE.
       END IF
       test_number=test_number+1
    END DO

  END FUNCTION is_prime

END PROGRAM pr3
2

Additional info about this:
When the test_number=5 the above operation should return 29, instead it returns 26 which isn’t equal to the given number 29

3

Your mistake is to think that test_numbernumber can be represented as a 64 bits integer. These numbers quickly become very large and the integer arithmetic then causes them to be truncated. Which leads to the wrong results. Try printing test_numbernumber and real(test_number)**number.

4

Oh, sorry about the bold text - the double asterisk gets in the way :frowning:

5

If you use back ticks (`) around the code snippet, the double asterisks are not treated as markdown syntax (i.e. bold markers).

6 
                     2   
  2.1474836E+09          
                     3   
  6.1767341E+14          
                     4   
  4.6116860E+18          
                     5   
  8.0333668E+18          
 The number was composite

The numbers are quite huge :frowning:

7

Yes, should have thought about that.

8

@Arjen
Unfortunately I don’t know how to remedy this . What do you recommend I do ?

9

One way would be to calculate mod(test_number**number, number) as:

residue = test_number
do i = 1,number-1
    residue = mod( testnumber*residue, number)
enddo

(Check the number of multiplications though!) This way you avoid the very large numbers while still computing the right modulo. And of course you could use a recursive algorithm as well to reduce the number of mulitplications to log(N) in stead of N. This under the condition that I got the maths right :slight_smile: .

1 Like
10

Apologies, I am a little confused. What is happening inside the do loop ?

11

This is just a language comment, not about the algorithm. In fortran, one would write this simply as

test_result = MOD(test_number**number,number) .EQ. test_number

There is no need for the if-then-else-endif construct.

1 Like
12

This is neat ! Never knew you could do this. Thank you