How to Find the First n Prime Numbers

A simple method to find a prime number is to test the divisibility of the number by all numbers that come before it. However, this is a tedious process. To make it faster, we must eliminate numbers with which we do not need to check the divisibility. We also need to eliminate even numbers after 2, which are not prime for sure. Let us see how to find primes easily.

  •  Firstly, if we have a test number, we only need to check it's divisibility by the prime numbers before it. Any composite factors will be detected therein. Not that we are looking for any, because the flag will tell the first time a number is divisible by some other and we will exit the loop.
  • Further, we do not need to check the number's divisibility by any number greater than its square root. This is because the square root (since we are talking about primes here) could be it's highest possible prime factor. If not, well we already said we are not interested in composites.
  • We need to specify the first three primes for which we will not test any divisibility. This makes the job easier.
  • If we detect a prime, that is if a number is not found to be divisible according to the tests in the previous steps, then we assign to it the next place in the array and increase the count of primes.
  • Notice that prime numbers are of the form, 6n+1,6n+5 (also called (6(n+1)-1) or simply 6N-1. This does not mean that all the numbers in this form are prime. The numbers in this form are not divisible by 2 and 3. And since the multiples of 2 and 3 occur more often than multiples of higher numbers,eliminating these multiples makes our 'prime hunt' easier.
Now for the Fortran Program, 

 !To find the first n prime numbers
program primenum
    implicit none
    integer prime(99999),i,nprime,testnum,incr,maxdiv,rem,flag,div,n
    write (*,*) "How many prime numbers do you wish to find?"
    read (*,*) n
    write (*,*) "The first ",n," prime numbers are"

    !To read the first three numbers of the array
    prime(1)=2
    prime(2)=3
    prime(3)=5

    !Initialize count
    nprime=3
  
    !To test divisibility of primes by numbers less than maxdiv, which is the sqrt of testnum
    testnum=7
    incr=2
  
    do   
        maxdiv=sqrt(real(testnum))
          flag=1
      
          do i=3,nprime
              div=prime(i)
             rem=mod(testnum,div) !To check divisibility only by primes
            if (rem==0) then
                flag=0
                  exit
            end if

            if (div<(maxdiv+1)) cycle
           end do
            
        !To increase count on detecting a prime
           if (flag==1) then
               nprime=nprime+1
            prime(nprime)=testnum
        end if
          
        !To terminate on fetching the desired number of primes
            if (nprime>n) exit
            
        !To test more odd numbers not divisible by 2,3
           incr=6-incr
  
        !To set the next number to test
        testnum=testnum+incr
    end do
            
    !To display the numbers. Don't get confused by the formatting. You can use your own style
    do i=1,n,5
      write(*,10) "p(",i,")=",prime(i),"p(",i+1,")=",prime(i+1),"p(",i+2,")=",prime(i+2),"p(",i+3,")=",prime(i+3),&
      &
"p(",i+4,")=",prime(i+4) !Ampersands are used to carry a long statement to the next line
      10 format (a2,i3,a2,i4,3x,a2,i3,a2,i4,3x,a2,i3,a2,i4,3x,a2,i3,a2,i4,3x,a2,i3,a2,i4,3x)
      end do
      end program primenum
   

Sample Output: The First 250 Prime Numbers

2 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. Hi! I tried exactly the same program in fortran 90, but it isn't working, Could you please help me? by here or by email? Thanks :)

    ReplyDelete