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?
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.
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():
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
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.
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.
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
Very nice implementation
Thank you very much
can i use this in my repository?
There is no copyright!
@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.
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.
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
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 )
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.
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.
Also, I agree simplicity and readability are very very important!
Yeah, KISS ! (Keep it simple, stupid: https://en.wikipedia.org/wiki/KISS_principle)
It was said before UNIX birth, but it is now considered a part of the UNIX philosophy:
Many interesting and beautiful principles than can inspire us.