Compute Nat.primeFactors

Stream: Is there code for X?

[Antoine Chambert-Loir]

I'm struggling to make Lean compute docs#Nat.primeFactors , even for small integers.

example : Nat.primeFactors 6 = {2, 3} := sorry

I can rewrite to docs#Nat.primeFactorsList, and have it apply the induction hypothesis step by step, but this is a bit painful, even for 6.

[Antoine Chambert-Loir]

Thinking about repeat suggested this:

import Mathlib.Data.Nat.Factors
example : Nat.primeFactorsList 6 = [2, 3] := by
    repeat
      change Nat.primeFactorsList (_ + 2) = _
      simp [Nat.primeFactorsList]
      norm_num

Is it the right way to do?

[Antoine Chambert-Loir]

Probably not — What if I don't want to compute myself the right hand side ?

[Aaron Liu]

I think this is an instance of #lean4 > opaque recursion definitions break mergeSort decidability

[Eric Wieser]

#15420

[Antoine Chambert-Loir]

Cool! (I had started writing a fueled version of Nat.primeFactors…)

[Eric Wieser]

That is probably still a good idea

[Eric Wieser]

If you do that the simproc can become a dsimproc

[Antoine Chambert-Loir]

Is it important that the fuel variable is initialized to something reasonable?
Finding all prime factors of n requires Nat.clog 2 n steps; I I have two versions, one that proves that it goes in that much steps, and another one where I initialize fuel as n.
But the first version doesn't work as rfl, while the second does…

[Antoine Chambert-Loir]

import Mathlib

theorem primeFactorsList2_aux (fuel : ℕ) (n : ℕ) (hn : n ≤ fuel + 1) (m : ℕ) (hm : 2 ≤ m) :
    n / m ≤ fuel := by
  apply Nat.le_of_lt_succ
  rw [Nat.div_lt_iff_lt_mul]
  apply lt_of_le_of_lt hn
  change fuel.succ < fuel.succ * m
  rw [Nat.lt_mul_iff_one_lt_right]
  exact hm
  exact Nat.zero_lt_succ fuel
  exact Nat.zero_lt_of_lt hm

def primeFactorsList3 : ℕ → List ℕ :=
  let rec go (fuel : ℕ) (n : ℕ) (hfuel : n ≤ fuel) : List ℕ :=
    match fuel with
    | 0 => []
    | fuel + 1 =>
      if hk : fuel = 0 ∨ n < 2 then []
      else
        let m := Nat.minFac (n)
        have hm : 2 ≤ m := by
            apply Nat.Prime.two_le
            simp [m]
            refine Nat.minFac_prime ?_
            simp only [not_or] at hk
            intro hn
            apply hk.2
            simp [hn]
        m :: (go fuel (n / m) (primeFactorsList2_aux _ _ hfuel _ hm))
  fun n ↦ go n n (le_refl _)

theorem primeFactorsList4_aux (fuel : ℕ) (n : ℕ) (h : ¬ (fuel = 0 ∨ n < 2)) (hn : n ≤ 2 ^ (fuel + 1)) (m : ℕ) (hm : 2 ≤ m) :
    n / m ≤ 2 ^ fuel := by
  apply Nat.le_of_lt_succ
  rw [Nat.div_lt_iff_lt_mul (Nat.zero_lt_of_lt hm)]
  apply lt_of_le_of_lt hn
  apply lt_of_lt_of_le (b := (2 ^ fuel).succ * 2)
  · simp only [Nat.succ_eq_add_one, add_mul, one_mul,
    ← pow_succ, Nat.lt_succ_iff, Nat.le_add_right]
  · exact Nat.mul_le_mul_left (2 ^ fuel).succ hm

def primeFactorsList4 : ℕ → List ℕ :=
  let rec go (fuel : ℕ) (n : ℕ) (hfuel : n ≤ 2 ^ fuel) : List ℕ :=
    match fuel with
    | 0 => []
    | fuel + 1 =>
      if hk : fuel = 0 ∨ n < 2 then []
      else
        let m := Nat.minFac (n)
        have hm : 2 ≤ m := by
          apply Nat.Prime.two_le
          simp [m]
          refine Nat.minFac_prime ?_
          simp only [not_or] at hk
          intro hn
          apply hk.2
          simp [hn]
        m :: (go fuel (n / m) (primeFactorsList4_aux fuel n hk hfuel _ hm))
  fun n ↦ go (Nat.clog 2 n) n (Nat.le_pow_clog Nat.one_lt_two n)


example : Nat.primeFactorsList3 105 = [3, 5, 7] := rfl
example : Nat.primeFactorsList4 105 = [3, 5, 7] := rfl -- doesn't work

[Eric Wieser]

I think the problem is that clog 2 n is itself missing a fuel argument?

[Antoine Chambert-Loir]

(Note that I do not understand the what in this way of writing a proof mkes this work.)

[Eric Wieser]

Indeed, example : Nat.clog 2 37 = (eval% Nat.clog 2 37) := rfl fails

[Antoine Chambert-Loir]

OK, so I have to fuel docs#Nat.clog as well. This makes sense.

[Eric Wieser]

I think you want your fuel argument to be as simple as possible, because how this method works is by reducing fuel to applications of succ

[Eric Wieser]

So we should fuel clog, but I doubt you want to use it in prime factors

[Antoine Chambert-Loir]

Hence the first version will be the right one.

[Eric Wieser]

I expect so, but perhaps the right approach is:

• Fuel clog
• Benchmark both versions, to see which is faster

[Antoine Chambert-Loir]

OK, I'll do that. (This is fun.)

[Antoine Chambert-Loir]

How do you benchmark something?
I tested both versions with big numbers, and crashed Lean.

[Eric Wieser]

set_option trace.profiler true in is a reasonable way

[Eric Wieser]

Pick numbers small enough not to crash, but big enough to take more than a second

[Antoine Chambert-Loir]

I wonder whether docs#Nat.minFac shouldn't be fueled as well. (I get fast results with #eval but Lean pretends the recursion is too high — there are less than 5 prime number.)

[Bhavik Mehta]

The change here would be to docs#Nat.minFacAux, which is the one containing the recursion. The intention is, as the docstring there describes, for norm_num or related to be the one computing this during proves. But (as I guess you're encountering just like I have), that's not always enough when the rest of a computation works by decide

[Bhavik Mehta]

As an aside, Nat.sqrt is also something I'd marked semireducible for the same reason, and that one should also be fueled

[Eric Wieser]

One thing to check here is that the weight of the extra fuel doesn't slow down #eval

[Eric Wieser]

If it does, then there is some extra work of writing a csimp lemma and keeping both versions

[Antoine Chambert-Loir]

I have a question about the docstring for docs#Nat.logC which says

An alternate definition for Nat.log which computes more efficiently. For mathematical purposes, use Nat.log instead, and see Nat.log_eq_logC.

Given function extensionality, I don't understand what is really meant here. Maybe, that it is easier to prove a mathematical lemma using the inductive definition of docs#Nat.log than that of docs#Nat.logC.

[Johan Commelin]

I think that's indeed the purpose of that comment. Some definitions are good for proofs, others are good for computations.

[Antoine Chambert-Loir]

i've made a preliminary fueled definition to replace docs#Nat.clog, which takes roughly the same time as the computation using docs#Nat.logC.
However, for large numbers, one needs to increase maxRecDepth to allow rfl to prove the desired equality:

set_option maxRecDepth 900 in
example : clog_fueled 2 123456789012345667890123456789012345667890123456789012345667890123456789012345667890 = 277 := rfl

This is roughly three times the size of the output, as if the tests that the iteration must perform were only done at the end. I had tried to make my recursion as much terminal as possible, but that is apparently not enough. Any idea?

[Johan Commelin]

@Antoine Chambert-Loir Is your fueled version on a branch somewhere, so that others can play with it?

[Antoine Chambert-Loir]

For the moment, it's just a plain file with iterative attempts. I'll make that at once. Should I put this in Mathlib.Data.Nat.Log or open a new file?

[Johan Commelin]

New file seems fine for experiments. We can merge it into the existing files later.

[Antoine Chambert-Loir]

#25557
branch#ACL/fuel_clog

[Johan Commelin]

@Antoine Chambert-Loir What do you think of

def clog' (b n : ℕ) : ℕ :=
  let rec go (fuel : ℕ) (b n : ℕ) : ℕ :=
  match fuel with
  | 0 => 0
  | fuel + 1 => if 1 < b ∧ 1 < n then go fuel b ((n + b - 1) / b) + 1 else 0
  go n b n

#reduce clog' 2 37 -- 6

[Antoine Chambert-Loir]

This was basically my first attempt. For large numbers, rfl doesn't work, it requires maxRecDepth to be roughly 10 times the size of the output.

[Antoine Chambert-Loir]

This is why I rewrote is using an accumulator (to avoid the final + 1).

[Johan Commelin]

I guess in the meta world it will be easy to make this much faster, because then we have access to the binary representation of n.

[Antoine Chambert-Loir]

I just noticed something bizarre: with my version, #reduce requires the same recursion depth as yours, but proving equality (example : ... = ... := rfl) requires less.

[Mario Carneiro]

This is the best version I could come up with in terms of keeping the recursion depth to clog n + O(1):

def clog_fueled (b n : ℕ) : ℕ :=
  if b ≤ 1 then 0 else
    n.rec (fun a _ => a) (fun _ ih a n =>
      (Nat.ble n 1).rec (ih (Nat.succ a) ((n + b - 1) / b)) a) 0 n

set_option maxRecDepth 286 in
#reduce clog_fueled 2
  123456789012345667890123456789012345667890123456789012345667890123456789012345667890

[Mario Carneiro]

note that this is not the kernel checker doing the reduction though

[Eric Wieser]

Antoine Chambert-Loir said:

I just noticed something bizarre: with my version, #reduce requires the same recursion depth as yours, but proving equality (example : ... = ... := rfl) requires less.

I would guess that reduce always expends all the fuel, whereas rfl can stop sooner?

[Mario Carneiro]

no, the recursion depth is absurdly low, I think it's either somehow tail recursive or the recursion isn't being tracked

[Mario Carneiro]

kernel evaluation here is weird, I don't think there is a way to avoid proper recursion

[Eric Wieser]

Does swapping the rec for a match cause trouble?

[Mario Carneiro]

for the bool or the nat?

[Eric Wieser]

The nat

[Mario Carneiro]

for the nat it's a disaster, it triggers equation compilation and you get a brecOn and a bunch of cruft

[Eric Wieser]

even with termination_by structural?

[Mario Carneiro]

that is structural recursion as lean understands it

[Eric Wieser]

oh

[Mario Carneiro]

Eric Wieser said:

Antoine Chambert-Loir said:

I just noticed something bizarre: with my version, #reduce requires the same recursion depth as yours, but proving equality (example : ... = ... := rfl) requires less.

I would guess that reduce always expends all the fuel, whereas rfl can stop sooner?

note that if you actually expended all the fuel it would not terminate in your lifetime

[Mario Carneiro]

since the fuel is ~2^277

[Antoine Chambert-Loir]

True, but you know a priori that you don't need that much fuel, just the log of it.

[Mario Carneiro]

I was thinking about that, but unfortunately it's not easy to exploit that, since you... need to compute the log in order to use that as the fuel instead

[Mario Carneiro]

it is a bit of a shame that the computation involves constructing the numbers 123...7890, 123...7889, 123...7888 and so on down to n-277 since these numbers are pointlessly big. If it was a megabytes large number this would be quite expensive when actually computing the base 2 log of the number would be trivial on a normal computational system

[Mario Carneiro]

and even computing the general base log of a huge number can slice the number down very quickly, resulting in a constant factor overhead over the input storage

[Johan Commelin]

But don't we have the binary representation of the number at hand, whenever we actually want to do this computation?

[Mario Carneiro]

yes? I'm assuming the number is being represented in binary here

[Mario Carneiro]

that's why log2 in particular is completely trivial even if the number is gigabytes large

[Antoine Chambert-Loir]

Indeed, that's almost a shame to have to make these computations to get the number of binary digits. Couldn't it be more directly accessible ?