Binet Formula

Stream: new members

[Anatole Dedecker]

I just spent a few hours defining the golden ratio and proving Binet formula for the Fibonacci sequence. Is this already in mathlib ? If not, should I PR it ?

[Anatole Dedecker]

Btw, ring_exp is awesome

[Bryan Gin-ge Chen]

I don't think we have this yet and I think it is worth PRing.

[Alex J. Best]

This is a codewars kata https://www.codewars.com/kata/5d2b89d4b90c0a001f4a6456/lean but I don't think that should stop us from having Binet in mathlib (however some peoples solutions on there which you can check out if you submit your solution might have cool methods worth considering), which proof did you do?

[Anatole Dedecker]

Oh I didn't know about this. My proof is based on this predicate :

def satisfy_fib_rec (u : ℕ → α) :=
  ∀ n, u (n+2) = u n + u (n+1)

which is stable under weighted sum (is this the right way to say that ?)

lemma satisfy_fib_rec_linear {u v : ℕ → α} {a b : α}
  (hu : satisfy_fib_rec u) (hv : satisfy_fib_rec v) : satisfy_fib_rec (λ n, a * u n + b * v n) :=
sorry

I also prove by induction that a sequence $u$ such that $u_0 = 0$ and $u_1 = 1$ which satisfies the above predicate is equal to the Fibonacci sequence :

lemma eq_fib_of_satisfy_fib_rec (u : ℕ → ℝ) (h₀ : u 0 = 0) (h₁ : u 1 = 1)
  (hₙ : satisfy_fib_rec u) :
  ∀ n, u n = n.fib :=
sorry

Then I just have to show that λ n, (φ^n) and λ n, ((-φ⁻¹)^n) satisfy this predicate, and what is left is just computation

[Anatole Dedecker]

I had a look to all answers. Mine was clearly not the shortest (but I wasn't really trying to make it short). Although maybe some of what I've proved about satisfy_fib_rec could be useful for studying e.g Lucas Number, I think I won't use it in a PR, there are much better proofs. @Alexey Solovyev @Junyan Xu do you mind if I use parts of your solutions to this kata for making the mathlib proof ?

[Alexey Solovyev]

@Anatole Dedecker You may use any part of my solution for the mathlib proof. In general, all Codewars solutions are licensed under the FreeBSD 2-Clause License (https://www.codewars.com/about/terms-of-service).

[Alex J. Best]

Yeah that would be my suggestion too it would be cool to do something a bit more general. Like Lucas numbers too or any length 2 linear recurrence, but I didnt want to hold up your PR!

[Anatole Dedecker]

Well, there's no hurry and that's not really complicated maths, I can definitely work on it if you think it's worth to have in mathlib

[Alex J. Best]

Having the version for any linear recurrence would be a great end goal, I think its definitely worth to have that in mathlib eventually. I.e. the solution of any https://en.wikipedia.org/wiki/Constant-recursive_sequence like in theorem 1 of this random pdf https://www.math.kth.se/math/GRU/2012.2013/SF1610/CINTE/mastertheorem.pdf . The case where the roots have multiplicity one is an important special case. These sorts of recurrences really do come up in modern number theory anyway :smile:

[Junyan Xu]

Glad you like my proof, and you are welcome to use it! Not sure about mathlib contribution policies but you may mention my github handle @alreadydone where appropriate.
Anatole Dedecker said:

I had a look to all answers. Mine was clearly not the shortest (but I wasn't really trying to make it short). Although maybe some of what I've proved about satisfy_fib_rec could be useful for studying e.g Lucas Number, I think I won't use it in a PR, there are much better proofs. Alexey Solovyev Junyan Xu do you mind if I use parts of your solutions to this kata for making the mathlib proof ?

[Anatole Dedecker]

Alex J. Best said:

Having the version for any linear recurrence would be a great end goal, I think its definitely worth to have that in mathlib eventually. I.e. the solution of any https://en.wikipedia.org/wiki/Constant-recursive_sequence like in theorem 1 of this random pdf https://www.math.kth.se/math/GRU/2012.2013/SF1610/CINTE/mastertheorem.pdf . The case where the roots have multiplicity one is an important special case. These sorts of recurrences really do come up in modern number theory anyway :)

Ok so I started working on this, and after a long time of hesitations I've settled down on the following definitions. Before I continue, could you tell me if they look usable ? (Same question about the statement of target theorem`)

import data.polynomial.ring_division
import linear_algebra.dimension
import algebra.polynomial.big_operators

open finset

open_locale big_operators

noncomputable theory

structure linear_recurrence (α : Type*) [field α] := (order : ℕ) (coeffs : fin order → α)

variables {α : Type*} [field α]

namespace linear_recurrence

variables (E : linear_recurrence α)

def is_solution (u : ℕ → α) :=
  ∀ n, u (n + E.order) = (∑ i, E.coeffs i * u (n + i))

def mk_sol (init : fin E.order → α) : ℕ → α
| n := if h : n < E.order then init ⟨n, h⟩ else
  ∑ k : fin E.order, have n - E.order + k.1 < n, by have := k.2; omega,
    E.coeffs k * mk_sol (n - E.order + k)

def is_sol_mk_sol (init : fin E.order → α) :
  E.is_solution (E.mk_sol init) :=
  λ n, by rw mk_sol; simp

def mk_sol_eq_init (init : fin E.order → α) :
  ∀ n : fin E.order, E.mk_sol init n = init n :=
  λ n, by rw mk_sol; simp [fin.coe_eq_val, n.2]

lemma eq_mk_of_is_sol_of_eq_init {u : ℕ → α} {init : fin E.order → α}
  (h : E.is_solution u) (heq : ∀ n : fin E.order, u n = init n) :
  ∀ n, u n = E.mk_sol init n
| n := if h' : n < E.order
  then by rw mk_sol; simp only [h', dif_pos]; exact_mod_cast heq ⟨n, h'⟩
  else begin
    rw mk_sol,
    rw ← nat.sub_add_cancel (le_of_not_lt h'),
    rw h (n-E.order),
    simp [h'],
    congr',
    ext k,
    exact have wf : n - E.order + k.1 < n, by have := k.2; omega,
      by rw eq_mk_of_is_sol_of_eq_init
  end

def sol_space : subspace α (ℕ → α) :=
{ carrier := {u | E.is_solution u},
  zero_mem' := λ n, by simp,
  add_mem' := λ u v hu hv n, by simp [mul_add, sum_add_distrib, hu n, hv n],
  smul_mem' := λ a u hu n, by simp [hu n, mul_sum]; congr'; ext; ac_refl }

lemma is_sol_iff_mem_sol_space (u : ℕ → α) :
  E.is_solution u ↔ u ∈ E.sol_space := by refl

def sol_space_linear_equiv_init_space :
  E.sol_space ≃ₗ[α] (fin E.order → α) :=
{ to_fun := λ u x, (u : ℕ → α) x,
  map_add' := λ u v, by ext; simp; norm_cast,
  map_smul' := λ a u, by ext; simp,
  inv_fun := λ u, ⟨E.mk_sol u, E.is_sol_mk_sol u⟩,
  left_inv := λ u, by ext n; symmetry; apply E.eq_mk_of_is_sol_of_eq_init u.2; intros k; refl,
  right_inv := λ u, function.funext_iff.mpr (λ n, E.mk_sol_eq_init u n) }

lemma sol_space_dim :
  vector_space.dim α E.sol_space = E.order :=
@dim_fin_fun α _ E.order ▸ E.sol_space_linear_equiv_init_space.dim_eq

def char_poly : polynomial α :=
  polynomial.monomial E.order 1 - (∑ i : fin E.order, polynomial.monomial i (E.coeffs i))

#check polynomial.eval₂_finset_sum

lemma geom_sol_iff_root_char_poly (q : α) :
  E.is_solution (λ n, q^n) ↔ E.char_poly.is_root q :=
begin
  rw [char_poly, polynomial.is_root.def, polynomial.eval],
  simp only [polynomial.eval₂_finset_sum, one_mul,
              ring_hom.id_apply, polynomial.eval₂_monomial, polynomial.eval₂_sub],
  split,
  { intro h,
    specialize h 0,
    simp only [zero_add] at h,
    rwa ← sub_eq_zero_iff_eq at h },
  { intros h n,
    rw sub_eq_zero_iff_eq at h,
    simp only [pow_add],
    rw [h, mul_sum],
    congr,
    ext,
    ring }
end

theorem target (h : ∑ x in E.char_poly.roots, E.char_poly.root_multiplicity x = E.order) (u : ℕ → α) :
  E.is_solution u ↔
    ∃ (a : α → polynomial α)
      (h : ∀ x, (a x).degree < E.char_poly.root_multiplicity x),
    u = λ n, ∑ x in E.char_poly.roots, (a x).eval n * x^n := sorry

end linear_recurrence

[Kevin Buzzard]

field \a should be something like comm_semiring \a. I don't see any subtraction or division in alpha here.

[Kevin Buzzard]

But this looks like a great set-up in general. Did you try using it to make sure it's usable for e.g. Binet?

[Anatole Dedecker]

Well yes but then I can't define the vector space of solutions

[Kevin Buzzard]

sure you can, it's the semimodule of solutions

[Kevin Buzzard]

the axioms for a vector space don't mention subtraction or division either

[Anatole Dedecker]

Kevin Buzzard said:

sure you can, it's the semimodule of solutions

What about dimensions ? I don't see them defined for anything else than a vector_space

[Kevin Buzzard]

only free semimodules have dimensions

[Anatole Dedecker]

Kevin Buzzard said:

But this looks like a great set-up in general. Did you try using it to make sure it's usable for e.g. Binet?

Oh I should try this indeed

[Kevin Buzzard]

but fin E.order → α is free

[Kevin Buzzard]

So now an interesting question is the following:

[Kevin Buzzard]

Say that the auxiliary polynomial (a monic polynomial with coefficients in alpha) is a product of linear factors (X-a1)(X-a2)...(X-an). When is (1,a1,a1^2,...,a1^{n-1}),(1,a2,a2^2,...,a2^{n-1}),...,(1,an,an^2,...an^{n-1}) a basis? And here you need a field.

[Kevin Buzzard]

Because even for the integers if you consider X^2-1, i.e. a_{n+2}=a_n, then the roots are +-1, so the power sequences are (1,1,1,...) and (1,-1,1,-1,...), and you cannot find a Z-linear combination of these which makes (1,0,1,0,...) because you can't divide by 2.

[Kevin Buzzard]

So I suggest that all the beginning stuff is done for semirings (you probably don't even need commutative if you say which order you are summing things in, but perhaps this is going too far), but if you want to try to prove the "first main theorem of recurrence relations", which is that if the aux poly has distinct roots in the ground ring then (1,a1,a1^2,...),(1,a2,a2^2,...) etc form a basis, then you should do this for a field.

[Anatole Dedecker]

Yep that looks good, I'll do this !

[Kevin Buzzard]

And then there are more fundamental theorems, for example with repeated roots there is some other solution, something like (0,1a,2a^2,3*a^3), and these give a basis, and probably you can also say something when the roots are not in the ground field too -- but actually I think that these more general theorems should be proved another way.

[Kevin Buzzard]

Let's stick to fields for this more subtle stuff. Then the vector space of solutions is n-dimensional because of the linear equiv you defined. There is also an operator sending (s0,s1,s2,s3,...) to (s1,s2,s3,...), and the point of the (1,a,a^2,a^3,...) solutions is that they are eigenvectors for this operator. The aux poly is the char poly of this operator. So actually you should prove the (1,a,a^2,..) stuff is a basis by deducing it from the general statement that an endomorphism of a fdvs (evdf) whose char poly has distinct linear factors has a basis of eigenvectors.

[Kevin Buzzard]

and the more general stuff comes from theorems like Jordan decomposition.

[Anatole Dedecker]

Hmmm I have read about this approach, but I thought we didn't have eigen stuff in mathlib :thinking:

[Kevin Buzzard]

Yes but it's not far away. I guess what I'm saying is that the natural thing to do is to wait until these theorems are in mathlib

[Anatole Dedecker]

Oh okay. So should I PR the infrastructure and then wait a bit before proving my target lemma ?

[Anatole Dedecker]

Btw is there a name for the "shift n consecutive terms" operator ?

[Anatole Dedecker]

Kevin Buzzard said:

only free semimodules have dimensions

I a free semimodule's dimension defined in mathlib ? I'm struggling to find it

[Kevin Buzzard]

I'm not sure that it will be there. You are proving that it's isomorphic to fin n -> alpha which is all you need.