How do I prove my prime sieve works?

Stream: new members

[Michal Wallace (tangentstorm)]

image.png

Hello,

I've implemented a prime sieve based on juggling arithmetic sequences around, and I see that it works. But now I would like to prove it.

https://github.com/tangentstorm/treesiv/blob/main/Main.lean

structure ASer where  -- arithmetic series (k + dn)
  k : Nat  -- constant
  d : Nat  -- difference
deriving Inhabited, Repr

instance : Ord ASer where  -- but there's no List.sort? !!
  compare s1 s2 :=
    match compare s1.k s2.k with
    | .eq => compare s1.d s2.d
    | ord => ord

-- mk k d = k + d*n
def mk (k d : Nat) : ASer := { k := k, d := d }

-- apply formula to n
def ap (s : ASer) (n : Nat) : Nat := s.k + s.d * n

-- apply one formula to another: r(n) := s(t(n))
def compose (s : ASer) (t : ASer) : ASer := mk (s.k + s.d * t.k) (s.d * t.d)

-- identity series
def id1 : ASer := mk 0 1

def partition (s : ASer) (n : Nat) : List ASer :=
  List.range n |>.map fun i => compose s $ mk i n

instance : ToString ASer where
  toString s :=
    if s.k == 0 then
      if s.d == 1 then "n"
      else s!"{s.d}n"
    else if s.d == 0 then s!"{s.k}"
    else if s.d == 1 then s!"n + {s.k}"
    else s!"{s.d}n + {s.k}"

def terms (s : ASer) (n : Nat) : List Nat := List.range n |>.map λ i =>ap s i

structure PrimeSieve where
  ps : List Nat         -- all primes we've used so far
  pr : Nat              -- current primorial (product of ps)
  np : Nat              -- next prime
  ss : List ASer        -- list of sequences
deriving Repr

instance : ToString PrimeSieve where
  toString s := s!"ps: {s.ps}, pr: {s.pr}, np: {s.np}, ss: {s.ss}"

def init : PrimeSieve := { ps := [], pr := 1, np := 2, ss := [id1] }

def step (s0 : PrimeSieve) : PrimeSieve :=
  let ps := s0.ps ++ [s0.np]
  let pr := s0.pr * s0.np
  let ss0 := (s0.ss.map fun s => partition s s0.np).join
  let ss := (ss0.filter fun s => s.k % s0.np != 0)  -- strip out multiples of np
  let np := (List.minimum? $ ss.map fun s => (let f1:=ap s 0; if f1 == 1 then ap s 1 else f1)).get! -- series with next prime
  { ps := ps, pr := pr, np := np, ss := ss }

def printStep (s : PrimeSieve) (n : Nat) : IO Unit := do
  IO.println s
  s.ss.forM fun s => do
    let width := 15
    let formula := (String.pushn (toString s) ' ' width).take width
    IO.println s!"{formula}: {(terms s n)}"

def main : IO Unit := do
  let mut sv := init
  let n := 10
  printStep sv n
  for _ in [0:5] do
    sv := step sv
    IO.println "----------------"
    printStep sv n

#eval main

But how do I even start?

Currently, it tracks a list of known primes and a next prime (np). I think I want to force things so that the next prime is of type Prime... Would I just bake a "prime" type directly into my data structure?

I think what I want to do is use the type system to assert:

• the "next prime" is a prime
• the "known primes" are all prime, and also contain ALL primes less than "next prime"

Then I think I to prove that I calculate the next prime correctly, I'd have to demonstrate that my list of arithmetic series together generate all natural numbers that are not multiples of a known prime, and (since they all have the same slope) the next prime is the minimal term in one of the series.

Does this general plan make sense? And if so, how do I even express it?

I'm going to start by just setting the type of "next prime" and see where it leads, but if anyone could give me a few hints, I'd really appreciate it.

Thanks!

[Jireh Loreaux]

We prefer #backticks to unreadable screenshots. Can you update?

[Mauricio Collares]

I think the link to the code matches what's in the screenshot (minus the infoview output)

[Michal Wallace (tangentstorm)]

Yes, but I'll post the code directly. :)

[Jireh Loreaux]

(oops, I'm sorry, I missed this SHOW MORE button this time, my apologies)

[Michal Wallace (tangentstorm)]

It seems like Mathlib.Data.Nat.Prime is a proposition, not just a Type.

I found this:

def NPrime : Type := { n: Nat // Nat.Prime n } deriving Repr
instance : ToString NPrime where
  toString s := s!"{s.val}"

Then I was able to declare my types like this:

structure PrimeSieve where
  ps : List NPrime      -- all primes we've used so far
  pr : Nat              -- current primorial (product of ps)
  np : NPrime           -- next prime
  ss : List ASer        -- list of sequences
deriving Repr

So far that required two changes: every time I assign the np field, I now need to pair it with a proof that the number is prime in special angle brackets.

The proof that the initial 2 is prime is provided by the Mathlib:

def init : PrimeSieve := { ps := [], pr := 1, np := ⟨2,Nat.prime_two⟩, ss := [id1] }

But for step, I'm on my own:

{ ps := ps, pr := pr, np := ⟨np,sorry⟩, ss := ss }

I'm not sure how I can even start to prove this without knowing that ps contains all the primes less than np, which won't be true unless I bake it into the type, and I don't know how to express that yet.

Possibly I can just add any proposition to the type with the magic // symbol?

[Luigi Massacci]

You can add fields to your structure that are propositions

[Michal Wallace (tangentstorm)]

Interesting.. I can almost see how that would work. I think the proof outline is something like this:

• theorem: np really is the next prime after each step, because
  • the items in ss together contain all nats > 2 that are not multiples of a current prime (and i need to prove that this holds at each partition step). this means:
    • the partition operation must actually partition the series (and not accidentally lose a value that might be the next prime)
    • the filter operation removes exactly the multiples of the current next prime
  • the sequences all have the same slope (.ss[i].d = .pr) so the lowest first term is the next prime (though currently I still have a 1 floating around as the first term in one series at each step)

[Michal Wallace (tangentstorm)]

On further reflection, I think what I really ought to do is prove a basic theorem about the idea of prime sieves:

• The smallest natural number that is coprime to all primes smaller than it is prime.

Then for my particular implementation, I need to show that:

• the list of known primes is complete up to the "next prime" variable
• the list of arithmetic progressions together represent every natural number except those coprime to the list of known primes.
• that from these, i am correctly identifying the minimal coprime number

I think the "representing all the coprimes" part is two more or less induction, except not using the "induction" tactic but an invariant on the type itself, so:

• the initial sieve value represents all natural numbers (maybe all > 2)
• at each step we remove exactly the multiples of the next prime (so we use the proof on the input as the induction hypothesis)

So the first part is an obvious-to-me mathematical statement about why prime sieves produce primes (so hopefully easy to prove in lean??)... And then once that framework is in place to describe formally what a prime sieve is, I can prove that this particular algorithm and data structure actually is a prime sieve.

I will continue to post here to document my progress. If anyone has observations/ideas for making any of this more approachable, please comment!

[Michal Wallace (tangentstorm)]

Woo! I managed to prove that if a natural number greater than 1 has no prime factors less than it, it must be prime:

theorem no_prime_factors_im_no_factors {c:ℕ} -- c is a candidate prime
  (hcg2: 2 ≤ c)                              -- c is greater than 2
  (hnpf: ∀ p < c, Nat.Prime p → ¬(p∣c))      -- c has no prime factors
  : Nat.Prime c := by
    have hno: ∀ m < c, m ∣ c → m=1 := by  -- c has no divisors at all
      intros m hmc hmdc -- give names to the above assumptions
      by_contra hnme1   -- assume ¬(m=1) and show contradiction
      have hmn1: m ≠ 1 := by assumption
      -- obtain a prime factor of m since it isn't 1
      obtain ⟨p, hpp, hpm⟩ : ∃ p, Nat.Prime p ∧ p ∣ m := Nat.exists_prime_and_dvd hmn1
      have h0: p∣c := Nat.dvd_trans hpm hmdc
      have h1: ¬(p∣c) := by
        have hzlc: 0 < c := by linarith
        have hzlm: 0 < m := Nat.pos_of_dvd_of_pos hmdc hzlc
        have hplm: p ≤ m := Nat.le_of_dvd hzlm hpm
        have hplc: p < c := by linarith
        exact hnpf p hplc hpp
      contradiction
    exact Nat.prime_def_lt.mpr ⟨hcg2, hno⟩

This was the main step that was worrying me, because it maps the conditions provided by a sieve to the conditions that mathlib uses in the definition of Nat.Prime.

[Bbbbbbbbba]

I guess if you are using Nat.exists_prime_and_dvd you may as well do this?

theorem no_prime_factors_im_no_factors {c:ℕ} -- c is a candidate prime
  (hcg2: 2 ≤ c)                              -- c is greater than 2
  (hnpf: ∀ p < c, Nat.Prime p → ¬(p∣c))      -- c has no prime factors
  : Nat.Prime c := by
    have hcn1 : c ≠ 1 := Ne.symm (Nat.ne_of_lt hcg2)
    obtain ⟨p, hpp, hpc⟩ : ∃ p, Nat.Prime p ∧ p ∣ c := Nat.exists_prime_and_dvd hcn1
    have hpec : p = c := by
      have hzlc: 0 < c := by linarith
      have hplc : p ≤ c := Nat.le_of_dvd hzlc hpc
      have hpgc : p ≥ c := Nat.le_of_not_lt (fun a => hnpf p a hpp hpc)
      exact Nat.le_antisymm hplc hpgc
    rw [← hpec]
    exact hpp

[Michal Wallace (tangentstorm)]

I am quite stuck at the moment. I want to refer to the infimum of the following set:

-- the set of naturals with no prime factors less than some c
-- "remaining set"? "residual set?"
def rs (c : Nat) : Set Nat := { n : Nat | ∀ p < c, Nat.Prime p → ¬(p∣n) }

I want to prove that :

• my algorithm produces the infimum of this set, and
• the infimum of the set is prime

I am pretty sure I can prove both of these statements, but I'm stuck trying to figure out how to express the idea of an infimum... I see that Inf $ rs c is an expression that produces a type, but I don't understand what to do with it.

How would I obtain a member of that type? What would I do with it? I think I'm missing something basic here.

[Michal Wallace (tangentstorm)]

Actually I might even be using the term wrong. I want the minimal element of the set. I know it exists. I can obtain it. I just don't know how to express the concept so I can prove I'm doing it. Any ideas?

[Yakov Pechersky]

You want docs#sInf, s for set

[Yakov Pechersky]

That will give you a Nat. Plain Inf is the binary operator of two elements

[Yakov Pechersky]

Regarding your definition, wouldn't 1 be in the set?

[Kevin Buzzard]

If c<=2 then the min is zero and if c>2 it's 1, so in fact the minimum is never prime. Is there a typo somewhere?

[Yakov Pechersky]

Perhaps adding the predicate that c < n to the set will rescue the definition?

[Daniel Weber]

It might be easier to use docs#Nat.find ?

[Mario Carneiro]

you realize you are reinventing docs#Nat.minFac ?

[Mario Carneiro]

no_prime_factors_im_no_factors also looks like a theorem that should already exist

[Michal Wallace (tangentstorm)]

I don't understand what you guys are saying. I'm trying to express the idea that rs(c) is the set of all natural numbers (n) such that n is not a multiple of a prime less than c. It's a prime sieve, so I'm trying to describe all the numbers that haven't been crossed out yet, and then I want to obtain the smallest number in that set so I can demonstrate that it is prime.

[Michal Wallace (tangentstorm)]

Is rs saying something other than that?

[Michal Wallace (tangentstorm)]

this already works:

-- if c is a member of rs c, then c is prime
lemma cprime (c : Nat) (h2: c≥2) (hcrc: c ∈ rs c ) : Nat.Prime c :=  by
  have hh : ∀ p < c, Nat.Prime p → ¬(p∣c) := by
    simp[rs] at hcrc
    exact hcrc
  exact no_prime_factors_im_no_factors h2 hh

-- show that 2 is part of rs 2 for initial case
lemma cprime2 : 2 ∈ rs 2 := by
  simp [rs]
  intro p hp hprime
  have h1 : p ≤ 1 := Nat.le_of_lt_succ hp
  have h2 : p ≥ 2 := (Nat.prime_def_lt.mp hprime).left
  have h3 : ¬(p ≤ 1) := not_le_of_gt h2
  contradiction

[Michal Wallace (tangentstorm)]

Oh ... I see my error now. Yes, I need to add c ≤ n... I think in this case the infimum is also the minimum but I'm not sure. I'll try to proceed with sInf. Thanks!

[Mario Carneiro]

do you have a #mwe of this?

[Michal Wallace (tangentstorm)]

Sort of. Here's what I have so far, and all of this typechecks.

import Mathlib.Data.Nat.Prime
import Mathlib.Tactic.Linarith.Frontend

lemma no_prime_factors_im_no_factors {c:ℕ} -- c is a candidate prime
  (h2lc: 2 ≤ c)                              -- c is at least 2
  (hnpf: ∀ p < c, Nat.Prime p → ¬(p∣c))      -- c has no prime factors
  : Nat.Prime c := by
    have : ∀ m < c, m ∣ c → m=1 := by  -- c has no divisors but 1 and c
      intro m hmc hmdc  -- give names to the above assumptions
      by_contra         -- assume ¬(m=1) and show contradiction
      obtain ⟨p, hpp, hpm⟩ : ∃ p, Nat.Prime p ∧ p ∣ m :=
        Nat.exists_prime_and_dvd ‹m ≠ 1›
      have : p∣c := Nat.dvd_trans hpm hmdc
      have : ¬(p∣c) := by
        have : 0 < c := by linarith
        have : 0 < m := Nat.pos_of_dvd_of_pos hmdc this
        have : p ≤ m := Nat.le_of_dvd this hpm
        have : p < c := by linarith
        exact hnpf p this hpp
      contradiction
    exact Nat.prime_def_lt.mpr ⟨‹2 ≤ c› , this⟩

def rs (c : Nat) : Set Nat := { n : Nat | c ≤ n ∧ ∀ p < c, Nat.Prime p → ¬(p∣n) }

lemma cprime (c : Nat) (h2: c≥2) (hcrc: c ∈ rs c ) : Nat.Prime c :=  by
  have hh : ∀ p < c, Nat.Prime p → ¬(p∣c) := by
    simp[rs] at hcrc
    exact hcrc
  exact no_prime_factors_im_no_factors h2 hh

-- show that 2 is part of rs 2 for initial case
lemma cprime2 : 2 ∈ rs 2 := by
  simp [rs]
  intro p hp hprime
  have h1 : p ≤ 1 := Nat.le_of_lt_succ hp
  have h2 : p ≥ 2 := (Nat.prime_def_lt.mp hprime).left
  have h3 : ¬(p ≤ 1) := not_le_of_gt h2
  contradiction

But I don't have anything about the minimal item, because I don't know how to express it.

[Mario Carneiro]

what happened to the sieve program?

[Michal Wallace (tangentstorm)]

Oh, the sieve is here: https://github.com/tangentstorm/treesiv/blob/main/Main.lean

[Mario Carneiro]

What is it you are trying to prove?

[Mario Carneiro]

I think what you want to prove by induction or execution is that minFac n >= a

[Mario Carneiro]

and if a gets to n then you know n is prime

[Michal Wallace (tangentstorm)]

Hang on. I'm trying to put the answer to your first question into words.

[Mario Carneiro]

-- but there's no List.sort? !!

There is List.mergeSort

[Mario Carneiro]

What is the purpose of this program? It seems like it does exponential work compared to the actual list of primes coming out in ps

[Michal Wallace (tangentstorm)]

A prime sieve (at least one that works like the sieve of eratosthenes) is some algorithm that at each step, has identified all the primes up to some prime P, and has "scratched out" all the multiples of those primes from the set of all natural numbers.

I have implemented a prime sieve, and want to prove that the algorithm works: that it always finds the next prime.

But I guess there are many ways you could implement a sieve like this. For example, the traditional implementation generates an array of booleans corresponding to a range of numbers, and then loops through and flips the bit for each multiple... The algorithm here is not as efficient. It's taking an initial arithmetic sequence (2 + n) and at each step is partitioning it into sub-sequences, using this to find the next prime, and removing the sequences that are multiples of a prime.

[Mario Carneiro]

The first step in proving an algorithm like this correct is to define the invariants. What are all the properties you expect to be true of the input and result of the step function?

[Michal Wallace (tangentstorm)]

Yes. There's no real practical point to it, other than as an exercise. Although it does have a nice property that it very quickly generates a bunch of huge candidate numbers that are somewhat likely to be prime.

[Mario Carneiro]

Some of them are in comments, and you can move them to fields of the structure like

pr : Nat              -- current primorial (product of ps)

as

pr_eq : pr = ps.foldr (·.1 * ·) 1

[Mario Carneiro]

Oh I wasn't questioning the practicality, only trying to understand what properties you are going for

[Michal Wallace (tangentstorm)]

From a doc I'm writing:

The sieve tracks this info:

• C is the “current prime”
• S is the set of all primes less than C.
• R (for “remaining”) is the set of naturals with no prime factors in S.
• Q is the set of sequences that model R

The smallest of these is the next prime. (May still need to prove this.)

I suppose I would need to show the following:

• At each step, I have a correct mapping to R.
• R’ := R - multiples of C
• C’ is the smallest one.

How do I prove that I’m tracking R correctly?
Again by “induction” using a proposition field.

-The field says I have a mapping to all naturals that are not multiples of found primes.

• I’m now removing multiples of one new prime.

• A sequence is of the form ( k + dn )
  - k∣d this means it’s k + (jk)n or  k(1+jn)]

• Same applies to any ik + (jk)n

• Every item in these series will be a multiple of K.
• How to show that only sequences whose constant terms are multiples of k yield multiples of k?
• K∣(X + JKN) iff X%K=0, because (JKN+X)%K reduces to X%K
• D is of the form JK, since D is the primorial whose largest factor is K
• So this is how we justify the step of partitioning the previous list of sequences by K: the goal is to make the difference coefficient a multiple of K. The way to do that is to split each sequence into K sequences, or rather to apply each sequence to the sequences  0+Kn, 1+Kn, 2+Kn… (K-1)+Kn

How do I prove that my algorithm finds the smallest item in R?

• Q partitions R into sequences of the form K+Dn
• The smallest item in R must be the smallest constant coefficient in one of these

[Mario Carneiro]

I guess you renamed everything in the lean code? I don't see C, S, R, Q

[Michal Wallace (tangentstorm)]

C S R Q are like an abstract description of what any prime sieve along the lines of eratosthenes must be doing... But the actual implementation can vary.

[Mario Carneiro]

well I'm asking more specifically what properties hold of the object you have called PrimeSieve

[Michal Wallace (tangentstorm)]

First, if it's not clear: I'm a software developer, not in any way a mathematician.... So I'm thinking like a software developer, and that might be a problem here. In my mind, I don't want to just directly prove that this prime sieve works... I want to make a general description of what a generic prime sieve does, and then demonstrate that this particular one does all those things.

So, I'm attempting to write an abstract proof in terms of these C S R Q things, and then I want to demonstrate that at each step, the prime sieve is somehow modelling these things.

[Mario Carneiro]

I think you end up with this list of invariants:

structure PrimeSieve where
  ps : List NPrime      -- all primes we've used so far
  pr : Nat              -- current primorial (product of ps)
  np : NPrime           -- next prime
  ss : List ASeq        -- list of sequences
  pr_eq : pr = ps.foldr (·.1 * ·) 1 -- pr is the product of ps
  ps_sorted : (ps ++ [np]).Pairwise (·.1 < ·.1) -- ps ++ [np] is sorted
  all_primes : ∀ x : NPrime, x.1 < np.1 → x ∈ ps -- ps ++ [np] contains the complete list of primes
  divisible : ∀ n : Nat, n > np.1 → (∀ x ∈ (ps ++ [np]), ¬x.1 ∣ 1) ↔ ∃ s ∈ ss, ∃ i, ap s i = n
    -- a number larger than np is contained in one of the sequences iff it does not divide all listed primes

[Mario Carneiro]

If you want to describe a "generic prime sieve", then write that down as a definition and prove that this prime sieve relates to that

[Michal Wallace (tangentstorm)]

Like, I'm also looking at another sieve someone else implemented: https://github.com/ykonstant1/esiv
And I want to show that that, too, has the behaviors I want in terms of C/S/R/Q.

But none of these algorithms actually have the "set of all natural numbers" lying around... You'd need infinite RAM. But instead, you can show that the algorithm somehow has a mapping to these numbers (a traditional sieve just generates it a bit at a time, this one models it as the list of sequences...)

[Mario Carneiro]

no worries, lean is good at holding "the set of all natural numbers" in its head :)

[Michal Wallace (tangentstorm)]

Yes, that's what I want to do. But I basically need to say "C' is the min item in R" and I don't know how to write that sentence.

[Mario Carneiro]

the simple way is to just write that out as a predicate: C' ∈ R ∧ ∀ x ∈ R, C' ≤ x

[Michal Wallace (tangentstorm)]

Let me try to actually write AbstractPrimeSieve in terms of these numbers.

[Michal Wallace (tangentstorm)]

... that... looks... pretty straightforward... :D

[Michal Wallace (tangentstorm)]

Hrm. Okay. I will try this. Thanks, @Mario Carneiro !!

[Michal Wallace (tangentstorm)]

It seems like the "generic PrimeSieve" ought to be defined as like class PrimeSieve and then I should rename the current thing and show instance : PrimeSieve MyPrimeSieve.

[Mario Carneiro]

That depends, the way you described a prime sieve was more like a concrete prime sieve but for the fact that it deals in noncomputable objects, and you would instead have some representation relation between your prime sieve and the abstract one

[Mario Carneiro]

this is basically the same situation as between Set A and Finset A: the former is an abstract mathematical thing and the latter is a concrete computable representation which can be interpreted to a Set

[Mario Carneiro]

It is possible to describe PrimeSieve as a class but in that case it's still not very abstract, more like parametric over representations

[Yakov Pechersky]

Michal Wallace (tangentstorm) said:

It seems like the "generic PrimeSieve" ought to be defined as like class PrimeSieve and then I should rename the current thing and show instance : PrimeSieve MyPrimeSieve.

Note, in lean, class and instance don't mean what they mean in Python, for example. Instead, class in Python is equivalent to structure in lean. And one creates an instance of a structure by a def foo : MyStruct := ... declaration.

class in lean defined a structure that _additionally_, for the parameters provided, should be found automatically by the typeclass search system, and as such, there should morally be only one instance for every MyStruct X.

[Michal Wallace (tangentstorm)]

So... I wanted to define actual sets, but whenever I refer to membership in the sets, I get errors about being unable to synthesize a "Membership" instance, and I don't know how to address that. It would have made the code look kind of nice:

class NonWorkingPrimeSieve (α : Type u) where
  C : α → NPrime
  S : (x:α) → { n: NPrime | n.val < (C x).val } -- primes < C
  R : (x:α) → { n: Nat | ∀ p, p ∈ (S x) → ¬(p.val ∣ n) }
-- ^ failed to synthesize  Membership ?m.512 ↑{n | ↑n < ↑(C x)}

[Michal Wallace (tangentstorm)]

So here is what I came up with. So far, it's just a definition, but it typechecks.

-- a generic prime sieve
class PrimeSieve (α : Type u) where
  C : α → NPrime
  -- nltC says a naturals less than C
  nltC (x:α) (n:Nat) : Prop := n < (C x).val
  -- primes less than C (the "known primes")
  pltC (x:α) (p:Nat) : Prop := Nat.Prime p ∧ p < (C x).val
  -- natural number coprime to all the known primes
  cpks (x:α) (n:Nat) : Prop := n ≥ 2 ∧ ∀ p, pltC x p → ¬(p∣ n)
  -- this is the invariant we need to have a sieve.
  -- if I can show that a particular instance produces a C
  -- that is the smallest natural that is coprime
  -- to all the known primes, then I should be able to
  -- prove the main result (hSiv) below.
  hMin : (x:α) →
    cpks x (C x).val                 -- C is in R
    ∧ ∀ n, cpks x n → n ≥ (C x).val  -- C is min of R
  init : α
  next : α → α
  -- the key property to prove:
  -- we generate the next consecutive prime at each step
  hSiv (x:α) :
    (x':α := next x) →        -- "after each step"
      (C x).val < (C x').val  -- "we have a new, bigger prime"
      ∧ ¬∃ p, Nat.Prime p →   -- "and there is no prime between them"
        ((C x).val<p ∧ p<(C x').val)
    := sorry
      -- I *think* I can prove this for all sieves,
      -- provided the sieve gives me the hMin fact

[Mario Carneiro]

When you write S : (x:α) → { n: NPrime | n.val < (C x).val }, you aren't declaring a set, you are declaring an element of a set. { n: NPrime | n.val < (C x).val } is one particular set and S is a function producing elements of that set

[Michal Wallace (tangentstorm)]

I am somewhat worried that the predicates here are like "default definitions" and that it's possible to override them with any proposition. I want a keyword like "final" in java.

[Mario Carneiro]

as a result, writing p ∈ (S x) later doesn't make sense because S x isn't a set

[Mario Carneiro]

I am somewhat worried that the predicates here are like "default definitions" and that it's possible to override them with any proposition. I want a keyword like "final" in java.

That's correct, this is not what you want. Try defining functions like nltC before the class, taking C as an argument

[Michal Wallace (tangentstorm)]

... ah, so the set comprehension is also a type signature. how would I declare a function that returns a set then?

[Mario Carneiro]

S : A -> Set NPrime or something

[Mario Carneiro]

and then you can put constraints on the set in later arguments

[Michal Wallace (tangentstorm)]

Hrm. Okay, I see how to do it as a separate proposition...

def pltC (c:NPrime) (p:NPrime) : Prop := p.val < c.val

class PrimeSieve (α : Type u) where
  C : α → NPrime
  S : α → Set NPrime
  hS: ∀ p, p ∈ S x → pltC (C x) p

When you say "later argument", are you saying I can define S and hS together in one definition?

[Mario Carneiro]

no, I was saying you can write hS

[Michal Wallace (tangentstorm)]

got it. thanks! this is very helpful!

[Mario Carneiro]

you can bundle some information inside S as well, but you will have to change the uses slightly since it won't be a plain Set NPrime anymore

[Mario Carneiro]

for example for that hS example you could also say

S : (x : α) -> { s : Set NPrime // ∀ p, p ∈ s → pltC (C x) p}

but then any later statements using S x should use (S x).1 instead

[Mario Carneiro]

or alternatively

S : (x : α) -> Set { p : NPrime // pltC (C x) p}

but now it's not a set of NPrime so you will have some proof obligations or so when talking about NPrimes being a member of the set

[Mario Carneiro]

I would usually opt for the version you wrote, with a plain set and then constraints on it, so that usage code isn't mucked up with casting things into the correct subtypes

[Michal Wallace (tangentstorm)]

It seems like lean will completely ignore a proof in a class definition, presumably until you instantiate that class.
is there a way to force it to check the proof without an instance of the class?

class PrimeSieve (α : Type u) where
  C : α → NPrime
  hMyAssertion (x:α) : (C x).val > 234234 := by
    intro you_can_say_anything_here
    rfl

[Mario Carneiro]

The syntax := by in class declarations is sneakily a separate thing, it means to run the provided tactic if the field is omitted

[Mario Carneiro]

if you put the by block in parentheses it will be checked at declaration time instead, as a regular optional value

[Michal Wallace (tangentstorm)]

aha! thanks!

[Mario Carneiro]

But I think you don't really want that either, especially for proof fields. If you can prove it outright from earlier fields then there is no point putting it in the structure

[Mario Carneiro]

This particular example of course can't be proved outright, since an arbitrary function C : α → NPrime need not have its values all greater than 234234

[Mario Carneiro]

if you think that C will normally be given in some explicitly computable way such that by decide or by omega or something will prove this side goal then you can use := by decide for that

[Michal Wallace (tangentstorm)]

Well, again, I'm thinking about this from a software perspective. I think of a class as something like an interface in java or C#, where there might be many sieves lying around. I want to define some hypotheses / contracts that any algorithm has to fufill, and then provide one proof that demonstrates that those conditions are enough to give you a working prime sieve.

[Michal Wallace (tangentstorm)]

(I was just making up an assertion that was obviously false for that example)

[Mario Carneiro]

Can this proof be conducted only using the assumptions in the structure?

[Mario Carneiro]

if so, the proof should not be in the structure itself

[Mario Carneiro]

the structure is where you put the initial data and assumptions about that data that every implementer must prove

[Michal Wallace (tangentstorm)]

here's what i have currently:

import Mathlib.Data.Nat.Prime

def NPrime : Type := { n: Nat // Nat.Prime n } deriving Repr
instance : ToString NPrime where
  toString s := s!"{s.val}"

-- naturals less than C
def nltC (c:NPrime) (n:Nat) : Prop := n < c.val
-- primes less than C
def pltC (c:NPrime) (p:NPrime) : Prop := p.val < c.val
-- natural numbers coprime to some known primes:
def cpks (ks:Set NPrime) (n:Nat) : Prop :=
  n ≥ 2  ∧ ∀ p ∈ ks, ¬(p.val ∣ n)

-- a generic prime sieve
class PrimeSieve (α : Type u) where
  -- implement your algorithm in terms of these:
  C : α → NPrime
  init : α
  next : α → α

  -- this class defines some sets for you to
  -- refer to when proving that it works:

  -- S: the set of "known primes", less than C
  S  (x:α) : Set NPrime := pltC (C x)
  hS (x:α) : ∀ p ∈ S x, pltC (C x) p := (by sorry)

  -- R is the remainding nats, coprime to all p∈S
  R  (x:α) : Set Nat  := cpks (S x)
  hR (x:α) : ∀ n ∈ R x, cpks (S x) n := (by sorry)

  -- these are tne steps you need to prove:
  -- apostrophe indicates result of the 'next' operation
  hCinR (x:α) : (C x).val ∈ R x             -- C is in R
  hCinS (x:α) : (C x) ∈ S (next x)          -- C is in S'
  hRmin (x:α) : ∀ n ∈ (R x), n ≥ (C x).val  -- C is min of R
  hSmax (x:α) : ∀ p ∈ (S x), p.val < (C x).val  -- C is max of S'
  hNewC (x:α) : (C $ next x).val > (C x).val -- C' > C

open PrimeSieve
-- demonstrate that (hS, hR, hMin, hNew) are enough to prove
-- that a sieve generates the next consecutive prime at each step
theorem hs_suffice (α : Type u) [PrimeSieve α]
  (x:α) (x':α) (hnx: x' = next x)
  : (C x').val > (C x).val  -- "we have a new, bigger prime"
  ∧ ¬∃ p:NPrime,   -- "and there is no prime between them"
    ((C x).val<p.val ∧ p.val<(C x').val)
  := by
    apply And.intro
    case left => -- "new bigger prime" is class invariant
      rw[hnx]; apply hNewC
    case right =>
      by_contra ex_p; let ⟨p, hp⟩ := ex_p -- assume p exists.
      have : cpks (S x') (C x').val := hR x' (C x').val $ hCinR x'
      simp[cpks] at this
      have ⟨hc'g2, hc'coprime⟩ := this
      -- if there's a prime between C and C', then
      -- C'
      have h: (C $ next x).val > (C x).val := hNewC x
      rw[←hnx] at h
      have : (C x).val < p.val := sorry
      have : p.val < (C x').val := sorry

[Michal Wallace (tangentstorm)]

I'm thinking that for any sieve implementation, you want to somehow talk about these sets... Otherwise you're tightly coupling the proof to the implementation.

[Michal Wallace (tangentstorm)]

but the sets depend on a value in the structure (C, the "current prime")

[Michal Wallace (tangentstorm)]

This does feel like a clunky way to structure it, though. Maybe a class isn't the right tool for the job?

[Michal Wallace (tangentstorm)]

(afk for a couple hours, will be back on this tonight)

[Mario Carneiro]

Here's how I would structure this class, since you have some fields, some definitions and then some more fields:

class PrimeSieveBase (α : Type u) where
  -- implement your algorithm in terms of these:
  C : α → NPrime
  init : α
  next : α → α
open PrimeSieveBase

/-- S: the set of "known primes", less than C -/
def S [PrimeSieveBase α] (x:α) : Set NPrime := { p | p.val < (C x).val }

/-- R is the remainding nats, coprime to all p∈S -/
def R [PrimeSieveBase α] (x:α) : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ S x, ¬(p.val ∣ n) }

-- a generic prime sieve
class PrimeSieve (α : Type u) extends PrimeSieveBase α where
  -- these are the steps you need to prove:
  -- apostrophe indicates result of the 'next' operation
  hCinR (x:α) : (C x).val ∈ R x             -- C is in R
  hCinS (x:α) : (C x) ∈ S (next x)          -- C is in S'
  hRmin (x:α) : ∀ n ∈ (R x), n ≥ (C x).val  -- C is min of R
  hSmax (x:α) : ∀ p ∈ (S x), p.val < (C x).val  -- C is max of S'
  hNewC (x:α) : (C $ next x).val > (C x).val -- C' > C

[Mario Carneiro]

Here's an alternative structuring which eliminates the type alpha entirely and just uses an abstract state:

/-- S: the set of "known primes", less than C -/
def S (C : Nat) : Set NPrime := { p | p.val < C }

/-- R is the remainding nats, coprime to all p∈S -/
def R (C : Nat) : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ S C, ¬(p.val ∣ n) }

structure PrimeSieveOK (C : Nat) : Prop where
  -- implement your algorithm in terms of these:
  hCp : Nat.Prime C             -- C is prime
  hCinR : C ∈ R C               -- C is in R
  hRmin : ∀ n ∈ R C, n ≥ C      -- C is min of R
  hSmax : ∀ p ∈ S C, p.val < C  -- C is max of S'

structure PrimeSieveNext (C : NPrime) (C' : Nat) : Prop where
  hCinS : C ∈ S C'       -- C is in S'
  hNewC : C' > C.val     -- C' > C

theorem PrimeSieveNext.no_prime_between
    (hC : PrimeSieveOK C) (next : PrimeSieveNext ⟨C, hC.hCp⟩ C') :
    ¬∃ p:NPrime, C<p.val ∧ p.val<C' := sorry

theorem PrimeSieveNext.next
    (hC : PrimeSieveOK C) (next : PrimeSieveNext ⟨C, hC.hCp⟩ C') :
    PrimeSieveOK C' := sorry

(I was going to have a structure for the state but in this formulation there is only one number in the state, the current prime value C, so it's overkill.) The idea is that in practice you would prove that your prime sieve relates to some abstract state C satisfying PrimeSieveOK, and next returns a value which is related via PrimeSieveNext

[Mario Carneiro]

At this point though PrimeSieveOK is just some predicate on natural numbers which probably has a name already. Indeed:

• hCinR is provable outright, because C, being prime, is coprime to all primes less than it
• hRmin is false because it contains (non-prime) natural numbers less than C, but if you add n >= C instead of n >= 2 to the definition of R then it becomes trivially provable outright
• hSmax is trivially true by definition

So all of the additional assumptions are unnecessary and PrimeSieveOK C is just another name for Nat.Prime C (and similarly PrimeSieveNext has two assumptions which are redundant with each other), suggesting a reformulation as:

theorem hCinR (h : Nat.Prime C) : C ∈ R C := sorry               -- C is in R
theorem hRmin (h : Nat.Prime C) : ∀ n ∈ R C, n ≥ C := sorry      -- C is min of R
theorem hSmax (h : Nat.Prime C) : ∀ p ∈ S C, p.val < C := sorry  -- C is max of S'

theorem PrimeSieveNext.next
    (hC : Nat.Prime C) (next : C < C') :
    Nat.Prime C' := sorry

although at this point we're clearly missing some assumption with which to prove C' is prime; most likely you want to say that C' is in R C or something along those lines

[Michal Wallace (tangentstorm)]

Oh, wow. Thanks for all this! Lots to process. :)

[Michal Wallace (tangentstorm)]

hRmin is false because it contains (non-prime) natural numbers less than C,

Are you sure? Because if you're right, I'm misunderstanding something about my definition:

def R (C : Nat) : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ S C, ¬(p.val ∣ n) }

I think this is a "set comprehension" that says "the set of natural numbers that greater than or equal to two, and are co-prime to every (prime) number in S"... So if C is 11, a non-prime number smaller than C (like 4) can't be in the set, because 4 is divisible by 2, and 2 ∈ S C.

[Michal Wallace (tangentstorm)]

structure PrimeSieveOK (C : Nat) : Prop where ... is warping my mind a little bit. It's a structure that can be applied like a function, and when applied, returns a Prop ?

[Mario Carneiro]

PrimeSieveOK C is a predicate, like Nat.Prime C. When your predicate is a big conjunction of things structure is a natural way to define it and give the sub-properties names

[Mario Carneiro]

Michal Wallace (tangentstorm) said:

I think this is a "set comprehension" that says "the set of natural numbers that greater than or equal to two, and are co-prime to every (prime) number in S"... So if C is 11, a non-prime number smaller than C (like 4) can't be in the set, because 4 is divisible by 2, and 2 ∈ S C.

You are right, that was my mistake. R will not have any numbers less than C in it

[Michal Wallace (tangentstorm)]

Mario Carneiro said:

PrimeSieveOK C is a predicate, like Nat.Prime C. When your predicate is a big conjunction of things structure is a natural way to define it and give the sub-properties names

Very cool! I guess that's similar to __call__ in python, but for some reason I didn't expect to find something like that in a functional language.

[Mario Carneiro]

I'm not sure what you mean: Lean does have an analogue of __call__ but this isn't it

[Mario Carneiro]

This is just "propositions are types"

[Mario Carneiro]

structure can be used to make data, and it can also be used to make propositions

[Mario Carneiro]

and structure can depend on parameters, which makes parametric data types when applied to data and predicates when used with propositions

[Mario Carneiro]

another aspect of this is the fact that you can have properties next to data fields in a structure, and lean doesn't differentiate between them, they are all just fields in the structure but some are data values and some are proofs

[Michal Wallace (tangentstorm)]

hSmax is trivially true by definition
I couldn't see that before, because that was buried way down in the definition of pltC... Your rewrite makes the available facts much clearer.

[Mario Carneiro]

By the way, the analogue of __call__ is CoeFun, for example if we wanted your ASeq to act like a function:

structure ASeq where  -- arithmetic sequence (k + dn)
  k : Nat  -- constant
  d : Nat  -- difference
deriving Inhabited, Repr

-- apply formula to n
def ap (s : ASeq) (n : Nat) : Nat := s.k + s.d * n

instance : CoeFun ASeq fun _ => Nat → Nat := ⟨ap⟩

def mySeq : ASeq := { k := 5, d := 3 }
#eval (mySeq 0, mySeq 1, mySeq 2) -- (5, 8, 11)

[Michal Wallace (tangentstorm)]

Oh, that's nice. That's definitely going in!

I see the difference between this and the previous code (namely this has an implementation), but in my mind, I read struct Foo (arg:Nat) : Prop as "Foo takes a Nat and returns a Prop"

It seems like the only reason you can't say structure ASeq (n:Nat) : Nat where... to do what you just did is because structure doesn't have anywhere to put the := ⟨ap⟩.

But... You can do that with Prop because prop doesn't require an implementation?

Or is : Prop somehow special and means "take the AND-sum of all the proposition fields in this structure" ?

[Mario Carneiro]

The reason you can't write structure ASeq (n:Nat) : Nat where is because structure declares a type, and the type of a type must be Prop or Type u for some u (these two cases are subsumed by the more general Sort u, where Prop = Sort 0 and Type u = Sort (u+1))

[Mario Carneiro]

That is, you should read Prop as a kind of type universe, similar to Type and unlike Nat

[Mario Carneiro]

You can also write structure Foo : Type 37 where if you want, although there aren't many reasons to use such big universes

[Michal Wallace (tangentstorm)]

Got it. So it is a bit like an AND of all the facts, but it isn't magic... If you have a value of this type, you must have instances of all those facts attached, or it wouldn't type check.

[Mario Carneiro]

(and if you want the OR of facts, you will want to look into inductive, which structure is actually syntax sugar for, and which also works to make both data and propositions)

[Michal Wallace (tangentstorm)]

Okay. I like what you've done here. I think no_prime_between still has to be part of the spec, because otherwise I think next $ StateWhereCis 5 might return next $ StateWhereCis 11 ... Like it might internally track 7, but you would never be able to observe C=7 because it got skipped over.

[Mario Carneiro]

right, the idea was that this goes into PrimeSieveNext, although it's not clear how you want to state this property in a way that doesn't make it tautological

[Mario Carneiro]

but if it's just one or two assumptions you might just not bother with a structure and just have a theorem that says everything you want from the abstraction

[Michal Wallace (tangentstorm)]

I also think maybe this spec ought to be called something like PrimeGenerator because it isn't necessarily a sieve, and whether it is or not, who cares? (Like, given C, you could just test whether each subsequent number is prime or not, without ever remembering any primes)

[Mario Carneiro]

I did not try very hard at the names :)

[Michal Wallace (tangentstorm)]

Oh, I didn't mean you. I'm the one calling it a Sieve the whole time. :)

[Mario Carneiro]

I agree at this level of abstraction it's not very sieve-like

[Mario Carneiro]

the concrete implementation was sieve like but this is just a sequence of primes

[Mario Carneiro]

in fact with my refactoring it's not even really a sequence, it's just a relation between adjacent primes

[Michal Wallace (tangentstorm)]

Do I even need sets? I was sort of imagining these prop fields as induction hypotheses... So at each step (S' := S ∪ {C}) and I just have to say something about removing the multiples of C from R. But maybe I can say this with types?

[Mario Carneiro]

sets are just sugar over predicates, but it can be nice to have the notation of e.g. adding an element like you just wrote

[Michal Wallace (tangentstorm)]

I don't know if anyone would be interested in hanging out for this, but I'm live-streaming some work on this today. https://youtube.com/live/MNz6R4g0pOU?feature=share

[Michal Wallace (tangentstorm)]

Well, that probably wasn't worth watching since it was just me fumbling around and looking dumb, but here's the main result: a second prime generator that just uses Nat.find.

section simple_gen

  def prime_gt (c : Nat) (p : Nat) : Prop :=
    Nat.Prime p ∧ c < p

  instance : Decidable (prime_gt c p) := by
    rw[prime_gt]; apply inferInstanceAs

  theorem ex_prime_gt (c:Nat) : ∃ p, prime_gt c p := by
    simp[prime_gt]
    let d := c + 1 -- because the line below has ≤ and we need <
    let ⟨p, hcp, hprime⟩ : ∃ (p : ℕ), d ≤ p ∧ Nat.Prime p :=
      Nat.exists_infinite_primes d
    use p; apply And.intro
    · exact hprime
    · linarith

  -- Use Nat.find to get the smallest n that satisfies P
  def next_prime (c:Nat) : Nat := Nat.find <| ex_prime_gt c

  def PrimeGt (c:Nat): Type := { p : Nat // prime_gt c p }

  -- !! seems like this would do the search twice?
  def next_primegt (c: Nat) : PrimeGt c :=
    ⟨Nat.find (ex_prime_gt c), Nat.find_spec (ex_prime_gt c)⟩

  def nprime_gt (p:PrimeGt c): NPrime :=
    ⟨p.val, (by
      have : prime_gt c p.val := p.property
      simp[prime_gt] at this
      exact this.left)⟩

  def next_nprime (c:Nat) : NPrime :=
    nprime_gt (next_primegt c)

  structure SimpleGen where
    c : NPrime

  instance : PrimeSieveState SimpleGen where
    C x := x.c
    init := ⟨2, Nat.prime_two⟩
    next x := ⟨nprime_gt (next_primegt x.c.val)⟩

end simple_gen

[Michal Wallace (tangentstorm)]

I broke the code up into separate files, too.. the prime generator spec is now here: https://github.com/tangentstorm/treesiv/blob/main/PrimeGen.lean

[Michal Wallace (tangentstorm)]

So. This one is not a sieve, and having it to compare kind of helps me think a little. The core thing I really want to show is that from a given state of the system, I can generate the next prime in sequence.

For this simple generator, the entire state of the algorithm is contained in C (the current prime). It doesn't remember any other primes that it has found, and it doesn't do any sieving, so there's no need to refer to those sets R and S.

So... I guess my goal now is to make a PrimeGenOk type that just demands a proof that C' is the next consecutive prime after C, and then prove that this simple generator succeeds here. (Looks like Nat.find_min will make this easy enough.)

From there, I think I make a general rule about prime sieves using the sets I've been talking about. The argument would be: if you track these three things C (a prime), S (primes less than C) and R (all naturals greater than 1 and coprime to C), and at each step, C' > C, then your prime generator is a sieve.

[Michal Wallace (tangentstorm)]

Does it matter if it's a sieve? Probably not, but I still have a hunch this framework is going to make it easier to prove that an arbitrary sieving algorithm generates the next prime from its internal state. In fact, maybe you only have to implement a proof that you're returning the smallest natural number in R, and the "prime sieve machinery" provides the proof that you have the next prime.

[Michal Wallace (tangentstorm)]

---

I feel like I'm down to one last hurdle, but I'm still a tiny bit stuck. Here's where I'm at:

• APrimeGenerator is an algorithm that generates a natural number P at each step, and proves:
  -- P is prime
  -- the next P (P') is greater than the current P
  -- There is no other prime number between P and P'

• A PrimeSieve is a PrimeGenerator that generates P' by producing some C : Nat and proving that C is the smallest Nat not divisible by any prime less than or equal to the original P.

Obviously, every PrimeGenerator will generate a number of this kind since that's one of the definitions of a prime. But some algorithms specifically use data structures that model the set of remaining coprimes, and specifically reduce the set at each step by filtering out multiples of known primes.

So now the problem is moving from math world to programming world.

I have two actual sieve algorithms in lean4 - mine and ykonstant1/esiv.

Esiv inspects a bitmask at each step so it potentially generates multiple primes at each step. It does not provide a proof of correctness, and I think with a little work I could refactor it into a form that fits my rule above.

Neither algorithm explicitly mentions "the set of naturals coprime to the primes found so far" but they both manipulate some data structure from which the minimal value of this set can be extracted.

I think the correctness proof proceeds by induction, but not using the induction keyword, because it's induction on function application rather than a recursive data type: in order to implement PrimeSieve, you have to provide init and next functions that produce a value of the type... And therefore have to provide all the necessary proofs.

Since next has type α → α, we can use any proof attached to the input value as the induction hypothesis.

I'm still unclear how to actually do it, though...

[Michal Wallace (tangentstorm)]

My instinct is to break it into several pieces, and prove:

1. there is some way to produce an element of R (the set of un-scratched-off numbers) from the data
2. there is a way to produce the minimal value of this set.
3. there is some change being made to the data that filters out the multiples of the newly found prime

[Mario Carneiro]

By the way, if you are collecting sieve implementations, I posted a prime sieve method (with a proof of correctness) here not too long ago: https://leanprover.zulipchat.com/#narrow/stream/270676-lean4/topic/Checking.20the.20Goldbach.20conjecture/near/429974821

[Michal Wallace (tangentstorm)]

have : buf[i] = true ↔ i.Prime := by

Aha!! This is exactly the sort of thing I need.

[Michal Wallace (tangentstorm)]

I'm aiming for the PrimeSieve class to actually provide the .Prime part based on some lower-level fact, but I think I see the invariant I want to track now.

• show that some operation (extract) on the state produces only numbers in R
• let R' be R without the multiples of P
• show that there is some operation called filter that changes the state so that extract now only produces numbers from R' (possibly using the previous proof as an induction hypothesis)
• show that there is some operation extractMin to extract the minimum value in R'

(off to study your proof in detail now...)

[Michal Wallace (tangentstorm)]

Oh, interesting. @Kevin Buzzard, you mentioned Oliveira e Silva's Goldbach verification project in the thread Mario just linked. He also has a similar project for the Collatz conjecture.

Me using arithmetic sequences to model a prime sieve is rather silly and inefficient, but the actual reason I was thinking about modelling arithmetic sequences was to illustrate how his Collatz verifier works.

It's essentially building a binary tree of arithmetic sequences in a top down, breadth-first manner, and then marking some of the nodes as dead before generating the next level.

[Jireh Loreaux]

My go-to reference for an actually efficient sieve of Erastosthenes is this. She also explains what makes most of the "standard implementations" rubbish.

[Michal Wallace (tangentstorm)]

Jireh Loreaux yeah, the distinction she's making is kind of what I'm trying to get at with the separate definitions of PrimeGenerator and PrimeSieve...

I find it interesting that there doesn't seem to be any way to define a proposition in Lean that would distinguish the different cases.

[Michal Wallace (tangentstorm)]

Her idea of using a priority queue of iterators is kinda neat.

[Michal Wallace (tangentstorm)]

I think I have a plan now:

---

Informal proof that we have a mapping from natural numbers to the residual set.

First, we show that the array of sequences ss represent a map from Nat to the residual set.

If q is the number of sequences, then we can map an input n to ss[q%n] (q/n)

In an arithmetic sequence A d k, if prime p | d, then p | A n iff p | k

Since we use a primorial for our delta, the individual sequences either never or always divide the primorial’s factors.

So the proof is by induction: first we know that in the input state, we had already filtered out all the other known primes. So all we have to do is filter out the “new” prime.

Informal proof that we can filter the set.

The filtering procedure is:

• First, “partition” the sequence by composing it with the sequences (i + p) for i in range(p). 

• This produces p new sequences, with a common delta divisible by p.

• Of these p sequences, exactly 1 will produce multiples of p. Remove this one.

Informal proof that we can find the minimal value:

We track a list/set of arithmetic sequences (ASeq) that all have the same delta.

Delta is a natural number, so term n <= term (n+1), and term 0 is the minimal term.

Since all sequences in our list have the same delta, the minimal value we can produce is term 0 of the sequence with smallest k.

[Michal Wallace (tangentstorm)]

---

It seems like it might be useful to have a notion of a Comb structure that represents a set or list of arithmetic sequences that all have the same delta. ("Comb" because if you plot the lines formed by the individual sequences, they’re parallel so it looks like running a comb through the sand).

Partitioning a sequence would produce a comb, and you should be able to merge combs and sequences that have the same delta. (Really, a sequence could just be a special case of a comb that has only one tooth) 

( "Rake" maybe? Apparently this might also be called a Delta-set in math?)

[Michal Wallace (tangentstorm)]

Progress: I created a Rake type, constructed an "identity rake" and showed that it maps the natural numbers to themselves by proving theorem idr_id (n:Nat) : (idr.term n = n).

[Michal Wallace (tangentstorm)]

The rake works like a function from Nat -> Nat. I suspect that in order to demonstrate that a particular Rake models a given subset of the naturals, I'm going to have to show that this function is a bijection.

[Michal Wallace (tangentstorm)]

(It occurs to me that a Rake might be exactly the same concept as what sieve theory calls a Wheel.)
edit: It is not, but can be used to make one.

[Michal Wallace (tangentstorm)]

I implemented aseq.gte n, which fiddles with the sequence to give you only the results greater than or equal to n:

#eval terms evens 10         -- [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
#eval terms (evens.gte 5) 10 -- [6, 8, 10, 12, 14, 16, 18, 20, 22, 24]

I proved that the results are >= n, but didn't yet prove that I haven't accidentally removed a result I shouldn't have.

[Michal Wallace (tangentstorm)]

I've also now got the concept of a RakeMap which makes the connection between a given Rake and a Prop that defines the corresponding set of natural numbers explicit:

So, for example, here's a proof that the idr (identity rake) models the set { n : Nat | True }

def idr : Rake := { -- the identity rake (maps n -> n)
  d := 1,
  seqs := [dseq 0 1],
  hsort := List.sorted_singleton _,
  hsize := Nat.zero_lt_one }

structure RakeMap (prop: Nat → Prop) where
  rake : Rake
  hbij : ∀ n, prop n ↔ ∃ m, rake.term m = n

-- proof that idr represents the identity function
def idrm : RakeMap (λ _ => True) := {
  rake := idr,
  hbij := by intro n; simp[Rake.term, idr]; unfold dseq; unfold term; simp }

[Michal Wallace (tangentstorm)]

The plan is to now implement RakeMap.gte n that filters out results less than n (using the ASeq.gte function that I just mentioned), and something like RakeMap.remove_multiples (p:NPrime), which will do what it says by using partition to split each sequence into p sequences with a common delta divisble by p, and then filter out the ones whose constant term is a multiple of p.

I'm pretty sure I will be able to show that both gte and remove_multiples preserve whatever property was there before, but also filter out the appropriate values. Since the identity idrm : RakeMap models the entire set of naturals, I will then be able to reduce that to those > 2, and filter out multiples of any prime, and (because the lists are sorted and I have a proof that the Rakemap.term 0 always returns the minimal element of the modeled set... Then I will have all the pieces necessary to demonstrate a formally verified prime sieve.

(yay!)

[Michal Wallace (tangentstorm)]

It looks like I'm going to have some minor hurdles to overcome with the "the tines of the rake are sorted" proof, but I will worry about that tomorrow.

[Michal Wallace (tangentstorm)]

Small update:

I got the issues with proving the list is sorted worked out in another thread.

I'm still working on RakeMap. I think I have a decent plan here. I already have the "identity rake" above. Now I need to show that for gte and remove_multiples, three facts hold:

section gte_lemmas

  variable (rm: RakeMap prop) (p n: Nat)

  lemma RakeMap.gte_drop -- gte drops terms < p
    : n < p → ¬(∃m, (rm.rake.gte p).term m = n) := by
    sorry

  lemma RakeMap.gte_keep -- gte keeps terms ≥ p
    : n≥p ∧ (∃pm, rm.rake.term pm = n) → (∃m, (rm.rake.gte p).term m = n) := by
    sorry

  lemma RakeMap.gte_same -- gte introduces no new terms
    : (∃pm, (rm.rake.gte p).term pm = n) → (∃pm, rm.rake.term pm = n) := by
    sorry

end gte_lemmas

From those three laws, I can prove the main result for the operation:

def RakeMap.gte_prop (prev : RakeMap prop) (p: Nat)
  : RakeMap (λ n => prop n ∧ n ≥ p) :=
  let rake := prev.rake.gte p
  let proof := by
    intro n; symm
    let hm : Prop := (∃m, rake.term m = n)
    let hpm : Prop := (∃pm, prev.rake.term pm = n )
    have : prop n ↔ hpm := prev.hbij n
    have : n<p → ¬hm := gte_drop prev p n
    have : n≥p ∧ hpm → hm := gte_keep prev p n
    have : hm → hpm := gte_same prev p n
    by_cases n<p; all_goals aesop
  { rake := rake, hbij := proof }

The logic for proving remove_multiples should be almost exactly the same. If there were more than two operations, I might consider making a RakeMapRefinement class, but since there are only two, I'll probably just write out both proofs.

[Michal Wallace (tangentstorm)]

Doing another live-proving stream in a bit if anyone cares to come hang out. https://youtube.com/live/CVJ7aZ6tRqY?feature=share

[Michal Wallace (tangentstorm)]

Mario Carneiro said:

• hSmax is trivially true by definition

This was quite tricky. You were correct. My initial formulation was:

  hCinR (x:α) : (C x).val ∈ R x             -- C is in R
  hCinS (x:α) : (C x) ∈ S (next x)          -- C is in S'
  hRmin (x:α) : ∀ n ∈ (R x), n ≥ (C x).val  -- C is min of R
  hSmax (x:α) : ∀ p ∈ (S x), p.val < (C x).val  -- C is max of S'
  hNewC (x:α) : (C $ next x).val > (C x).val -- C' > C

The problem was that hSmax didn't express what my comment says it did. It should have said S <| next x.

My goal was to prove that from laws like these:

theorem hs_suffice (α : Type) [PrimeGen α] [PrimeSieve α] (x:α)
  :  (C' x > C x)  ∧ (¬∃ p:NPrime, C x < p ∧  p < C' x) := by _

I was unable to prove this until I finally re-added the informal rule in the comment. The final rules look more or less like this:

class PrimeGen (α : Type u) where
  C : α → NPrime
  init : α
  next : α → α

open PrimeGen
def C' {α : Type u} [PrimeGen α] (x: α) : NPrime := C <| next x

section prime_sieve
  variable {α : Type} [PrimeGen α]  (x: α)

  -- S: the set of "known primes", less than C
  def S : Set NPrime := { p | p < (C x) }

  -- R: the set of "remaining" numbers, coprime to all known primes
  -- describe the set of naturals with no prime factors less than some c
  def R : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ (S x), ¬(p.val ∣ n) }

  def S' : Set NPrime := S <| next x
  def R' : Set Nat := R <| next x
end prime_sieve

class PrimeSieve (α : Type) [PrimeGen α] : Prop where
  hRmin  (x:α) : ∀n ∈ R x, C x ≤ n    -- C is min of R
  hCmax  (x:α) : ∀ p ∈ S' x, C x ≥ p  -- C is the max of S'
  hNewC  (x:α) : C x < C' x           -- C < C'

[Michal Wallace (tangentstorm)]

I may break the definition of PrimeGen up so that it generates C x : Nat and separately provides a proof of Nat.Prime C x ... I believe have enough scraps lying around to prove that C x is still prime for allPrimeSieve instances, given the above definitions.

[Michal Wallace (tangentstorm)]

In trying to implement instance : PrimeSieve RakeSieve I realized I had to re-add most of the "trivial" statements I threw away before... Only instead of defining the current prime C in terms of the sets I'm tracking, the correct approach was to define the sets in relation to some highest known Nat.Prime P, and then have the PrimeSieve implementation produce a C:Nat together with a proof that C has the generic properties a prime sieve would look for.

The PrimeSieve class itself provides the generic proof that C is actually the next prime.

So... Now PrimeGen just generates a prime P at each step: (the change was just renaming C to P):

class PrimeGen (α : Type) where
  P : α → NPrime
  init : α
  next : α → α

And then:

section prime_sieve
  variable {α : Type} [PrimeGen α]  (x: α)
  -- S: the set of "known primes", ≤ P
  def S : Set NPrime := { p | p ≤ (P x) }
  -- R: the set of "remaining" numbers, coprime to all known primes
  -- describe the set of naturals with no prime factors less than some c
  def R : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ (S x), ¬(p.val ∣ n) }
end prime_sieve

/--
A PrimeSieve is a PrimeGen that uses a sieve to generate primes.
In practice, this means that it somehow models the set of natural
numbers greater than 1 that are coprime to a list of known primes,
and then repeatedly:
  - obtains the minimum of this set as the next prime
  - eliminates multiples of the new prime. -/
class PrimeSieve (α : Type) [PrimeGen α] where
  -- P (g:α) : NPrime
  C (g:α) : Nat -- the "current" or "candidate" prime
  hSmax  (g:α) : ∀ p ∈ S g, P g ≥ p   -- P is the max of S
  hCinR  (g:α) : C g ∈ R g            -- C is an element of R
  hRmin  (g:α) : ∀ n ∈ R g, C g ≤ n   -- C is min of R
  hCgtP  (g:α) : C g > P g            -- C > P
  hP'C   (g:α) : P' g = C g           -- P' = C

[Michal Wallace (tangentstorm)]

and from those definitions, i was able to prove:

-- demonstrate that g generates the next consecutive prime at each step.
-- "we have a new, bigger prime, and there is no prime between them".
theorem hs_suffice (α : Type) [PrimeGen α] [PrimeSieve α] (g:α)
  :  (Nat.Prime (C g)) ∧ (C g > P g) ∧ (¬∃ q:NPrime, P g < q ∧  q < C g) := by
    split_ands
    · exact c_prime α g
    . exact PrimeSieve.hCgtP g
    · exact no_skipped_prime α g

[Michal Wallace (tangentstorm)]

I really like this because it formally verifies the basic concept of a Prime Sieve, completely independent of the specific data structure or algorithm that's doing the work.

[Michal Wallace (tangentstorm)]

Current snapshot: https://github.com/tangentstorm/treesiv/blob/856e41a57f7ab5897d63b532934b4c666b7060af/PrimeSieve.lean

[Michal Wallace (tangentstorm)]

Actually hSmax looks like it might be trivial again, but the proofs are using it. I'll see if I can get rid of it tomorrow. Then I just need to expose the right parts of the RakeSieve algorithm to prove that it meets these criteria.

[Michal Wallace (tangentstorm)]

an unexpected problem

Uh-oh. I may need to rethink something here. In the previous design, the sieve provided its own proof that C was prime at each step.

Now, implementing class PrimeSieve and proving that your C:Nat has certain properties is meant to be sufficient to show that C is prime.

But this all depends on the sieve having access to a maximum known prime P. The idea was that each step, C is proven prime and becomes the new P.

Currently, in order to be a PrimeSieve, a type has to first be a PrimeGen that generates a prime P.

In order to implement PrimeGen, I just need to produce any prime. Okay.

But in order to implement PrimeSieve, that prime has to be the C from the previous step.

But I can't prove that C is prime without implementing PrimeSieve.

This is backwards.

The sieve should not have to be a PrimeGen first. Instead:

• class PrimeSieve a should be structure PrimeSieve a [SieveState a]
• C and all the certifications about C move to this SieveState type.
• SieveState's next function should take as an argument a proof that C is prime

Then maybe make PrimeGen a bit stronger so that it actually generates the primes in sequence, and have the final proof be that PrimeSieve implements PrimeGen.

[Michal Wallace (tangentstorm)]

:check: The fix is done:

class SieveState (α : Type) where
  P     (s:α) : NPrime
  C     (s:α) : Nat
  next  (s:α) (hC: Nat.Prime (C s)) : α
  hNext (s:α) (hC: Nat.Prime (C s)) : (P (next s hC) = C s)

• class PrimeSieve α, S, R require only SieveState α
• there is no longer a need for hSmax
• PrimeGen requires a proof that there are no gaps between generated primes
• SimpleGen (a PrimeGen based on Nat.find) supplies such a proof (this was quite educational)
• there is a basic outline for RakeSieve which implements SieveState.

Current snapshot is here: https://github.com/tangentstorm/treesiv/blob/6f1ab7145c1fc1e5e6f773cbdf8b6a033b1d907b/PrimeSieve.lean

[Michal Wallace (tangentstorm)]

What's left to do

• prove that I can model a bijection to the set R by removing multiples of each new prime from the RakeMap
• that will require showing that the existing functions can isolate the arithmetic sequences that produce multiples of the prime (they can) and that I correctly detect and remove them.
• implement instance PrimeGen PrimeSieve(all the important lemmas are proven, but I still have to write the instance).
• prove 3 lemmas about RakeMap.gte n (that it can filter out terms less than n)
• prove 3 lemmas about RakeMap.rem p (that it removes multiples of a prime)
• implement instance PrimeSieve RakeSieve

It's possible PrimeSieve might wind up having to be a structure rather than a class, but the basic proof requirements are the same either way.

After that, I have some optimizations in mind, but the more important thing is to write and animate the video I'm making about all this for 3blue1browns' Summer of Math Exposition contest (which, aside from just wanting to learn lean, was the whole point of this project, and why I keep posting a devlog here).

Anyway, I may do a live-proving stream tonight after my kids go to bed.

[Michal Wallace (tangentstorm)]

Doing live stream for a bit. Currently rebooting mid stream because this 10 iterations of this sieve crashed my whole machine :)

[Michal Wallace (tangentstorm)]

back online and streaming for an hour or two at https://www.youtube.com/live/wU3RisZt9tk

[Michal Wallace (tangentstorm)]

The major result of last night's stream was to remove the need for a SieveState to provide the proof that candidate prime C is greater than highest known prime P:

hCgtP (g:α) : C g > P g -- C > P

I still need this fact, but I was able to prove it from the other definitions.

[Michal Wallace (tangentstorm)]

This brings the "burden of proof" for a SieveState down to this:

class SieveState (α : Type) where
  P     (s:α) : NPrime
  C     (s:α) : Nat
  next  (s:α) (hC: Nat.Prime (C s)) : α
  hNext (s:α) (hC: Nat.Prime (C s)) : (P (next s hC) = C s)

class PrimeSieve (α : Type) [SieveState α] where
  hCinR  (g:α) : C g ∈ R g            -- C is an element of R
  hRmin  (g:α) : ∀ n ∈ R g, C g ≤ n   -- C is min of R

hNext is just a contract that you have to always make your (state.next hC) set the new P to the old C, so that's trivial.

Now this is looking a bit messy to me:

  variable {α : Type} [SieveState α]  (x: α)
  -- S: the set of "known primes", ≤ P
  def S : Set NPrime := { p | p ≤ (P x) }
  -- R: the set of "remaining" numbers, coprime to all known primes
  -- describe the set of naturals with no prime factors less than some c
  def R : Set Nat := { n | n ≥ 2 ∧ ∀ p ∈ (S x), ¬(p.val ∣ n) }

I think I will simplify R to a single function R from Nat to Set Nat, (like mario suggested), then just inline S (because it's no longer used anywhere), and then maybe provide a SieveState.Rp that just returns (R x.P).

I'm a little worried that refactoring will break all my proofs, but I guess I can just try it and if it's too hard to change, I'll roll it back.

(tracking as issue 1)

[Michal Wallace (tangentstorm)]

Actually, this seems like it might actually be a blocker.

It's the same sort of circular dependency mistake I made with PrimeGen and PrimeSieve... This time, in order to prove my RakeSieve is a SieveState, I have to make statements about the set R, but R is un-necessarily defined in terms of SieveState, so I can't use the definition of R.

[Michal Wallace (tangentstorm)]

What is this mistake I keep making?

I suppose it's a form of circular reasoning. I define an abstraction in terms of some property I want it to have, but then either:

• assume i have the property when i should instead be proving it, or
• tightly couple the property's definition to the abstraction's definition, so I can't prove anything has the property without it already being the thing I want to prove that it is.

Well that's a mouthful, and maybe there are two different mistakes here, but at least I can watch out for this in the future.

[Michal Wallace (tangentstorm)]

---

Decoupling R and getting rid of S was easy.

-- R: the "residue", or "remaining" nats coprime to all primes < p.
-- In a sieve, these are the numbers that haven't yet been "sifted out."
def R (P:NPrime): Set Nat := { n | n ≥ 2 ∧ ∀ p ≤ P, ¬(p.val ∣ n) }

Most of the proofs were fixed by just replacing R g, withR(P g) everywhere.

[Michal Wallace (tangentstorm)]

---

Every time I try to demonstrate that RakeSieve is a SieveState, I find a new problem or learn something about Lean.

For example, this is not possible to prove:

theorem RakeSieve.hCinR  (g:RakeSieve) : C g ∈ R (P g)            -- C is an element of R
  := sorry

Why not? Well, because it doesn't follow from anything specific about the RakeSieve structure. Rather it's a statement that would have to cover every possible way of constructing a RakeSieve, and this is an open set.

I could prove this theorem for an individual function... Something like g' = next g → C g' ∈ R (P g'), but I don't like that, and it wouldn't let me prove that it holds for every RakeSieve object, so RakeSieve still couldn't be a PrimeSieve.

So instead, the proposition has to be attached to the structure itself, and each individual function has to supply its own proof.

I already knew how to do this, and it's obvious to me in retrospect, but I still had the mistake sitting in the stub file this whole time, and it was still surprising to me when I sat down to prove it. :)

[Michal Wallace (tangentstorm)]

Well, even that's not true. I still need the standalone theorem to prove instance : PrimeSieve RakeSieve... But I can't prove it unless there's a corresponding proof on the structure. And C and P aren't even part of the structure... Those are defined as part of SieveState.

instance : SieveState RakeSieve where
  P x := x.p
  C x := x.c
  next x hC := RakeSieve.next x hC
  hNext := by aesop

It seems like the general rule is that if you want ensure that a property holds in a class, you have to also prove a very similar thing for the structure itself. (But in this case, instead of writing it in terms of C x and P x, I have to write it in terms of x.c and x.p... and then hopefully derive the same facts for C x and P x...

[Michal Wallace (tangentstorm)]

Live-proving again tonight, starting in a few minutes.
https://youtube.com/live/jQ_b5-y623o?feature=share

[Michal Wallace (tangentstorm)]

almost there.

In the latest version of RakeSieve includes proofs for all of the following properties:

structure RakeSieve where
  prop : Nat -> Prop
  rm: RakeMap prop
  p : NPrime               -- the current prime
  c : Nat                  -- the candidate for next prime
  hprop : ∀ n:Nat, prop n ↔ n ∈ R p
  hCinR : c ∈ R p
  hRmin : ∀ r ∈ R p, c ≤ r

... except it doesn't yet prove hRmin for the next operation. That proof, and everything else in the whole system, is now waiting on me to verify the properties of RakeMap itself, which is just a list of arithmetic sequences that doesn't need to know anything about primes.

The remaining issues until the proof is complete:

• demonstrate that term 0 of a Rakemap is the smallest term
• deal with the possibility that rem (removing multiples of a number) could remove all the sequences, which means that the set is empty and term can't return a value (therefore, rem has to return Option RakeMap instead of RakeMap)
• fill in the six lemmas to verifygte and rem accurately model the set operations they're supposed to reflect.
• write the (hopefully tiny amount of) code to show instance : PrimeSieve PrimeGen

... And then it's done. Still kind of a long list, but it at least feels like it's all downhill from here. :)

(I still want to add some laziness to make it not so embarrassingly slow, but that'll have to wait until I've made enough progress on the SoME video.)

[Michal Wallace (tangentstorm)]

(still) almost there

The only proofs left to do are in Rake itself. I'm currently doing some refactoring work to eliminate the DSeq type, which was just a subtype of ASeq where the d was baked into the type (so I could assert that all the sequences in the Rake had the same delta). Instead, the new version of Rake just stores a list of Nats (ks) for the different constant terms, and has only one copy of d, and just instantiates aseq r.ks[i] r.d whenever it needs the sequence.

I've also discovered I have a problem that will not actually arise in a prime sieve, but that complicates my type: basically, a rake is a list of parallel arithmetic sequences, and if that list were empty, then it would be meaningless to talk about the nth term. But one of my core operations is to remove multiples of a number.

Since the number I'm passing in is always prime, and conditions in the prime sieve are set up so that each arithmetic sequence I'm dealing with always contains numbers coprime to it, I will never actually create an empty list. But proving all that would be a lot of work, so if I don't prove it, lean won't know for certain that the list of sequences is non-empty. That makes me want to return an Option Rake... But that's a pain. Instead, I think what I will do is basically what Nat does for 2 - 5... return something like a 0.

So, I'm thinking I allow d=0, and then I can create rakes that return arbitrary repeating patterns of constants... And in particular, the rake { d := 0, ks := [0] } would just give you 0 for every term. So, rather than Option Rake, I'll just return this predefined rake0... I probably will have to pass in another hypothesis or two to the rem operation to get away with this, but I'm cross that bridge when I come to it.

[Michal Wallace (tangentstorm)]

Rake invariants

I'm down to the last 6 lemmas... Three each to say that the gte and rem operations remove, retain, and don't add extra elements of the set...

They're all looking like they're going to be very unpleasant. Why? Because for both operations, I'm first calling out to a lower level operation on arithmetic sequences in ASeq, and then after the important work is done, I'm doing a bunch of extra work to make sure the sequences within the rake are sorted.

But the thing is, it's easy to sort the sequences at any time, so I think instead of forcing the thing to always be sorted, I'm going to add a Boolean field that says whether it's sorted, and then the proof field will contain an implication: sorted → List.Sorted (·<·) ks...

That would also allow rakes to produce a wider variety of sequences. (For example, with d=0 it could produce arbitrary cyclic patterns).

The reason for all this is that the proofs seem to require me to break the pipeline down into small steps and prove that the same properties hold at each step. You can see it in the implementation of gte:

def Rake.gte (r: Rake) (n: Nat) : Rake  :=
  let f : ℕ → ℕ := (λk => let s := aseq k r.d; (ASeq.gte s n).k)
  let ks₀ := r.ks.map f
  let ks₁ := ks₀.dedup
  let ks₂ := sort_nodup ks₁ ks₀.nodup_dedup
  have hsize : 0 < ks₂.val.length := by
    have h₀ : 0 < ks₀.length := Nat.lt_of_lt_of_eq r.hsize (r.ks.length_map _).symm
    have h₁ : 0 < ks₁.length := length_pos_of_dedup h₀
    have h₂: ks₂.val.length = ks₁.length := ks₁.length_mergeSort _
    exact Nat.lt_of_lt_of_eq h₁ h₂.symm
  { d := r.d, ks := ks₂, hsort := ks₂.prop.left, hsize := hsize }

If I don't break this apart, I don't see how I can avoid duplicating that same kind of logic for each of the 6 proofs. (And possibly for 3 more proofs that the partition step called by rem works).

The other alternative I was considering was writing a macro or tactic, but so far I'm having a hard time just writing one proof by hand to see what it is I'd need to generate...

[Michal Wallace (tangentstorm)]

update : making that change simplified things tremendously.
https://github.com/tangentstorm/treesiv/commit/205aef38d7bb377d2ebcc9bc8322c471e767a219?diff=split&w=0

[Michal Wallace (tangentstorm)]

Here's how that same function looks now:

def Rake.gte (r: Rake) (n: Nat) : Rake  :=
  let f : ℕ → ℕ := (λk => let s := aseq k r.d; (ASeq.gte s n).k)
  { d := r.d, ks := r.ks.map f, sorted := false,
    hsize := Nat.lt_of_lt_of_eq r.hsize (r.ks.length_map _).symm }

[Michal Wallace (tangentstorm)]

the other 90%....

I thought these last few proofs would be easy, but I've had a hard time even getting started on them.

I think I may actually need 8 proofs: 2 for the partition operation, to show that every generated term in the old rake is in the new rake and vice versa.... Whereas the gte and rem operations have to show:

• every term in the new rake was in the old rake
• every term in the old rake that we want to keep is in the new rake
• every term in the old rake that we want to discard got discarded.

rem calls partition under the hood, so i imagine I'll need to verify it along the way.

But how do I make the proofs? The statements I wrote weeks ago all look like this:

variable (rm: RakeMap prop) (p n: Nat)

lemma RakeMap.gte_drop -- gte drops terms < p
    : n < p → ¬(∃m, (rm.rake.gte p).sort.val.term m = n) := sorry

lemma RakeMap.gte_keep -- gte keeps terms ≥ p
  : n≥p ∧ (∃pm, rm.rake.term pm = n) → (∃m, (rm.rake.gte p).sort.val.term m = n) := sorry

lemma RakeMap.gte_same -- gte introduces no new terms
  : (∃m, (rm.rake.gte p).sort.val.term m = n) → (∃pm, rm.rake.term pm = n) := by

For the keep and same lemmas, it seems like I have to demonstrate that a particular input produced a particular output before and after the transformation.

In my other Rake proofs, I've completely unfolded Rake.term and do a bunch of modular arithmetic to map a specific input to a specific position in the sequence... But I think I can avoid doing this for every proof if I just had a theorem that says it's sufficient to show that one of the constants in rake.ks plus some m * rake.d produces the desired value. In other words, it shouldn't really matter which term it is, or where it is in the sequence.

For the proof of gte_drop, I guess I need to show that the minimum term is greater than or equal to the given number.

For rem_drop... I think I have to show that partition can produce sequences that either always produce or never produce multiples of the given prime, and that those are the ones I remove.

It all seems easy enough in my head, but I'm still struggling a bit to express this in lean...

[Michal Wallace (tangentstorm)]

.... It occurs to me that gte is only used once at the moment, to remove 0 and 1 (since these can't be prime). In the interest of time, I think I'll just cut it out for now, and manually prove that { d := 1, ks := [2] } corresponds to { n : Nat // n ≥ 2 }.

[Michal Wallace (tangentstorm)]

Rake.term_iff

if I just had a theorem that says it's sufficient to show that one of the constants in rake.ks plus some m * rake.d produces the desired value

Welp. It took me 5 hours, but I proved it (Rake.lean)

theorem Rake.term_iff (r:Rake)
  : ∀m, (∃n, r.term n = m) ↔ (∃ k∈r.ks,  ∃x:Nat, k+x*r.d = m) := sorry

[Michal Wallace (tangentstorm)]

I live-streamed that session, btw. I don't really expect anyone to want to watch me flounder around, but on the off chance you do, there's now six of these videos from the past month in the playlist is here:

https://www.youtube.com/playlist?list=PLwI69hwjdBsY-hn-JTA7fG-v3gWWVmTxq
hqdefault.webp

[Michal Wallace (tangentstorm)]

I think the effort I put into term_iff is paying off. I was able to prove that sort preserves the set of terms the rake produces by applying the exact same argument in both directions simultaneously using a long chain of tactics under all_goals. The only differences at the end were the direction a few rewrite rules were applied.

theorem sort_term_iff_term  (r:Rake) (n:Nat)
  : (∃m, r.term m = n) ↔ (∃m', r.sort.val.term m' = n) := by
  if h: r.sorted then unfold sort; aesop
  else
    let ⟨rs, ⟨hd, hsort, hperm⟩⟩ := r.sort; simp
    have : ∀k, k∈rs.ks ↔ k∈r.ks.dedup := fun k => List.Perm.mem_iff hperm
    have hmem : ∀k, k∈r.ks ↔ k∈rs.ks := by aesop
    apply Iff.intro
    all_goals
      intro h; rw[term_iff] at h; obtain ⟨k, hk₀, hk₁⟩ := h
      specialize hmem k
      rw[term_iff]
      use k
      apply And.intro
    · rwa[←hmem]
    · rwa[hd]
    · rwa[hmem]
    · rwa[←hd]

(commit fd5c45d)

[Michal Wallace (tangentstorm)]

even more "asmost there"

I've made some pretty decent progress, after a long night of fighting with natural number division and temporarily giving up on my proof that partition retains the same set of terms, I moved on to proving the properties of the rem operation (removing multiples of a number).

I proved all 3 properties of rem, that I needed, but there's a problem: I have to do something if a rake contains nothing but multiples of x and you tell it to remove all multiples of x... I could have rem return Option Rake but I knew that would cause a ton of headache down the line, so instead I just have it return a zer that has the formula 0 + 0*n (so an infinite stream of zeros).

Unfortunately, this also causes headaches for the caller, so I've realized what I really need to pass in is a proof that the Rake can produce numbers that are not divisible by X.

In fact, Dirichlet's prime number theorem tells us that any arithmetic sequence k + d *n where k and d are co-prime will yield an infinite number of primes, and there is an elementary proof of this due to Atle Selberg that I think might be approachable in Lean.... At some point.

But I don't think I'll need anything like that for my purposes, since I've already got a prime sieve, and I already know that the set of "numbers I haven't crossed off yet" contain all the primes. In fact, term 0 is the next prime, and multiples of all smaller primes are already crossed out, so the next compound number must be (term 0)^2... and there must be a prime in between because a prime always exists between n and 2*n... a fact which is already in mathlib...

All together, that means either term 1 is also a new prime. (I might have to prove that the nth prime is always smaller than the nth primorial, which seems obvious but might be hard (???).... But even if I can't prove that, the next smallest term after term 0 would be term (ks.length) so either way it should be good.)

So, in short, I now need to:

• finish the proof that partition leaves the terms unchanged
• patch the rem proofs to eliminate the zer case using proof of a second prime in the rake
• create and supply that proof

[Michal Wallace (tangentstorm)]

QED