Skip to content

Commit

Permalink
Attempt p7 in fortran (18)
Browse files Browse the repository at this point in the history
  • Loading branch information
LivInTheLookingGlass committed Oct 3, 2024
1 parent 1a008f6 commit 80912dd
Showing 1 changed file with 13 additions and 11 deletions.
24 changes: 13 additions & 11 deletions fortran/src/include/primes.f90
Original file line number Diff line number Diff line change
@@ -1,14 +1,16 @@
module primes
use constants

implicit none
integer, parameter :: bits_per_int = 32
integer, allocatable :: is_prime(:)
integer :: current_n = 0
integer, parameter :: bits_per_int = 64
integer(i18t), allocatable :: is_prime(:)
integer(i18t) :: current_n = 0
logical :: initialized = .false.

contains

subroutine sieve_up_to(n)
integer, intent(in) :: n
integer(i18t), intent(in) :: n
integer :: p, i

do p = 2, n
Expand All @@ -20,9 +22,9 @@ subroutine sieve_up_to(n)
end do
end subroutine sieve_up_to

recursive integer function next_prime(last) result(ret)
integer, intent(in) :: last
integer :: next
recursive integer(i18t) function next_prime(last) result(ret)
integer(i18t), intent(in) :: last
integer(i18t) :: next

do next = last + 1, current_n
if (get_prime_bit(next)) then
Expand All @@ -36,8 +38,8 @@ recursive integer function next_prime(last) result(ret)
end function next_prime

subroutine expand_sieve(potential_n)
integer, intent(in), optional :: potential_n
integer :: new_size, new_n, i
integer(i18t), intent(in), optional :: potential_n
integer(i18t) :: new_size, new_n, i

if (present(potential_n)) then
new_n = max(potential_n, current_n)
Expand Down Expand Up @@ -65,7 +67,7 @@ end subroutine expand_sieve

! Function to set a bit to 0
subroutine clear_prime_bit(num)
integer, intent(in) :: num
integer(i18t), intent(in) :: num
integer :: i, b
i = (num / bits_per_int) + 1
b = mod(num, bits_per_int)
Expand All @@ -75,7 +77,7 @@ end subroutine clear_prime_bit

! Function to check if a bit is set
logical function get_prime_bit(num) result(bit)
integer, intent(in) :: num
integer(i18t), intent(in) :: num
bit = logical(btest(is_prime(num / bits_per_int + 1), mod(num, bits_per_int)))
end function get_prime_bit

Expand Down

0 comments on commit 80912dd

Please sign in to comment.