Working with Nat.factorization

Stream: mathlib4

[Eric Wieser]

In lean 4.7.0, this proof worked just fine:

import Mathlib

#eval (Nat.factorization 10485760000000000000000000) 5

example : (Nat.factorization 10485760000000000000000000) 5 = 19 := by
  rfl -- fails in latest lean

but since docs#Nat.find is now (transitively) irreducible, it no longer goes through. What's the preferred way to prove things like the statement above? Should Nat.find be globally overriden to be semireducible?

[Yaël Dillies]

@Kyle Miller suggested a reimplementation of Nat.find in the past

[Eric Wieser]

here I assume

[Eric Wieser]

In the meantime, what's the preferred way to prove the above?

[Matthew Ballard]

Eric Wieser said:

since docs#Nat.find is now (transitively) irreducible, it no longer goes through.

Shouldn't this mean that with_unfolding_all rfl goes through?

[Eric Wieser]

This works, but seems pretty silly:

theorem Nat.factorization_list_prod {l : List ℕ} (hS : ∀ x ∈ l, x ≠ 0) :
    l.prod.factorization = (l.map Nat.factorization).sum := by
  induction l with
  | nil => simp
  | cons head tail ih =>
    rw [List.prod_cons, List.map_cons, Nat.factorization_mul, List.sum_cons, ih]
    · aesop
    · aesop
    · aesop


example : (Nat.factorization 10485760000000000000000000) 5 = 19 := by
  change (Nat.factorization (eval% Nat.primeFactorsList 10485760000000000000000000).prod) 5 = 19
  rw [Nat.factorization_list_prod]
  · simp [List.map_cons, Nat.Prime.factorization, Nat.prime_two, Nat.prime_five, List.sum_cons]
  · decide

[Michael Stoll]

It does not look like docs#Nat.find is irreducible. Trying

import Mathlib

#eval (Nat.factorization 10485760000000000000000000) 5

unseal Nat.find in
example : (Nat.factorization 10485760000000000000000000) 5 = 19 := by
  rfl -- fails in latest lean

gives the error

failed to set `[local semireducible]`, `Nat.find` is currently `[semireducible]`, `[irreducible]` expected
use `set_option allowUnsafeReducibility true` to override reducibility status validation

and indeed, there is no @[irreducible] in the source.

[Ruben Van de Velde]

There was a recent change where definitions based on well-founded recursion became irreducible without being tagged explicitly (if my memory serves me right)

[Matthew Ballard]

Ruben Van de Velde said:

There was a recent change where definitions based on well-founded recursion became irreducible without being tagged explicitly (if my memory serves me right)

Does this break with_unfolding_all?

[Eric Wieser]

I think the issue here is that well-founded recursion is used?

[Michael Stoll]

But why does Lean claim that `Nat.find` is currently `[semireducible]` ?

[Michael Stoll]

Same with docs#Nat.findX, which is the actual definition using well-founded recursion.

[Eric Wieser]

They're all semireducible, but wrap Acc.rec which isn't

[Michael Stoll]

...and apparently one cannot make Acc.rec (semi)reducible; at least unseal Acc.rec in says

failed to set reducibility status, `Acc.rec` is not a definition

[Matthew Ballard]

#15409

[Matthew Ballard]

I didn't test whether Eric's example now works

[Matthew Ballard]

Matthew Ballard said:

I didn't test whether Eric's example now works

It doesn't

[Eric Wieser]

I suspect the path forward here includes writing a norm_num extension for padicValNat and factorization; I tried to make a start with

@[norm_num DFunLike.coe (Nat.factorization _) _]
def evalNatFactorization : NormNumExt where eval {u α} e := do
  match u, α, e with
  | 0, ~q(ℕ), ~q(Nat.factorization $m $p) =>
    trace[debug] "Ok"
    let ⟨mlit, hm⟩ ← deriveNat m q(inferInstance)
    let ⟨plit, hp⟩ ← deriveNat p q(inferInstance)
    trace[debug] m!"Testing primality of {p}"
    match ← deriveBool q(Prime $plit) with
    | ⟨true, (hp : Q(Prime $plit))⟩ =>
      trace[debug] "Found a prime"
      -- let nlit : Q(ℕ) := mkRawNatLit <| padicValNat plit.natLit! mlit.natLit!
      failure
    | ⟨false, (hp : Q(¬ Prime $plit))⟩ =>
      trace[debug] "Not a prime"
      return .isNat q(inferInstance) q(nat_lit 0) q(⟨factorization_eq_zero_of_non_prime _ $hp⟩)
  | _, _, _ => failure

set_option trace.debug true

-- test-cases

example : Nat.factorization 24 4 = 0 := by
  norm_num

example (x) (hx : x = 3) : Nat.factorization 24 2 = 3 := by
  norm_num
  rw [hx]

but deriveBool doesn't seem to work

[Matthew Ballard]

Eric Wieser said:

In lean 4.7.0, this proof worked just fine:

Is there specific commit where this works because I went all the way back to v4.2.0 and nothing worked

[Eric Wieser]

Hmm, perhaps I minimized badly here :(.

[Matthew Ballard]

set_option maxRecDepth 10000 in
example : (Nat.factorization 10485760000000000000000000) 5 = 19 := by
  simp (config := { ground := true })

might do it if can finish. It can at least do smaller numbers.

[Matthew Ballard]

example : (Nat.factorization 10485760000000000000000000) 5 = 19 := by
  have : Nat.Prime 5 := by decide
  unfold Nat.factorization; simp
  split_ifs
  unfold padicValNat; simp
  unfold multiplicity; simp
  unfold Nat.find; simp
  unfold Nat.findX; simp
  rw [WellFounded.fix_eq]; dsimp; split_ifs with h'; swap; exact h' (by decide)
  iterate 18 (rw [WellFounded.fix_eq]; dsimp; split_ifs with h'; swap; simp; exact h' (by decide))
  rw [WellFounded.fix_eq]; dsimp; split_ifs with h'; swap; rfl
  exact False.elim <| (by decide : ¬ _ ∣ _) h'

is unpleasant but works.

[Michael Stoll]

Would it help to have a version of docs#Nat.find that takes an explicit upper bound? E.g.

def Nat.find' (P : Nat → Prop) [DecidablePred P] : Nat → Nat
| 0 => 0
| n + 1 => let r := find' P n
           if r < n then r else
           if P n then n
           else n + 1

such that Nat.find P bound is the smallest n < bound that satisfies P if such an n exists, and is bound otherwise.

I would think that in many uses of Nat.find, an explicit upper bound is available.

[Matthew Ballard]

Probably not I take that back. The following is a more reasonable toy model and works ok

def foo (n : ℕ) (hn : n ≤ 100) : ℕ :=
  go n 0 hn
where
  go (a b : ℕ) (ha : a + b ≤ 100) :=
  if h : a = 0 then b
  else go (a - 1) (b + 1) (by omega)
  termination_by 100 - b
  decreasing_by
    simp_wf
    omega

#eval foo 5 (by decide)

example : foo 5 (by decide) = 5 := by
  rfl

[Michael Stoll]

def Nat.find'_aux (P : Nat → Prop) [DecidablePred P] (a n : Nat) : Nat → Nat
| 0 => n
| k + 1 => if P a then a else find'_aux P (a + 1) n k

def Nat.find' (P : Nat → Prop) [DecidablePred P] (n : Nat) : Nat :=
  find'_aux P 0 n n

def Nat.mult_aux (k m n : Nat) : Prop := k ^ m ∣ n ∧ ¬ k ^ (m + 1) ∣ n

instance : ∀ k n : Nat, DecidablePred (Nat.mult_aux k · n) :=
  fun _ _ _ => instDecidableAnd

def Nat.multiplicity' (m n : Nat) : Nat :=
  if n = 0 then 0 else Nat.find' (mult_aux m · n) n

example : Nat.multiplicity' 5 10485760000000000000000000 = 19 := rfl

(using 187 heartbeats)

[Michael Stoll]

Here is a more fleshed-out (and somewhat improved) version that includes a proof that Nat.multiplicity' agrees with docs#multiplicity:

import Mathlib.RingTheory.Multiplicity

/-- A missing API lemma for `PartENat.find`. -/
lemma PartENat.find_eq_natFind {P : Nat → Prop} [DecidablePred P] (h : ∃ n, P n) :
    PartENat.find P = Nat.find h := by
  apply Part.ext'
  · simp only [find_dom _ h, dom_natCast]
  · simp only [dom_natCast, find_get, get_natCast', imp_self, implies_true]


namespace Nat

/-- Helper function for `Nat.findWithBound`. -/
@[simp]
def findWithBound_aux (P : Nat → Prop) [DecidablePred P] (n : Nat) : Nat → Nat
| 0 => n
| k + 1 => if P n then n else findWithBound_aux P (n + 1) k

/-- The defining property of `Nat.findWithBound_aux`. -/
lemma findWithBound_aux_spec (P : Nat → Prop) [DecidablePred P] (n k m : Nat) :
    findWithBound_aux P n k = m ↔
      n ≤ m ∧ m ≤ k + n  ∧ (∀ l, n ≤ l → l < m → ¬ P l) ∧ (P m ∨ m = k + n) := by
  induction k generalizing n with
  | zero =>
      simp only [findWithBound_aux, zero_add]
      refine ⟨fun H ↦ ?_, fun ⟨H₁, H₂, _⟩ ↦ le_antisymm H₁ H₂⟩
      simpa only [H, le_refl, zero_add, or_true, and_true, true_and]
        using fun _ hl₁ hl₂ _ ↦ (hl₂.trans_le hl₁).false
  | succ k ih =>
      simp only [findWithBound_aux]
      by_cases hP : P n <;> simp only [hP, ↓reduceIte]
      · refine ⟨fun H ↦ ?_, fun ⟨H₁, _, H₃, _⟩ ↦ ?_⟩
        · simpa only [H, le_refl, le_add_iff_nonneg_left, _root_.zero_le, self_eq_add_left,
            add_eq_zero, one_ne_zero, and_false, or_false, true_and]
            using ⟨fun _ hl₁ hl₂ _ ↦ (hl₂.trans_le hl₁).false, H ▸ hP⟩
        · exact le_antisymm H₁ <| le_of_not_lt <| (imp_not_comm.mp <| H₃ n le_rfl) hP
      · specialize ih (n + 1)
        rw [show k + (n + 1) = k + 1 + n by rw [add_comm n, add_assoc]] at ih
        refine ⟨fun H ↦ ?_, fun ⟨H₁, H₂, H₃, H₄⟩ ↦ ?_⟩
        · obtain ⟨ih₁, ih₂, ih₃, ih₄⟩ := ih.mp H
          refine ⟨(le_add_right n 1).trans ih₁, ih₂, fun l hl₁ hl₂ ↦ ?_, ih₄⟩
          exact (le_iff_lt_or_eq.mp hl₁).elim (fun h ↦ ih₃ l h hl₂) fun h ↦ h ▸ hP
        · have H : n ≠ m := by
            rintro rfl
            exact (zero_ne_add_one k).symm <| self_eq_add_left.mp <| H₄.resolve_left hP
          exact ih.mpr ⟨lt_of_le_of_ne H₁ H, H₂, fun l h ↦ H₃ l ((le_add_right n 1).trans h), H₄⟩


/-- `Nat.findWithBound P n` returns the smallest `m ≤ n` such that `P m` holds if such an `m`
exists; otherwise, `Nat.findWithBound P n = n`. -/
def findWithBound (P : Nat → Prop) [DecidablePred P] (n : Nat) : Nat :=
  findWithBound_aux P 0 n

/-- The defining property of `Nat.findWithBound`. -/
lemma findWithBound_eq_iff (P : Nat → Prop) [DecidablePred P] (n m : ℕ) :
    findWithBound P n = m ↔ m ≤ n ∧ (∀ k < m, ¬ P k) ∧ (P m ∨ m = n) := by
  simpa only [_root_.zero_le, add_zero, true_implies, true_and] using findWithBound_aux_spec P 0 ..

lemma findWithBound_spec (P : Nat → Prop) [DecidablePred P] (n : Nat) :
    P (findWithBound P n) ∨ findWithBound P n = n :=
  ((findWithBound_eq_iff P n _).mp rfl).2.2

lemma findWithBound_le' (P : Nat → Prop) [DecidablePred P] (n : Nat) :
    findWithBound P n ≤ n :=
  ((findWithBound_eq_iff P n _).mp rfl).1

lemma findWithBound_le {P : Nat → Prop} [DecidablePred P] (n : Nat) {m : Nat} (h : P m) :
    findWithBound P n ≤ m :=
  le_of_not_lt <| (imp_not_comm.mp <| ((findWithBound_eq_iff P n _).mp rfl).2.1 m) h

lemma findWithBound_min {P : Nat → Prop} [DecidablePred P] {n m : Nat}
    (hm : m < findWithBound P n) :
    ¬P m :=
  ((findWithBound_eq_iff P n _).mp rfl).2.1 m hm

/-- `Nat.findWithBound P n` agrees with `Nat.find (_ : ∃ m, P m)` if `P m` holds
for some `m ≤ n`. -/
lemma findWithBound_eq_find_of_witness {P : Nat → Prop} [DecidablePred P] {m n : Nat} (h : m ≤ n)
    (hP : P m) :
    findWithBound P n = Nat.find ⟨m, hP⟩ := by
  rw [findWithBound_eq_iff]
  refine ⟨(find_le_iff ⟨m, hP⟩ n).mpr ⟨m, h, hP⟩, ?_, .inl <| Nat.find_spec ⟨m, hP⟩⟩
  exact fun k hk ↦ (lt_find_iff ⟨m, hP⟩ k).mp hk k le_rfl

/-- A constructive version of `multiplicity` for natural numbers. -/
def multiplicity' (m n : Nat) : PartENat :=
  if n = 0 ∨ m = 1 then ⊤ else findWithBound (fun k ↦ ¬ m ^ (k + 1) ∣ n) n

/-- `Nat.multiplicity' m n` agrees with `multiplicity m n` for `m n : Nat`. -/
lemma multiplicity'_eq_multiplicity (m n : Nat) : multiplicity' m n = multiplicity m n := by
  rcases eq_or_ne n 0 with rfl | hn
  · simp only [multiplicity', true_or, ↓reduceIte, multiplicity.zero]
  rcases eq_or_ne m 1 with rfl | hm
  · simp only [multiplicity', or_true, ↓reduceIte, isUnit_one, multiplicity.isUnit_left]
  simp only [multiplicity', hn, hm, or_self, ↓reduceIte]
  have key : ¬ m ^ (n + 1) ∣ n := by
    intro H
    rcases eq_or_ne m 0 with rfl | hm₀
    · exact hn <| zero_dvd_iff.mp <| (zero_pow (M₀ := Nat) <| succ_ne_zero n) ▸ H
    · exact lt_irrefl n <|
        calc n
          _ < n + 1 := lt_add_one n
          _ < m ^ (n + 1) := lt_pow_self (by omega) _
          _ ≤ n := le_of_dvd (zero_lt_of_ne_zero hn) H
  rw [findWithBound_eq_find_of_witness le_rfl key]
  exact (PartENat.find_eq_natFind ⟨_, key⟩).symm

count_heartbeats in -- 280
example : multiplicity 5 10485760000000000000000000 = 19 := by
  rw [← multiplicity'_eq_multiplicity]
  rfl

[Michael Stoll]

Some questions that came to my mind:

• Would something along the lines of Nat.findWithBound be a reasonable addition to Mathlib?
• And if so, should docs#padicValNat (and possibly similar and maybe further definitions) be refactored to use it?
• What is the rationale for having both docs#ENat and docs#PartENat ?

[Michael Stoll]

@Eric Wieser @Matthew Ballard @Kyle Miller :up: ?

[Matthew Ballard]

I have been ruminating