eratosthenes sieve

Stream: general

[Asei Inoue]

I implemented the sieve of Eratosthenes with Lean. It has only been implemented and the proof has not yet been written. I want to prove that this function really returns a list of prime numbers less than or equal to n, but how do I prove this for a function defined in do notation?

def eratosthenes (n : Nat) : Array Bool := Id.run do
  let mut isPrime := Array.mkArray (n + 1) true

  isPrime := isPrime.set! 0 false
  isPrime := isPrime.set! 1 false

  for p in [2 : n + 1] do
    if not isPrime[p]! then
      continue

    let mut q := 2 * p
    while q <= n do
      isPrime := isPrime.set! q false
      q := q + p

  return isPrime

def runTest (n : Nat) : IO Unit := do
  let l := eratosthenes n
  for i in [:l.size] do
    if l[i]! then
      IO.print s!"{i} "

#eval runTest 19

[Asei Inoue]

my attempt is here:

theorem prime_iff_eratos (n p : Nat) (hp : p ≤ n) : (eratosthenes n)[p]! ↔ Nat.Prime p := by
  constructor <;> intro h

  case mp =>
    rw [@Nat.prime_def_lt]

    constructor
    case left =>
      sorry

    case right =>
      sorry
  sorry

[Mauricio Collares]

Related: https://leanprover.zulipchat.com/#narrow/stream/270676-lean4/topic/Checking.20the.20Goldbach.20conjecture/near/429974821

[Harald Helfgott]

Hi - just two obvious remarks: it is enough for p to go up to sqrt(n), and we can start wtrh q=p*p. The former point will be important when we go to the segmented sieve.

[Asei Inoue]

Thanks. How about this?

def eratosthenes (n : Nat) : Array Bool := Id.run do
  let mut isPrime := Array.mkArray (n + 1) true

  isPrime := isPrime.set! 0 false
  isPrime := isPrime.set! 1 false

  for p in [2 : n + 1] do
    if not isPrime[p]! then
      continue

    if p ^ 2 > n then
      break

    let mut q := p * p
    while q <= n do
      isPrime := isPrime.set! q false
      q := q + p

  return isPrime

def runTest (n : Nat) : IO Unit := do
  let l := eratosthenes n
  for i in [:l.size] do
    if l[i]! then
      IO.print s!"{i} "

#eval runTest 3190

[Harald Helfgott]

Also, I am a beginner - how do we check that the Lean compiler has realized that we are referring to eacha changed array entry at most once (within the function), so that the entire array isn't copied?

[Harald Helfgott]

The code looks good. Let me run it.

[Asei Inoue]

how do we check that the Lean compiler has realized that we are referring to eacha changed array entry at most once

I wonder if we could also make implementations that are obviously bad and compare their execution times by #time.

[Henrik Böving]

Harald Helfgott said:

Also, I am a beginner - how do we check that the Lean compiler has realized that we are referring to eacha changed array entry at most once (within the function), so that the entire array isn't copied?

You cannot reasily and reliably do that right now. The most accurate thing I can offer is to read the IR which is hard for newcomers obviously, the next best thing is dbgTraceIfShared but it is possible that the compiler will optimize your code differently in the presence of it.

There are a number of ideas floating around on how one might address this like type systems or opt-in runtime checks but all of those have issues and we do not currently have time to explore them.

[Asei Inoue]

:sad: oh

[Jason Rute]

Harald Helfgott said:

Also, I am a beginner - how do we check that the Lean compiler has realized that we are referring to eacha changed array entry at most once (within the function), so that the entire array isn't copied?

See this discussion for how to tell if an object is being duplicated: https://leanprover.zulipchat.com/#narrow/stream/270676-lean4/topic/Efficient.20hashmap.20counter.20from.20list

[Jason Rute]

(But Hendrik likely knows better than me, so take it with his warnings in mind.)

[Henrik Böving]

Yes that id dbgTraceIfShared but it can occasionally give you false positives or negatives. But in general what marc describes there is the closest approximation we have beyond "read the IR"

[Asei Inoue]

it outputs same addresses. it would be really the eratosthenes sieve.

unsafe def eratosthenes (n : Nat) : Array Bool := Id.run do
  let mut isPrime := Array.mkArray (n + 1) true

  isPrime := isPrime.set! 0 false
  isPrime := isPrime.set! 1 false

  for p in [2 : n + 1] do
    dbg_trace ptrAddrUnsafe isPrime
    if not isPrime[p]! then
      continue

    if p ^ 2 > n then
      break

    let mut q := p * p
    while q <= n do
      isPrime := isPrime.set! q false
      q := q + p

  return isPrime

#eval eratosthenes 131