Efficient modular fibonacci

Stream: general

[Neil Strickland]

I want to prove that the 44th Fibonacci number is divisible by 89. I wrote:

import tactic.norm_num

def FP_step : ℕ × ℕ → ℕ × ℕ :=
 λ ⟨a,b⟩, ⟨b,(a + b) % 89⟩

def FP : ℕ → ℕ × ℕ
| 0 := ⟨0,1⟩
| (n + 1) := FP_step (FP n)

def F (n : ℕ) : ℕ := (FP n).fst

#eval F 44

lemma L : F 44 = 0 := sorry

The #eval gives 0 with no visible delay. If I try to prove L by rfl then I get a deterministic timeout, which I find surprising; I am only asking for 44 additions and 44 reductions mod 89 of numbers that are < 2 * 89, which does not seem too bad even in Peano-world. It also seems that norm_num is not useful here (unless there is some non-obvious way to invoke it). Am I misunderstanding something about the count of operations? Is there something different that I should be doing?

[Chris Hughes]

This works. Takes a little while though.

lemma L : F 44 = 0 := by unfold F FP FP_step; norm_num

[Neil Strickland]

Thanks, that's much better.

[Kevin Buzzard]

Related: https://github.com/kbuzzard/lean-squares-in-fibonacci

[Kevin Buzzard]

Unfinished proof that the largest square in the Fibonacci sequence is 144.

[Kevin Buzzard]

The #eval gives 0 with no visible delay.

Of course it's #reduce that is the problem...

[Kevin Buzzard]

import tactic.norm_num

def FP_step : ℕ × ℕ → ℕ × ℕ :=
 λ ⟨a,b⟩, ⟨b,(a + b) % 89⟩

def FP : ℕ → ℕ × ℕ
| 0 := ⟨0,1⟩
| (n + 1) := FP_step (FP n)

def F (n : ℕ) : ℕ := (FP n).fst

lemma L : F 44 = 0 := begin
  unfold F FP FP_step,
  -- worth putting your cursor here and inspecting the goal
  norm_num
  end

[Kevin Buzzard]

How do you switch those ...'s off again?

[Kevin Buzzard]

I think it's a set_option

[Kevin Buzzard]

Still some left:

import tactic.norm_num

def FP_step : ℕ × ℕ → ℕ × ℕ :=
 λ ⟨a,b⟩, ⟨b,(a + b) % 89⟩

def FP : ℕ → ℕ × ℕ
| 0 := ⟨0,1⟩
| (n + 1) := FP_step (FP n)

def F (n : ℕ) : ℕ := (FP n).fst

set_option pp.max_depth 10000
set_option pp.max_steps 10000
set_option pp.indent 0


lemma L : F 44 = 0 := begin
  unfold F FP FP_step,
  -- worth putting your cursor here and inspecting the goal
  norm_num
  end

[Kevin Buzzard]

That looks to me like a few more than 44 additions...

[Chris Hughes]

I guess it's not evaluating it in a sensible order. After unfolding, you probably end up with F 10 a lot of times.

[Chris Hughes]

I wonder what the fastest proof is, or how to solve this with bigger numbers.

[Mario Carneiro]

You can solve it in linear time if you craft appropriate lemmas

[Mario Carneiro]

that's how norm_num itself works, of course

[Kevin Buzzard]

The fastest proof is to do it by pure thought ;-)

[Kevin Buzzard]

But if you want a computational proof, just work out F n mod 89 for the first 45 naturals

[Mario Carneiro]

isn't that what he's doing?

[Kevin Buzzard]

It didn't look like it from the prettyprinter output.

[Chris Hughes]

F 44 = F 42 + (F 42 + F 41), so you end up computing F 42 twice, and it get worse the smaller you get.

[Kevin Buzzard]

^^^

[Mario Carneiro]

but that's not how it was defined... I don't think this has exponential blowup

[Kevin Buzzard]

My 16 year old son says that he once saw some way of defining Fibonacci numbers in a functional language which didn't suffer from this blow-up.

[Kevin Buzzard]

You're right though, this isn't like what Chris says, this looks like a more sensible implementation.

[Kevin Buzzard]

The point somehow seems to be that each time you iterate, the term gets about 1.5 times as complicated

[Kevin Buzzard]

and Lean never simplifies anything until the end

[Mario Carneiro]

I have a solution but I'm waiting on lean to compile to test it -_-

[Kevin Buzzard]

F 8 looks like this:

⊢ (((((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89 +
   ((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89) %
  89 +
  (((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89 +
   (((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89 +
    ((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89) %
   89) %
  89) %
89,
((((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89 +
   (((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89 +
    ((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89) %
   89) %
  89 +
  ((((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89 +
    ((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89) %
   89 +
   (((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89 +
    (((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89 +
     ((1 + (0 + 1) % 89) % 89 + ((0 + 1) % 89 + (1 + (0 + 1) % 89) % 89) % 89) % 89) %
    89) %
   89) %
  89) %
89).fst

[Kevin Buzzard]

FP 8 looks like a mess because FP 7 is a mess and Lean doesn't tidy it up before making FP 8

[Chris Hughes]

So it is only ever adding 0 and 1. 1.5 times bigger each time is about right for the really bad implementation given that the golden ratio is ~1.6.

[Kevin Buzzard]

Right, each term is about 1.6 times bigger than the one before in terms of e.g. number of additions mentioned in the term.

[Kevin Buzzard]

Similarly for the number of % 89's. If a mentions X of them and b mentions Y, then (a+b)%89 mentions X+Y+1, and this is growing like Fibonacci

[Mario Carneiro]

this works:

lemma FP_succ {n a b c : ℕ}
  (h : FP n = (a, b)) (h2 : (a + b) % 89 = c) : FP (nat.succ n) = (b, c) :=
by rw [← h2, FP, h]; refl

lemma L : F 44 = 0 := begin
  have : FP 0 = (0, 1) := rfl,
  iterate 44 { replace := FP_succ this (by norm_num; refl) },
  exact congr_arg prod.fst this
end

[Kevin Buzzard]

  iterate 44 { replace := FP_succ this (by norm_num; refl) }, -- put your cursor here for some goal fun

[Mario Carneiro]

it's possible to reduce the succ's but it's going to generate 44 add-one's anyway so I didn't bother

[Neil Strickland]

Thanks, that's interesting. I see roughly how it works but will need to digest a bit further.

[Kevin Buzzard]

I guess replace := ... is shorthand for replace this := ...which means "replace the hypothesis this with a new hypothesis also called this.

[Kevin Buzzard]

lemma L : F 44 = 0 := begin
  have this : FP 0 = (0, 1) := rfl,
  iterate 5 { replace this := FP_succ this (by norm_num; refl) },
  -- this : FP (nat.succ (nat.succ (nat.succ (nat.succ 1)))) = (5, 8)
  sorry
end

[Kevin Buzzard]

A hypothesis without a name defaults to this and Mario is abusing this a bit :-)

[Kevin Buzzard]

But replace is essential to the method because we don't have the machinery to make 44 names without going into meta land (as far as I know)

[Kevin Buzzard]

This works:

lemma L : F 44 = 0 := begin
  have H : FP 0 = (0, 1) := rfl,
  iterate 44 { have H := FP_succ H (by norm_num; refl) },
  exact congr_arg prod.fst H,
end

but you end up with 44 hypotheses all called H. The reason it works is that Lean uses the most recently generated one at each stage.

[Kenny Lau]

I feel compulsory to note that although this modified algorithm is O(n), which is better than the naive algorithm of O(2^n), there is an algorithm of O(log n)

[Mario Carneiro]

reduce mod 44?

[Kenny Lau]

no, Fibonacci

[Kenny Lau]

it only dawned on me that F 11 = 89 which would imply that F 44 is divisible by 89

[Kenny Lau]

import tactic.norm_num

constant F : ℕ → ℕ
axiom F_0 : F 0 = 0
axiom F_1 : F 1 = 1
axiom F_2 : F 2 = 1
axiom F_succ_succ {n : ℕ} : F (n.succ.succ) = (F n + F n.succ) % 89
axiom F_bit0 {n} : F (bit0 n) = ((F n * F n.succ) % 89 * 2 % 89 + 89 - (F n * F n) % 89) % 89
axiom F_bit1 {n} : F (bit1 n) = ((F n * F n) % 89 + (F n.succ * F n.succ) % 89) % 89
section
variables {n : ℕ} (m : ℕ) (H : n.succ = m)
theorem F_bit0' : F (bit0 n) = ((F n * F m) % 89 * 2 % 89 + 89 - (F n * F n) % 89) % 89 := H ▸ F_bit0
theorem F_bit1' : F (bit1 n) = ((F n * F n) % 89 + (F m * F m) % 89) % 89 := H ▸ F_bit1
end
theorem F_succ_succ' (n m k : ℕ) (H1 : n.succ = m) (H2 : m.succ = k) : F k = (F n + F m) % 89 := H2 ▸ H1 ▸ F_succ_succ

-- (44, 45) -> (22, 23) -> (10, 11) -> (4, 5) -> (2, 3) -> (0, 1)
theorem F0 : F 0 = 0 := F_0
theorem F1 : F 1 = 1 := F_1
theorem F2 : F 2 = 1 := F_2
theorem F3 : F 3 = 2 := by rw [F_bit1' 2 rfl, F1, F2]; refl
theorem F4 : F 4 = 3 := by rw [F_bit0' 3 rfl, F2, F3]; norm_num
theorem F5 : F 5 = 5 := by rw [F_bit1' 3 rfl, F2, F3]; refl
theorem F6 : F 6 = 8 := by rw [F_succ_succ' 4 5 6 rfl rfl, F4, F5]; refl
theorem F10 : F 10 = 55 := by rw [F_bit0' 6 rfl, F5, F6]; norm_num
theorem F11 : F 11 = 0 := by rw [F_bit1' 6 rfl, F5, F6]; norm_num
theorem F12 : F 12 = 55 := by rw [F_succ_succ' 10 11 12 rfl, F10, F11]; norm_num
theorem F22 : F 22 = 0 := by rw [F_bit0' 12 rfl, F11, F12]; norm_num
theorem F23 : F 23 = 88 := by rw [F_bit1' 12 rfl, F11, F12]; norm_num
theorem F44 : F 44 = 0 := by rw [F_bit0' 23 rfl, F22, F23]; norm_num
theorem F45 : F 45 = 1 := by rw [F_bit1' 23 rfl, F22, F23]; norm_num

[Kenny Lau]

O(log 2) algorithm ^

[Kevin Buzzard]

The fastest proof is to do it by pure thought ;-)

Periods of Fibonacci numbers modulo primes are well understood.

[Kenny Lau]

is it?

[Kevin Buzzard]

A slight shortcut in what you did Kenny -- you can show that $F_{2n}$ is a multiple of $F_n$ saving you from having to work out $F_{12}$ and $F_{23}$

[Kevin Buzzard]

It's just about the order of $\frac{1+5^{1/2}}{2}$ modulo $p$ right?

[Kevin Buzzard]

$(1+5^{1/2})^{11}=(1-5^{1/2})^{11}$ and we're done, right?

[Kevin Buzzard]

The pure thought proof says that the order divides $p-1$ (note that $5$ is a square mod $89$) so $F_{p-1}$ is automatically 0 mod $p$, and $F_{(p-1)/2}$ will be $0$ mod $p$ if $(1+5^{1/2})/2$ is a quadratic residue mod $p$, which it is when $p=89$.

[Kevin Buzzard]

I'm also using that $p=1$ mod 4 here, because then it follows that $(1-5^{1/2})/2$ is also a QR

[Kevin Buzzard]

so their 44th powers are both 1.

[Kenny Lau]

well I want F45 also to show that the period is 44

[Kenny Lau]

but interesting observation