primes.f90

View source code here on GitHub!

function  primes/next_prime(last)
Parameters:

last [integer(i18t)] :: The last prime you received (or 1)

Return:

answer [integer(i18t)] :: The next prime in the sequence

Call to:

expand_sieve()

subroutine  primes/prime_factor(num, factor)
Parameters:

  • num [integer(i18t),intent(inout)] :: The number you want to factor

  • factor [integer(i18t),intent(out)] :: The factor produced

Call to:

next_prime(), expand_sieve()

function  primes/is_composite(n)
Parameters:

n [integer(i18t)] :: The number you want to determine the compositeness of

Return:

answer [integer(i18t)] :: The first prime factor if composite, or 0 otherwise

Call to:

prime_factor()

function  primes/is_prime(n)
Parameters:

n [integer(i18t)] :: The number you want to determine the primality of

Return:

answer [logical]

Call to:

prime_factor()

subroutine  primes/expand_sieve(potential_n)

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.

Note

If you call this with a value lower than the current size, there is no cost, and I highly encourage you to do so.

Options:

potential_n [integer(i18t),optional] :: The size you wish to expand the sieve to

Called from:

next_prime(), prime_factor()

  1module primes
  2    use constants
  3    implicit none
  4    integer(i18t), parameter :: bits_per_int = 64
  5    integer(i18t), allocatable :: prime_sieve(:)
  6    integer(i18t) :: current_n = 0
  7contains
  8    subroutine sieve_up_to(n)
  9        integer(i18t), intent(in) :: n
 10        integer(i18t) :: p, i
 11
 12        do p = 2, n
 13            if (get_prime_bit(p)) then
 14                do i = p * p, n, p
 15                    call clear_prime_bit(i)
 16                end do
 17            end if
 18        end do
 19    end subroutine sieve_up_to
 20
 21    recursive integer(i18t) function next_prime(last) result(ret)
 22        integer(i18t), intent(in) :: last
 23        integer(i18t) :: next
 24
 25        do next = last + 1, current_n
 26            if (get_prime_bit(next)) then
 27                ret = next
 28                return
 29            end if
 30        end do
 31
 32        call expand_sieve()
 33        ret = next_prime(next)
 34    end function next_prime
 35
 36    subroutine expand_sieve(potential_n)
 37        integer(i18t), intent(in), optional :: potential_n
 38        integer(i18t) :: new_size, new_n
 39
 40        if (present(potential_n)) then
 41            new_n = potential_n
 42        else
 43            new_n = current_n * 2
 44        end if
 45
 46        if (new_n <= current_n) then
 47            return
 48        endif
 49
 50        new_size = (new_n / bits_per_int) + 1
 51        ! recalculate n to fit boundaries
 52        new_n = bits_per_int * (new_size - 1)
 53        if (allocated(prime_sieve)) then
 54            deallocate(prime_sieve)
 55        end if
 56        allocate(prime_sieve(new_size))
 57        if (.not. allocated(prime_sieve)) then
 58            print *, "Memory allocation failed for prime sieve. Exiting."
 59            stop ERROR_ALLOCATE_FAILED
 60        end if
 61        prime_sieve = -1
 62        call clear_prime_bit(0_i18t)
 63        call clear_prime_bit(1_i18t)
 64        call sieve_up_to(new_n)
 65        current_n = new_n
 66    end subroutine expand_sieve
 67
 68    ! Function to set a bit to 0
 69    subroutine clear_prime_bit(num)
 70        integer(i18t), intent(in) :: num
 71        integer(i18t) :: i
 72        i = (num / bits_per_int) + 1
 73        prime_sieve(i) = iand(prime_sieve(i), ieor(-1_i18t, 2**int(mod(num, bits_per_int), kind(prime_sieve(1)))))
 74    end subroutine clear_prime_bit
 75
 76    ! Function to check if a bit is set
 77    pure logical function get_prime_bit(num) result(bit)
 78        integer(i18t), intent(in) :: num
 79        bit = logical(btest(prime_sieve(num / bits_per_int + 1), int(mod(num, bits_per_int), kind(prime_sieve(1)))))
 80    end function get_prime_bit
 81
 82    subroutine prime_factor(num, factor)
 83        integer(i18t), intent(inout) :: num
 84        integer(i18t), intent(out) :: factor
 85        if (num < 0) then
 86            factor = -1
 87            num = -num
 88        elseif (num <= 1) then
 89            factor = num
 90        else
 91            factor = 2
 92            call expand_sieve(int(sqrt(real(num)), i18t) + 1)
 93            do while (factor <= num)
 94                if (mod(num, factor) == 0) then
 95                    num = num / factor
 96                    return
 97                end if
 98                factor = next_prime(factor)
 99            end do
100        endif
101    end subroutine
102
103    integer(i18t) function is_composite(num) result(factor)
104        integer(i18t), intent(in) :: num
105        integer(i18t) :: localnum
106        localnum = num
107        call prime_factor(localnum, factor)
108        if (factor == num) then
109            factor = 0
110        endif
111    end function
112
113    logical function is_prime(num) result(ip)
114        integer(i18t), intent(in) :: num
115        integer(i18t) check
116        if (num < 2) then
117            ip = .false.
118        elseif (num <= 3) then
119            ip = .true.
120        else
121            check = mod(num, 6_i18t)
122            ip = (check == 1 .or. check == 5) .and. is_composite(num) == 0
123        endif
124    end function
125end module primes

Tags: fortran-iterator, prime-number