HRÓÐÓLFR
HRÓÐÓLFR

Reputation: 5992

Prime generator in Fortran

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

Answers (2)

Ondřej Čertík
Ondřej Čertík

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

Henrik
Henrik

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

Related Questions