Sieve of Atkin implementation on Elixir

Why an Atkin?

On Animus, we have the tradition of publishing a weekly challenge that keeps us sharp in logical thinking, and treating it like a game. One week, someone throws out the challenge, “Find prime numbers within a given value”. Cool!

Joe, a great Animus Rust programmer, presented a solution with the Sieve of Eratosthenes and I was impressed. So my nerd-curiosity was triggered and I started to look for another sieve, founding the Atkin solution.

A bit of context

The Sieve of Atkin's paper was published in 2003 by A. O. L. Atkin and Daniel J. Bernstein. It filters composite numbers similar to Eratosthenes, but this algorithm operates in terms of modulo-60 remainders. Additionally, it has a better theoretical complexity than others:

Sieve of Eratosthenes: O(N*log(log(N)))
Sieve of Sundaram: O(N*log N)
Sieve of Atkin: (N / (log (log N)))

Furthermore, the Sieve of Atkin is a better solution than Eratosthenes’s because It performs some preliminary work and then marks off multiples of squares of primes. This itself represents an improvement in thetime complexity.

The rules

1. 2 and 3 are special cases. So those are added to the primes list without processing.
2. 4x² + y² = n is odd and modulo-12 remainder is 1 or 5
3. 3x² + y² = n is odd and modulo-6 remainder is 1
4. 3x² — y² = n is odd and modulo-12 remainder is 11

If any number in our list is square multiple, then remove it.

Let’s code it!

• Convert the given number into a range and then into a list so that we can iterate over it.
• Prepare a map to associate true/false with the integers inside.
  *To avoid O(n) complexity, we need to use maps instead of keyword lists.
• Create helper functions to detect and add 2 and/or 3 to the sieve.

Check if 2 or 3 are within the range and add them.
*Note that we are delegating some logic to ‘update_sieve/2’, and it is crucial to understand that this function reverses the current status of the prime. For example, certain conditions might set a non-prime number like 25 as true, only to be modified later by another condition which might set it back to false.

defmodule SieveOfAtkin do
  def run(limit) do
    sieve_range = 1..limit |> Enum.to_list()
    sieve = 1..limit |> Enum.into(%{}, fn i -> {i + 1, false} end)

    sieve = maybe_add_2_or_3(limit, sieve)
  end

  defp update_sieve(sieve, n) do
    # Invert the current value
    Map.update(sieve, n, false, &(!&1))
  end
  
  defp maybe_add_2_or_3(limit, sieve) when limit >= 3 do
    update_sieve(sieve, 2)
    |> update_sieve(3)
  end

  defp maybe_add_2_or_3(limit, sieve) when limit == 2, do: update_sieve(sieve, 2)
  defp maybe_add_2_or_3(_limit, sieve), do: sieve
end

We have now to iterate x and y until x * x > limit or y * y > limit. Let's see this in Python first to fully understand the idea.

while x * x <= limit
  y = 1
  while y * y <= limit
    ...
    y += 1
  x += 1

Elixir version:

defmodule SieveOfAtkin do  
  def run(limit) do
    sieve_range = 1..limit |> Enum.to_list()
    sieve = 1..limit |> Enum.into(%{}, fn i -> {i + 1, false} end)
  
    sieve = maybe_add_2_or_3(limit, sieve)
  
    iterate_x(sieve_range, limit, sieve, sieve_range)
  end
  
  defp iterate_x([x | _tail], limit, sieve, _sieve_range) when x * x >= limit, do: sieve
  defp iterate_x([], _limit, sieve, _sieve_range), do: sieve
  
  defp iterate_x([x | tail], limit, sieve, sieve_range) do
    iterate_y(sieve, x, tail, sieve_range, limit, sieve_range)
  end
  
  defp iterate_y(sieve, _x, tail_x, [], limit, sieve_range),
    do: iterate_x(tail_x, limit, sieve, sieve_range)
  
  defp iterate_y(sieve, _x, tail_x, [y | _tail], limit, sieve_range) when y * y >= limit,
    do: iterate_x(tail_x, limit, sieve, sieve_range)
  
  defp iterate_y(sieve, x, tail_x, [y | tail], limit, sieve_range) do
    # Check conditions 1, 2 and 3
  end
  
  
  defp maybe_add_2_or_3(limit, sieve) when limit >= 3 do
    update_sieve(sieve, 2)
    |> update_sieve(3)
  end
  
  defp maybe_add_2_or_3(limit, sieve) when limit == 2, do: update_sieve(sieve, 2)
  defp maybe_add_2_or_3(_limit, sieve), do: sieve
end

The conditions:

defp condition1(sieve, x, y, limit) do
  # 4x² + y² = n
  n = 4 * x * x + y * y

  if(n <= limit and rem(n, 12) in [1, 5]) do
    update_sieve(sieve, n)
  else
    sieve
  end
end

defp condition2(sieve, x, y, limit) do
  # 3x² + y²
  n = 3 * x * x + y * y

  if(n <= limit and rem(n, 12) == 7) do
    if x == 5 and y == 10 do
    end

    update_sieve(sieve, n)
  else
    sieve
  end
end

defp condition3(sieve, x, y, limit) do
  # 3x² - y²
  n = 3 * x * x - y * y

  if(x > y and n <= limit and rem(n, 12) == 11) do
    update_sieve(sieve, n)
  else
    sieve
  end
end

Perfect Squares

The last condition is to check that the numbers we have added to the primes list are Squares-Free, which means prime numbers can’t be multiples of squares. Again, I will show this logic in Python first to fully understand the idea, and then the implementation in Elixir.\

Python:

while r * r <= limit:
  if sieve[r]:
      for i in range(r * r, limit+1, r * r):
          sieve[i] = False
  r += 1

Elixir:

defp is_perfect_square(sieve, r, limit) when r * r <= limit do
  case sieve do
    %{^r => true} ->
      sieve =
        Stream.unfold(r * r, fn x ->
          if x <= limit, do: {x, x + r}
        end)
        |> Enum.reduce(sieve, fn i, sieve ->
          Map.put(sieve, i, false)
        end)

      is_perfect_square(sieve, r + 1, limit)

    _ ->
      is_perfect_square(sieve, r + 1, limit)
  end
end

Run the conditions in order:

def run(limit) do
  sieve_range = 1..limit |> Enum.to_list()
  sieve = 1..limit |> Enum.into(%{}, fn i -> {i + 1, false} end)

  sieve = maybe_add_2_or_3(limit, sieve)

  iterate_x(sieve_range, limit, sieve, sieve_range)
  |> is_perfect_square(5, limit)
  |> Enum.filter(fn {_i, value} -> value end)
  |> Enum.map(fn {i, _value} -> i end)
  |> Enum.sort()
end

Complete code:

Final thoughts:

What I’ve learned from this exercise is that:

• This code can (and will) be improved. The use of the Enum module is not the best practice for performance and can be replaced with recursive functions.
• The use of “if” statements should be replaced as a best practice.
• Currently, this script takes one second to process how many prime numbers may be found in one million while the Rust version takes only 0.3 seconds. We’re not attempting to match Rust’s power here; Elixir isn’t designed for that purpose. However, I believe that the processing time can be further improved.

iex(1)> SieveOfAtkin.run(1_000_000)
Executed in 1.310634s
Restults: 78498

Credits:

Author: Julián Somoza

Editor: Laura Arcuri

https://www.geeksforgeeks.org/sieve-of-atkin/