sieve_of_eratosthenes Function

    pure function sieve_of_eratosthenes(N) result(primes)
        integer(isp), intent(in) :: N
        integer(isp), dimension(:), allocatable :: primes
        logical, dimension(:), allocatable :: is_prime
        integer(isp) :: i, j, count_primes, limit, idx

        allocate (is_prime(0:N))
        is_prime = .true.
        is_prime(0:1) = .false.

        limit = int(sqrt(real(N, kind=sp)))

        do i = 2, limit
            if (is_prime(i)) then
                do j = i * i, N, i
                    is_prime(j) = .false.
                end do
            end if
        end do

        count_primes = count(is_prime)
        allocate (primes(count_primes))
        idx = 1
        do i = 2, N
            if (is_prime(i)) then
                primes(idx) = i
                idx = idx + 1
            end if
        end do
        deallocate (is_prime)
    end function sieve_of_eratosthenes