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
!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 |
This comment has been removed by the author.
ReplyDeleteHi! 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