Zulip Chat Archive

Stream: general

Topic: Efficient modular fibonacci


view this post on Zulip Neil Strickland (Mar 24 2019 at 17:03):

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?

view this post on Zulip Chris Hughes (Mar 24 2019 at 17:16):

This works. Takes a little while though.

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

view this post on Zulip Neil Strickland (Mar 24 2019 at 17:23):

Thanks, that's much better.

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:03):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:03):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:08):

The #eval gives 0 with no visible delay.

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:10):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:11):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:11):

I think it's a set_option

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:17):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:19):

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

view this post on Zulip Chris Hughes (Mar 24 2019 at 18:24):

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

view this post on Zulip Chris Hughes (Mar 24 2019 at 18:27):

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

view this post on Zulip Mario Carneiro (Mar 24 2019 at 18:47):

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

view this post on Zulip Mario Carneiro (Mar 24 2019 at 18:48):

that's how norm_num itself works, of course

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:52):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:52):

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

view this post on Zulip Mario Carneiro (Mar 24 2019 at 18:55):

isn't that what he's doing?

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:57):

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

view this post on Zulip Chris Hughes (Mar 24 2019 at 18:57):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 18:57):

^^^

view this post on Zulip Mario Carneiro (Mar 24 2019 at 18:59):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:02):

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.

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:03):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:04):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:04):

and Lean never simplifies anything until the end

view this post on Zulip Mario Carneiro (Mar 24 2019 at 19:05):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:05):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:06):

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

view this post on Zulip Chris Hughes (Mar 24 2019 at 19:07):

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.

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:07):

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.

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:09):

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

view this post on Zulip Mario Carneiro (Mar 24 2019 at 19:10):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 19:21):

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

view this post on Zulip Mario Carneiro (Mar 24 2019 at 19:22):

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

view this post on Zulip Neil Strickland (Mar 24 2019 at 20:07):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 20:13):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 20:14):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 20:14):

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

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 20:15):

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)

view this post on Zulip Kevin Buzzard (Mar 24 2019 at 20:17):

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.

view this post on Zulip Kenny Lau (Mar 24 2019 at 21:56):

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)

view this post on Zulip Mario Carneiro (Mar 24 2019 at 22:45):

reduce mod 44?

view this post on Zulip Kenny Lau (Mar 24 2019 at 22:56):

no, Fibonacci

view this post on Zulip Kenny Lau (Mar 24 2019 at 23:20):

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

view this post on Zulip Kenny Lau (Mar 24 2019 at 23:24):

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

view this post on Zulip Kenny Lau (Mar 24 2019 at 23:24):

O(log 2) algorithm ^

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:33):

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

Periods of Fibonacci numbers modulo primes are well understood.

view this post on Zulip Kenny Lau (Mar 25 2019 at 00:36):

is it?

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:37):

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

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:37):

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

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:38):

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

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:41):

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

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:43):

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

view this post on Zulip Kevin Buzzard (Mar 25 2019 at 00:43):

so their 44th powers are both 1.

view this post on Zulip Kenny Lau (Mar 25 2019 at 01:40):

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

view this post on Zulip Kenny Lau (Mar 25 2019 at 01:40):

but interesting observation


Last updated: May 18 2021 at 16:25 UTC