Hi all,
I would like to implement sieve of eratosthenes algorithm to generate prime numbers in a specific interval. To do this we need to delete some of elements of array in each step in order to remain an array with all prime numbers.
Is there a way to do this?
1 Like
Hi @ELNS!
For this and other similar problems I would recommend using a masking approach whereby you maintain a logical array of equivalent length and you use the boolean value at each position in the logical array to indicate whether you have ‘deleted’ an item in the original array.
You can then iterate over the two arrays at the end to reveal only the non-masked elements.
As you might know, there are data structures that allow element deletion easily like linked lists, however they have performance overheads and they give no benefits for your particular task IMO.
2 Likes
Hi @ELNS
you will generally use a LOGICAL array, initialized with .TRUE., and eliminate the multiples of each number by writing .FALSE. where needed.
Be aware that in Fortran, a LOGICAL will probably use 4 bytes on your PC. If you want to work on an array with a size similar to the RAM of your PC, you should use a LOGICAL(1), which will probably be stored in one byte.
If you want to go higher, you will have to optimize the memory usage by using instructions to manipulate bits in INTEGER. See for example IBSET():
https://gcc.gnu.org/onlinedocs/gfortran/IBSET.html
Also: https://stackoverflow.com/questions/14143677/how-to-implement-a-bit-array-in-modern-fortran
2 Likes
I wrote code below:
program main
implicit none
integer :: i , j, k
logical , dimension(999) :: a = [(.true., i = 2, 1000)]
do i = 2, size(a)
if (i ** 2 > size(a)) then
exit
elseif (a(i) .eqv. .true.) then
do j = 2, size(a)
if (a(i * j) .eqv. .true.) then
a(i * j) = .false.
else
cycle
end if
end do
else
cycle
end if
end do
do k = 2, size(a)
if (a(k) .eqv. .true.) then
print *, k
else
cycle
end if
end do
end program main
But after i an that, i got this error:
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x7F4E6D5D2CB0
#1 0x7F4E6D5D1EB0
#2 0x7F4E6CE4A3DF
#3 0x4011E7 in MAIN__ at Exercise_11.f90:?
Segmentation fault (core dumped)I think it would be clearer and so safer to simply write:
logical , dimension(999) :: a
a = .true.
(is 1 a prime ? Yes, no… 
)
And be careful, if you put a print *, i*j you will see that i*j can be 1000.
Another thing, the three pieces of code here are useless because at the end of the loops:
else
cycle
end if
And personally I would compute the sqrt of size(a) to avoid the first if.
I think you can also write if (a(i)) then instead of if (a(i) .eqv. .true.) then, a being a boolean.
Try to write the simplest possible code, it will be easier to find the bugs. Here your code is most complicated than it should and so more difficult to read.
1 Like
Here is my solution:
program main
implicit none
integer :: i, j, k, root
logical(1), dimension(1000) :: a
a = .true.
a(1) = .false.
root = int(sqrt(real(size(a))))
do i = 2, root
if (a(i)) then
do j = i*i, size(a), i
a(j) = .false.
end do
end if
end do
do k = 1, size(a)
if (a(k)) then
print *, k
end if
end do
end program main
Just a doubt about root = int(sqrt(real(size(a)))) (will it be OK in every case ?)
Note that with do j = i*i, size(a), i you are sure you won’t go out of bounds.
1 Like
Seems to run until 2e9.
If you want to go beyond the maximum positive value of integer(4), you will need some thinking.
Note also that you can simply write:
do k = 1, size(a)
if (a(k)) print *, k
end do
1 Like
There is no copyright! 
1 Like
@ELNS I believe this: https://godbolt.org/z/yfqhQ8 online compiler explorer would be very useful to you, and as a matter of fact, to everyone. I refer to it almost on a daily basis. The compiler’s output is updated live as you type the code. On the right upper box, you see the assembly generated code and on the right lower box the output. You can also insert flags. By just adding -fcheck=all you immediately get a much more useful error message. There are other features that you can play with. Most importantly it supports all gfortran versions along with other compilers that you can test against each other.
ps. Fortran is just one of the many languages it supports.
3 Likes
Thank Stavros,
This is extremely useful.
I guess I am bit late to this party, but here is my 5 cents!
I did this quite a while ago in my little project. I will just copy my implementation here, which is almost the same with @vmagnin’s version but using allocatable arrays instead.
subroutine sieve_of_Eratosthenes_int32(n, prime_arr)
integer, intent(in) :: n
logical, allocatable, dimension(:) :: prime_arr
integer :: i, j
allocate( prime_arr(0:n) )
prime_arr = .true.
prime_arr(0:1) = .false.
do i = 2, floor( sqrt( real(n, dp) ) )
if ( prime_arr(i) ) then
j = i * i
do while ( j <= n )
prime_arr(j) = .false.
j = j + i
end do
end if
end do
end subroutine sieve_of_Eratosthenes_int32
Then you can write an int64 version and wrap them with an interface. But Sieve of Eratosthenes can be slow as well (at least that’s my experience when solving this problem), when dealing big prime numbers try the Miller-Rabin test. But then you will either have to write a big integer library yourself or write a wrapper to use GMP/MPFR (which is why Fortran standard library is sooooo needed!).
One word about your initialization style
logical , dimension(999) :: a = [(.true., i = 2, 1000)]
this will likely create an temporary array and slow down your code, for example, if you compile it by gfortran with -Warray-temporaries. You can do this as well
! declaration
logical :: a(2:1000)
! initialization
a = .true.
3 Likes
I too am a bit late to the party here. I guess I’d prefer not to use a mask and to use the pack intrinsic which is a way of deleting elements in an array (actually it creates a new array from a mask). I’m not claiming this will give you the best performance but it is I think a bit simpler. I also prefer not to use the square root but just exit the loop when there’s nothing more to do.
program sieve
! Tested with gfortran 10.1.0
! gfortran -o sieve sieve.f90 -std=f2018 -fimplicit-none -fdefault-real-8 -Wall -Wextra
implicit none
integer, parameter :: n = 30
integer :: i, j
integer :: numbers(n) = [(i,i=1,n)]
integer, allocatable :: primes(:)
numbers(1) = 0
i = 2
do
if (numbers(i) > 0) then
do j = i*i, n, i
numbers(j) = 0
end do
end if
i = i+1
if (i > n) exit
end do
primes = pack(numbers,numbers>0)
print *, primes
stop
end program sieve
5 Likes
There is a decrease in performance, but from the table below (quoted from this book), the larger the size of the data the smaller the difference between two methods. Also, I agree simplicity and readability are very very important! (People in this community should have been suffered more or less from some poorly written legacy FORTRAN codes 
)
2 Likes
Welcome to the Discourse @simong and thanks for your post!
My interpretation of your solution is that it is still a masking approach but a more memory-efficient one by simply not using an additional logical array. The pack intrinsic is a natural choice for filtering the masked array at the end.
2 Likes
You’re quite right this is a masking method - I’ve just used the numbers array as its own mask so that the algorithm is not as obscured by the implementation. Of course you can’t delete an element of an array in Fortran any more than you can in C so some sort of mask will always be needed. To do away with that a linked-list could be used but that does seem like a lot more effort.
2 Likes
Also, I agree simplicity and readability are very very important!
Yeah, KISS ! (Keep it simple, stupid: KISS principle - Wikipedia)
It was said before UNIX birth, but it is now considered a part of the UNIX philosophy:
Unix philosophy - Wikipedia
Many interesting and beautiful principles than can inspire us.
Beauty matters!
5 Likes