Reputation: 5992
I'm just trying to get a little familiar with Fortran (because I can) so I wrote this small program that uses the Sieve of Eratosthenes to generate primes. Here is the program:
program prime
implicit none
integer num_primes, at, found, i
logical is_prime
integer, allocatable, dimension(:) :: primes ! array that will hold the primes
print *, "How many primes would you like to find?"
read (*, *) num_primes
allocate (primes(num_primes))
primes(1) = 2
at = 2
found = 1
do
is_prime = .true. ! assume prime
do i = 1, found
if (modulo(at, primes(i)) == 0) then ! if divisible by any other element
is_prime = .false. ! in the array, then not prime.
at = at + 1
continue
end if
end do
found = found + 1
primes(found) = at
print *, at
at = at + 1
if (found == num_primes) then ! stop when all primes are found
exit
endif
end do
end program prime
Running this will show what the error is, for example, trying to find 10 primes will yield the following numbers: 3, 5, 7, 11, 13, 16, 17, 19, 23. Obviously 16 is not prime. What could be wrong?
Upvotes: 0
Views: 8524
Reputation: 830
Here is the working program by implementing Henrik suggestion (and then I did a few formatting changes):
program prime
implicit none
integer :: num_primes, at, found, i
logical :: is_prime
integer, allocatable, dimension(:) :: primes ! array that will hold the primes
print *, "How many primes would you like to find?"
read(*, *) num_primes
allocate(primes(num_primes))
primes(1) = 2
at = 2
found = 1
do
is_prime = .true. ! assume prime
do i = 1, found
if (modulo(at, primes(i)) == 0) then ! if divisible by any other element
is_prime = .false. ! in the array, then not prime.
exit
end if
end do
if (is_prime) then
found = found + 1
primes(found) = at
print *, at
end if
at = at + 1
if (found == num_primes) then ! stop when all primes are found
exit
end if
end do
end program prime
When you run this, you'll get:
How many primes would you like to find?
20
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
and the primes
array will contain all the prime numbers. Is that what you wanted?
Upvotes: 2
Reputation: 23324
This is not the Sieve of Eratosthenes, which doesn't need modulo
.
It's sort of trial division, but when you find a divisor, you increment at
by 1 and start testing with the next prime where you should re-start.
When you find 15 is divisable by 3, at
is incremented to 16 and tested against 5, 7, an 11, which are of course not divisors of 16 and so 16 is printed.
I'm not a FORTRAN programmer, but I think you should replace the lines
at = at + 1
continue
by
exit
and execute the lines
found = found + 1
primes(found) = at
print *, at
only if is_prime is true.
Upvotes: 3