CODEWITHSUNDEEP
algorithmfortranprimes
Andro
Andro

Reputation: 2232

searching a prime number

I hope I am not duplication any question but the suggested topic did not provide with any similar problem. I have a function that check if a number is a prime one. Now this is the slowest possible way to search for a prime.

subroutine is_prime_slow(num, stat)
        implicit none
        logical :: stat
        integer :: num
        integer :: i
        if ((num .le. 3) .and. (num .gt. 1)) then
            stat = .true.
            return
        end if

        ! write(*,*) 'limit = ',limit
        do i = 2,num - 1
            ! write(*,*) 'mod(',limit,i,') = ',mod(limit,i)
            if (mod(num,i) == 0) then
                stat = .false.
                return
            end if
        end do
        stat = .true.   
        return     
    end

Now let's say that I do some improvement to it. 

subroutine is_prime_slow(num, stat)
    implicit none
    logical :: stat
    integer :: num
    integer :: i
    if ((num .le. 3) .and. (num .gt. 1)) then
        stat = .true.
        return
    end if
    ! IMPROVEMENT
    if ((mod(num,2) == 0) .or. (mod(num,3) == 0) .or. (mod(num,5) == 0) .or. (mod(num,7) == 0)) then
        stat = .false.
        return
    end if

    ! write(*,*) 'limit = ',limit
    do i = 2,num - 1
        ! write(*,*) 'mod(',limit,i,') = ',mod(limit,i)
        if (mod(num,i) == 0) then
            stat = .false.
            return
        end if
    end do
    stat = .true.   
    return     
end

Now the second version excludes more than half of numbers e.g. all that are divisible by 2,3,5,7. How is it possible that when I time the execution with the linux 'time' program, the 'improved' version performs just as slowly? Am I missing some obvious trick?

Searching the first 900000 numbers:
1st: 4m28sec
2nd  4m26sec

Upvotes: 2

Views: 1129

Answers (2)

Joni
Joni

Reputation: 111339

The multiples of 2, 3, 5, and 7 are quickly rejected by the original algorithm anyway, so jumping over them does not improve the performance at all. Where the algorithm spends most of its time is in proving that numbers with large prime factors are composite. To radically improve the performance you should use a better primality test, such as Miller-Rabin.

A simpler improvement is testing factors only up to sqrt(num), not num-1. If that doesn't make immediate sense, think about how big the smallest prime factor of a composite number can be. Also, if you are looking for primes from 1 to N, it may be more efficient to use a prime number sieve, or testing divisibility only by primes you have already found.

Upvotes: 7

Alexander Vogt
Alexander Vogt

Reputation: 18108

I just recently coded something similar ;-)

! Algorithm taken from https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
subroutine eratosthenes_sieve(n, primes)
  implicit none
  integer,intent(in)    :: n
  integer,allocatable,intent(out)   :: primes(:)
  integer               :: i, j, maxPrime, stat
  logical               :: A(n)

  maxPrime = floor(sqrt(real(n)))
  A = .true.

  do i=2,maxPrime
    j = i*i
    do
      A(j) = .false.
      j = j + i ; if ( j .gt. n ) exit
    enddo
  enddo !i

  allocate( primes( count(A)-1 ), stat=stat )
  if ( stat /= 0 ) stop 'Cannot allocate memory!'

  j = 1
  do i=2,n ! Skip 1
    if ( .not. A(i) ) cycle
    primes( j ) = i
    j = j + 1 ; if ( j > size(primes) ) exit
  enddo
end subroutine

This algorithm gives you all prime numbers up to a certain number, so you can easily check whether your prime is included or not: 

if ( any(number == prime) ) write(*,*) 'Prime found:',number

Upvotes: 2

Related Questions