Library additions.fib

author Yves Bertot, Inria

From Coq Require Import Extraction ZArith Lia.

From mathcomp Require Import all_ssreflect all_algebra.

Set Implicit Arguments.

Fixpoint fib (n : nat) :=
  if n is S (S p as p') then fib p + fib p' else 1.

Fixpoint fibt (n : nat) (acc1 acc2 : nat) : nat :=
  if n is S p then fibt p (acc2 + acc1) acc1 else acc1.

Fixpoint Zfibt (n : nat) (acc1 acc2 : Z) : Z :=
  if n is S p then Zfibt p (Z.add acc2 acc1) acc1 else acc1.

Lemma fibt_aux (n k: nat) : fibt n (fib k.+1) (fib k) = fib (n + k).+1.

Lemma fibtP (n : nat) : fibt n 1 0 = fib n.

Fixpoint bits p acc : list bool :=
  match p with
  | xHtrue :: acc
  | xI pbits p (true :: acc)
  | xO pbits p (false :: acc)
  end.

Lemma bits_cat p a : bits p a = bits p [::] ++ a.

Lemma bitsP p :
  Pos.to_nat p =
  \sum_(i < size (bits p [::])) nth false (bits p [::])
      (size (bits p [::]) - 1 - i) × 2 ^ i.

Section with_matrices.

Variable R : comUnitRingType.

Import GRing.Theory.
Open Scope ring_scope.

Definition ZtoR (x : Z) : R :=
  if (Z.ltb x 0) then
    -(Z.abs_nat (-x))%:R else (Z.abs_nat x)%:R.

Lemma ZtoRD (x y : Z) : ZtoR (x + y) = ZtoR x + ZtoR y.

Lemma ZtoRM x y : ZtoR (x × y) = ZtoR x × ZtoR y.

Lemma iter_mul (n : nat) (m v : 'M[R]_2) :
  iter n [eta *%R m] v = m ^+ n × v.

Definition fibm : 'M[R]_2 :=
  \matrix_(i, j) if (val i == 1%N) && (val j == 1%N) then 0%R else 1%R.

Lemma fibmP n : fibm ^+ n.+2 =
  \matrix_(i, j)
     if (val i == 0%N) && (val j == 0%N) then (fib n.+2)%:R else
     if ((val i + val j == 1)%N) then (fib n.+1)%:R else
      (fib n)%:R.


Definition m4lval (l : seq Z) (i j : nat) :=
  nth Z0 l (j + 2 × i).

Definition m4lmx (l : seq Z) : 'M[R]_2 :=
  \matrix_(i, j) (ZtoR (m4lval l i j)).

Definition m4lmul (l1 l2 : seq Z) :=
  (Z.add (m4lval l1 0 0 × m4lval l2 0 0) (m4lval l1 0 1 × m4lval l2 1 0)) ::
  (Z.add (m4lval l1 0 0 × m4lval l2 0 1) (m4lval l1 0 1 × m4lval l2 1 1)) ::
  (Z.add (m4lval l1 1 0 × m4lval l2 0 0) (m4lval l1 1 1 × m4lval l2 1 0)) ::
  (Z.add (m4lval l1 1 0 × m4lval l2 0 1) (m4lval l1 1 1 × m4lval l2 1 1)) :: nil.

Open Scope Z_scope.
Definition m4lfib :=
  [:: 1; 1;
      1; 0].

Close Scope Z_scope.


Definition m3lmx (l : seq Z) : 'M[R]_2 :=
  \matrix_(i, j)
   if (i == j) then
     ZtoR (nth Z0 l i)
   else
     ZtoR (nth Z0 l 2).

Definition m3lmul (l1 l2 : seq Z) :=
  (Z.add (nth Z0 l1 0 × nth Z0 l2 0) (nth Z0 l1 2 × nth Z0 l2 2)) ::
  (Z.add (nth Z0 l1 2 × nth Z0 l2 2) (nth Z0 l1 1 × nth Z0 l2 1)) ::
(Z.add (nth Z0 l1 0 × nth Z0 l2 2) (nth Z0 l1 2 × nth Z0 l2 1)) :: nil.

Definition m3lfib := [:: Z.pos xH ; Z0; Z.pos xH].

Fixpoint fastexp (m : list Z) (p : positive) : list Z :=
  match p with
  | xHm
  | xI p
    let r := fastexp m p in
    m3lmul (m3lmul r r) m
  | xO p
    let r := fastexp m p in
      m3lmul r r
  end.

Fixpoint fastexp2 (m : list Z) (p : positive) (acc : list Z) : list Z :=
  match p with
  | xHacc
  | xO pfastexp2 m p (m3lmul acc acc)
  | xI pfastexp2 m p (m3lmul m (m3lmul acc acc))
  end.

Definition fastexp3 {A : Type} (mul : A A A)
  (m : A) :=
  fix f (l : list bool) (acc : A) : A :=
  match l with
  | nilacc
  | true :: lf l (mul (mul acc acc) m)
  | false :: lf l (mul acc acc)
  end.

Definition my_pow {A : Type} (mul : A A A) (m : A)
  (p : positive)
  : A :=
  fastexp3 mul m (behead (bits p nil)) m.

Definition m3lid := [:: Z.pos xH; Z0; Z0; Z.pos xH].

Definition slowexp (m : list Z) p :=
  Pos.iter (m3lmul m) m3lid p.




Definition binary_power_mult (mul : list Z list Z list Z) :
  list Z list Z positive list Z :=
  fix f (x a : list Z) (p : positive) : list Z :=
  match p with
  | xHmul a x
  | xO qf (mul x x) a q
  | xI qf (mul x x) (mul a x) q
end.

Definition fastexp4 (mul : list Z list Z list Z)
  : list Z positive list Z :=
  fix f (x : list Z) (p : positive) :=
  match p with
  | xHx
  | xO qf (mul x x) q
  | xI qbinary_power_mult mul (mul x x) x q
  end.

Lemma m3lmulP l1 l2:
  GRing.comm (m3lmx l1) (m3lmx l2)
  m3lmx (m3lmul l1 l2) = m3lmx l1 × m3lmx l2.

Lemma m3lfibP : m3lmx m3lfib = fibm.

Lemma iter_comm {A : Type} (op : A A A) (e1 e2 : A)
  (assoc: associative op)(cm : op e1 e2 = op e2 e1) n :
  op e1 (iter n (op e1) e2) = op (iter n (op e1) e2) e1.

Lemma iter_combine {A : Type} (op : A A A) (e1 e2 : A)
  (assoc: associative op)(cm : op e1 e2 = op e2 e1) n m :
  op (iter n (op e1) e2) (iter m (op e1) e2) =
  (iter (n + m) (op e1) (op e2 e2)).

Lemma fastexp3P {A : Type} (op : A A A) (id e : A) (h l : nat)
  (v : list bool) (assoc : associative op) (cm : op id e = op e id)
  (idn : a, op id a = a):
  l = (\sum_(i < size v) nth false v (size v - 1 - i) × 2 ^ i)%N
  iter (h × 2 ^ size v + l) (op e) id =
  fastexp3 op e v (iter h (op e) id).

Lemma headbits p l : bits p l = true :: behead (bits p l).

Lemma my_powP m p :
  @my_pow 'M[R]_2 *%R m p = m ^ Pos.to_nat p.

Lemma my_pow_m3lmul m p :
  m3lmx (my_pow m3lmul m p) = my_pow *%R (m3lmx m) p.

Definition m2lmul : list Z list Z list Z :=
 fun l1 l2
 let a := nth Z0 l1 0 in
 let b := nth Z0 l1 1 in
 let c := nth Z0 l2 0 in
 let d := nth Z0 l2 1 in
   [:: Z.add (a × (c + d)) (b × c) ; Z.add (a × c) (b × d)].

Definition m2lmx (l : list Z) : 'M[R]_2 :=
  let a := nth Z0 l 0 in
  let b := nth Z0 l 1 in
  \matrix_(i, j)
      if ((val i == 0%N) && (val j == 0%N)) then
         ZtoR (a + b)
      else if (i == j) then
         ZtoR b
      else ZtoR a.

Definition m2lfib := [:: Zpos xH; Z0].

Lemma m2lfibP : m2lmx m2lfib = fibm.

Lemma m2lmulP m1 m2 :
  m2lmx (m2lmul m1 m2) = m2lmx m1 × m2lmx m2.

Lemma my_pow_m2lmul m p :
  m2lmx (my_pow m2lmul m p) = my_pow *%R (m2lmx m) p.

Lemma fibZ3P p :
  ZtoR (nth Z0 (my_pow m3lmul m3lfib p) 0) = (fib (Pos.to_nat p))%:R.

Lemma fibZ2P p :
  ZtoR (nth Z0 (my_pow m2lmul m2lfib p) 0 +
        nth Z0 (my_pow m2lmul m2lfib p) 1) = (fib (Pos.to_nat p))%:R.

Definition m3lpow (m : list Z) (n : positive) :=
  fastexp3 m3lmul m (behead (bits n nil)) m.

Definition m2lpow (m : list Z) (n : positive) :=
  fastexp3 m2lmul m (behead (bits n nil)) m.

End with_matrices.

Definition bigarg := 30000%positive.

Extraction "bigfib" Z.ltb Z.div_eucl Pos.iter Z.log2 fastexp4
  my_pow m4lmul m3lmul m2lmul m4lfib m3lfib m2lfib bigarg.