On the search of a prime sieve

Stream: lean4

[Bernardo Borges]

Hello everybody! I am tackling a programming problem that involves the Sieve of Eratosthenes, meaning I want to precalculate the prime numbers up to 2000000 (2M). The first naive approach was the pure functional one:

def fromTo (l r : Nat) : List Nat := map (·+l) $ range (r-l+1)

partial def primes : List Nat := go (fromTo 2 2000000)
  where
    go : List Nat → List Nat
    | [] => []
    | p :: xs => p :: go (xs.filter  (· % p != 0))

No, it would'nt even run.
To make it more feasible, i diverged to a imperative approach. Now I will try to generate a is_prime : Array Bool that tells me what I need to know. I arrived at this code, based on a Rust submission that worked:

/-- Range with `[l,r)` step `s` -/
def fromToStep (l r s : Nat) : List Nat := map (·+l) $ map (·*s) $ range ((r-l+1)/s)

def LIM := 2000000

def sieve : IO (Array Bool) := do
  let mut p := Array.mk $ replicate (LIM+1) true
  p := p.set! 0 false
  p := p.set! 1 false
  for i in fromToStep 4 (LIM+1) 2 do
    p := p.set! i false
  let mut i := 3
  while i * i <= LIM + 1 do
    if p.get! i then
      for j in fromToStep (i*i) (LIM+1) (2*i) do
        p := p.set! j false
    i := i + 1
  return p

This attempt gave me a

Stack overflow detected. Aborting.

What can be done next? Can I use the heap? Happy to hear any pointers!

[Jireh Loreaux]

Well, I can't help but say that if you're implementing the Sieve of Erastosthenes, then you should really have a look at this paper: https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf

[Jireh Loreaux]

Barring that I have nothing meaningful to say .

[Number Eighteen]

My best attempt at implementing naive Eratosthenes:

@[specialize, inline]
partial def iterfun {α : Type u} (x y : Nat) (a : α) (f : α → α) : α :=
  cond (y ≤ x) a (iterfun x (y-1) (f a) f)

partial def simpleSieve (n : Nat) : Array Nat :=
  let s := n.toFloat.sqrt.ceil.toUInt64.toNat
  let rec loop i A :=
    if i ≥ s then
      A.foldl (fun (R, j) b ↦ (cond b (R.push j) R, j.succ)) (#[], 0) |>.fst
    else
      loop i.succ <| cond A[i]!
        (iterfun i (n/i).succ (A, i)
          (fun (R, j) ↦ (R.set! (j * i) false, j.succ)) |>.fst)
        A
  loop 2 <| Array.mkArray n.succ true

This can sift up to a couple of billion integers in 30sec/1min. All optimizations I tried on this code made it significantly slower, which is very frustrating.