From ea546abe005ec131950e97ba9f0a173b74544266 Mon Sep 17 00:00:00 2001
From: Olivia Appleton <olivia.kay.appleton@gmail.com>
Date: Mon, 30 Sep 2024 23:02:53 -0500
Subject: [PATCH] Attempt p7 in fortran (2)

---
 docs/src/fortran/lib/primes.rst |  4 +--
 fortran/src/include/primes.f90  | 43 ++++++++++++++++++---------------
 fortran/src/p0007.f90           |  2 +-
 3 files changed, 27 insertions(+), 22 deletions(-)

diff --git a/docs/src/fortran/lib/primes.rst b/docs/src/fortran/lib/primes.rst
index 7debce75..bec38539 100644
--- a/docs/src/fortran/lib/primes.rst
+++ b/docs/src/fortran/lib/primes.rst
@@ -11,9 +11,9 @@ View source code :source:`fortran/include/primes.f90`
     :returns answer: The next prime in the sequence
     :rtype answer: integer
 
-.. f:function:: expand_sieve(potential_n)
+.. f:subroutine:: expand_sieve(potential_n)
 
-    A helper function that lets you pre-initialize the cache of primes to a high number.
+    A helper routine that lets you pre-initialize the cache of primes to a high number.
     Otherwise you incur a cost of ``O(n loglogn)`` each time ``n`` doubles.
 
     :p integer potential_n: The size you wish to expand the sieve to
diff --git a/fortran/src/include/primes.f90 b/fortran/src/include/primes.f90
index d0551726..9d864821 100644
--- a/fortran/src/include/primes.f90
+++ b/fortran/src/include/primes.f90
@@ -3,34 +3,38 @@ module primes
     integer, parameter :: bits_per_int = 32
     integer, allocatable :: is_prime(:)
     integer :: current_n = 0
+    logical :: initialized = .false.
 
 contains
 
-    function initialize_sieve(n)
+    subroutine initialize_sieve(n)
         integer, intent(in) :: n
         integer :: size
         ! Calculate size of the bit array (integer array to hold bits)
-        size = (n / bits_per_int) + 1
-        allocate(is_prime(size))
-        is_prime = -1  ! Initialize all bits to 1 (all numbers are initially prime)
-        clear_prime_bit(0)
-        clear_prime_bit(1)
-        current_n = n
-        sieve_up_to(current_n)
-    end function initialize_sieve
+        if (.not. initialized) then
+            size = (n / bits_per_int) + 1
+            allocate(is_prime(size))
+            is_prime = -1  ! Initialize all bits to 1 (all numbers are initially prime)
+            call clear_prime_bit(0)
+            call clear_prime_bit(1)
+            current_n = n
+            initialized = .true.
+            call sieve_up_to(current_n)
+        end if
+    end subroutine initialize_sieve
 
-    function sieve_up_to(n)
+    subroutine sieve_up_to(n)
         integer, intent(in) :: n
         integer :: p, limit
         limit = sqrt(real(n))
         do p = 2, limit
             if (get_prime_bit(p)) then
                 do i = p * p, n, p
-                    clear_prime_bit(i)
+                    call clear_prime_bit(i)
                 end do
             end if
         end do
-    end function sieve_up_to
+    end subroutine sieve_up_to
 
     recursive integer function next_prime(last)
         integer, intent(in) :: last
@@ -43,14 +47,15 @@ recursive integer function next_prime(last)
             end if
         end do
 
-        expand_sieve()
+        call expand_sieve()
         next_prime = next_prime(next)  ! Recursively find the next prime
     end function next_prime
 
-    function expand_sieve(potential_n)
+    subroutine expand_sieve(potential_n)
         integer, intent(in), optional :: potential_n
         integer :: new_size, new_n, i
 
+        call initialize_sieve()
         if (present(potential_n)) then
             new_n = max(potential_n, current_n)
         else
@@ -60,17 +65,17 @@ function expand_sieve(potential_n)
         new_size = (new_n / bits_per_int) + 1
         allocate(is_prime(new_size))
         is_prime = -1  ! All bits set to 1
-        clear_prime_bit(is_prime, 0)
-        clear_prime_bit(is_prime, 1)
+        call clear_prime_bit(is_prime, 0)
+        call clear_prime_bit(is_prime, 1)
         sieve_up_to(new_n)
         current_n = new_n
-    end function expand_sieve
+    end subroutine expand_sieve
 
     ! Function to set a bit to 0
-    function clear_prime_bit(num)
+    subroutine clear_prime_bit(num)
         integer, intent(in) :: num
         is_prime(num / bits_per_int) = bitand(is_prime(num / bits_per_int), ieor(0xFFFFFFFF, 1 << (mod(num, bits_per_int))))
-    end function clear_prime_bit
+    end subroutine clear_prime_bit
 
     ! Function to check if a bit is set
     logical function get_prime_bit(num) result(bit)
diff --git a/fortran/src/p0007.f90 b/fortran/src/p0007.f90
index 7b265e7a..d0b7ec31 100644
--- a/fortran/src/p0007.f90
+++ b/fortran/src/p0007.f90
@@ -15,7 +15,7 @@ integer function p0007() result(answer)
         integer :: i
 
         answer = 1
-        expand_sieve(2**17)
+        call expand_sieve(2**17)
         do i = 1, 10001
             answer = next_prime(answer)
         end do