Skip to content

Instantly share code, notes, and snippets.

@alecloudenback
Last active December 28, 2022 03:33
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save alecloudenback/5c3a4d2511cfbaa158c9556e6871ac70 to your computer and use it in GitHub Desktop.
Save alecloudenback/5c3a4d2511cfbaa158c9556e6871ac70 to your computer and use it in GitHub Desktop.
List of primes

Primes Benchmark

I got nerd-sniped from this LinkedIn post. On my computer, the R screenshot above ran in 1.564 minutes on my laptop for primes up to 1 million (Macbook Air M1 with arm64 binaries for Julia/R).

I had to give it a shot for comparison's sake. A basic 1:1 translation in Julia was 13.43 seconds (7x faster), and following the same conceptual algorithm (check for remainders up to sqrt(x)) but without as many intermediate allocations was 0.017 seconds (over 5,250x faster).

Prime divisibility check algorithim

R - 93.8 seconds

library(magrittr)
library(purrr)
library(dplyr)
library(lubridate)

primes_less_than <- 1000000

v.primes <- c(2,3)

start <- now()
walk((2:ceiling(primes_less_than/2)) * 2 -1,
  function(x) {
    if ((map_dbl(v.primes[v.primes <= sqrt(x)],
                 ~ x %% .) %>% min) > 0) {
                   v.primes[length(v.primes) + 1] <<- x
    }
    }) %>%
  suppressWarnings

print(now() - start)

Julia - 13.43 sec (7x faster)

function primes_up_to(k) 
	primes = [2,3]
	for n in 3:2:k 
		remainders = n .% primes[primes .<= (n)]
		if  minimum(remainders;init=1) > 0 
			push!(primes,n)
		end
	end
	return primes
end

Julia - 0.029 seconds (3,250x faster)

This follows the same idea (check for remainders up to sqrt(k)) but does not have as many intermediate allocations.

function primes_up_to(k)
	primes = [2,3]
	for n in 3:2:k 
		d = false
		i = 1
		while (~d) && (i < length(primes)) && (primes[i] <= sqrt(n)) 
			
			if iszero(n % primes[i])
				d = true
			end
			i += 1
		end
		if ~d
			push!(primes,n)
		end
			
			
	end
	return primes
end

Julia - 0.0175 seconds (5,500x faster)

By splitting into sub-functions, it becomes more readable (IMO) and facilitates some further optimizations.

function primes_up_to(k)
	primes = [2,3]
	for n in 3:2:k 
		if ~is_divisible(primes,n)
			push!(primes,n)
		end
	end
	return primes
end

function is_divisible(prime_list,n)
	# integer sqrt keeps everything integers instead of comparing float with int
	limit = isqrt(n) 
	
	for i  eachindex(prime_list)
		p = prime_list[i]
		
		# short circuit and return if the conditionals are met
		p > limit && return false
		iszero(n % p) && return true
					
	end
	
	# we've exhausted the search and can return false
	return false	
end

Other methods

These aren't directly comparable, as they have different algorithms. Two other options:

Primes.jl - 0.0013 seconds

using Primes
primes(1_000_000)

FastPrimeSieve - 0.0003 seconds

using FastPrimeSieve
collect(FastPrimeSieve.SmallSieve(1_000_000))
@bkamins
Copy link

bkamins commented Dec 27, 2022

I did some minor tweaks of your Julia code (without changing the algorithm used) and wonder how it compares:

function primes_up_to(k)
	primes = UInt32[2, 3]
	for n in UInt32.(3:2:k) 
		is_divisible(primes, n) || push!(primes, n)
	end
	return Int.(primes) # still return vector of Int as originally requested
end

function is_divisible(prime_list, n)
	limit = isqrt(n) 
	
	for p ∈ prime_list
		p > limit && return false
		iszero(n % p) && return true
	end
	return false	
end

(I hardcoded UInt32, but alternatively one could compare k to typemax(UInt32) and choose wider type if k is large).

@alecloudenback
Copy link
Author

@bkamins your version is indeed faster:

Main solution:

julia> @benchmark primes_up_to(1_000_000)
BenchmarkTools.Trial: 321 samples with 1 evaluation.
 Range (min … max):  15.367 ms …  21.083 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     15.521 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   15.582 ms ± 394.369 μs  ┊ GC (mean ± σ):  0.05% ± 0.56%

    ▁▅▇██▆▆▄▃▂▁
  ▄████████████▆▇▆▄▆▆▄▄▄▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▄▄▁▁▁▁▄ ▇
  15.4 ms       Histogram: log(frequency) by time      16.8 ms <

 Memory estimate: 1.83 MiB, allocs estimate: 11.

UInt32 version:

julia> @benchmark primes_up_to(1_000_000)
BenchmarkTools.Trial: 361 samples with 1 evaluation.
 Range (min … max):  13.486 ms …  15.491 ms  ┊ GC (min … max): 0.00% … 11.32%
 Time  (median):     13.738 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.851 ms ± 364.746 μs  ┊ GC (mean ± σ):  0.65% ±  2.20%

      ▁▄▇▇█▆▅▃▁▁
  ▄▄▇▇██████████▇▆▄▁▄▆▆▆▁▄▄▁▄▁▁▁▆▁▁▄▁▁▄▄▄▄▆▆▁▄▁▄▆▁▄▇▁▆▄▄▄▁▄▆▄▆ ▇
  13.5 ms       Histogram: log(frequency) by time      15.2 ms <

 Memory estimate: 3.42 MiB, allocs estimate: 15.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment