Require Import List.
Require Export Complex.
Require Export Quantum.
Require Import FTA.

Ltac Hhelper :=
   unfold Mmult;
   unfold big_sum;
   unfold I;
   simpl;
   C_field_simplify;
   try lca;
   C_field.

Lemma Y_eq_iXZ : σy = Ci .* σx × σz.
Lemma H_eq_Hadjoint : hadamard† = hadamard.

Hint Rewrite Y_eq_iXZ H_eq_Hadjoint : Q_db.

Lemma ItimesIid : I 2 × I 2 = I 2.
Lemma XtimesXid : σx × σx = I 2.
Lemma YtimesYid : σy × σy = I 2.
Lemma ZtimesZid : σz × σz = I 2.
Lemma HtimesHid : hadamard × hadamard = I 2.

Hint Rewrite ItimesIid XtimesXid YtimesYid ZtimesZid HtimesHid : Q_db.

Lemma ZH_eq_HX : σz × hadamard = hadamard × σx.
Lemma XH_eq_HZ : σx × hadamard = hadamard × σz.
Lemma SX_eq_YS : Sgate × σx = σy × Sgate.
Lemma SZ_eq_ZS : Sgate × σz = σz × Sgate.
Lemma cnotX1 : cnot × (σx ⊗ I 2) = (σx ⊗ σx) × cnot.
Lemma cnotX2 : cnot × (I 2 ⊗ σx) = (I 2 ⊗ σx) × cnot.
Lemma cnotZ1 : cnot × (σz ⊗ I 2) = (σz ⊗ I 2) × cnot.
Lemma cnotZ2 : cnot × (I 2 ⊗ σz) = (σz ⊗ σz) × cnot.

Hint Rewrite ZH_eq_HX XH_eq_HZ SX_eq_YS SZ_eq_ZS cnotX1 cnotX2 cnotZ1 cnotZ2 : Q_db.

Local Open Scope nat_scope.

Definition orthogonal {n m} (S : Matrix n m) : Prop :=
  ∀ i j, i ≠ j → inner_product (get_vec i S) (get_vec j S) = C0.

Definition orthonormal {n m} (S : Matrix n m) : Prop :=
  orthogonal S ∧ (∀ (i : nat), i < m → norm (get_vec i S) = 1%R).

Definition WF_Orthonormal {n m} (S : Matrix n m) : Prop :=
  WF_Matrix S ∧ orthonormal S.

Lemma inner_product_is_mult : ∀ {n} (i j : nat) (S : Square n),
  inner_product (get_vec i S) (get_vec j S) = (S† × S) i j.

Lemma inner_product_comm_conj : ∀ {n} (v u : Vector n),
  inner_product v u = Cconj (inner_product u v).

Definition MatrixP (m n : nat) := nat → nat → Polynomial.

Notation SquareP n := (MatrixP n n).

Definition eval_matP {n m} (A : MatrixP n m) (c : C) : Matrix n m :=
  fun x y ⇒ (A x y)[[c]].

Definition reduceP {n} (A : SquareP (S n)) (row col : nat) : SquareP n :=
  fun x y ⇒ (if x <? row
              then (if y <? col
                    then A x y
                    else A x (1+y))
              else (if y <? col
                    then A (1+x) y
                    else A (1+x) (1+y))).

Lemma reduceP_eval_mat : ∀ {n} (A : SquareP (S n)) (c : C) (x y : nat),
  reduce (eval_matP A c) x y = eval_matP (reduceP A x y) c.

Fixpoint DeterminantP (n : nat) (A : SquareP n) : Polynomial :=
  match n with
  | 0 ⇒ [C1]
  | S 0 ⇒ A 0 0
  | S n' ⇒ (big_sum (fun i ⇒ [(parity i)] *, (A i 0) *, (DeterminantP n' (reduceP A i 0)))%C n)
  end.

Arguments DeterminantP {n}.

Lemma DetP_simplify : ∀ {n} (A : SquareP (S (S n))),
  DeterminantP A =
  (big_sum (fun i ⇒ [(parity i)] *, (A i 0) *, (DeterminantP (reduceP A i 0)))%C (S (S n))).

Lemma Peval_Det : ∀ {n} (A : SquareP n) (c : C),
  Determinant (eval_matP A c) = (DeterminantP A)[[c]].

Definition prep_mat {n} (A : Square n) : SquareP n :=
  (fun x y ⇒ if (x =? y) && (x <? n) then [A x y; -C1] else [A x y]).


Definition deg_elem_leq_1 {n} (A : SquareP n) : Prop :=
  ∀ i j, degree (A i j) ≤ 1.

Lemma del1_reduce : ∀ {n} (A : SquareP (S n)) (i j : nat),
  deg_elem_leq_1 A → deg_elem_leq_1 (reduceP A i j).

Lemma bound_deg_matP : ∀ {n} (A : SquareP n),
  deg_elem_leq_1 A → degree (DeterminantP A) ≤ n.

Lemma del1_prep_mat : ∀ {n} (A : Square n),
  deg_elem_leq_1 (prep_mat A).

Lemma reduce_prep_mat : ∀ {n} (A : Square (S n)),
  reduceP (prep_mat A) 0 0 = prep_mat (reduce A 0 0).

Lemma detP_deg : ∀ {n} (A : Square n),
  degree (DeterminantP (prep_mat A)) = n.

Lemma connect : ∀ (n : nat) (A : Square (S n)),
  ∃ (p : Polynomial), (Polynomial.degree p) > 0 ∧
    (∀ c : C, Determinant (A .+ (-c .* I (S n))) = p[[c]]).

Lemma connect2 : ∀ (n : nat) (A : Square (S n)),
  ∃ (c : C), det_eq_c C0 (A .+ (-c .* I (S n))).

Lemma exists_eigenvector : ∀ (n : nat) (A : Square (S n)),
  WF_Matrix A →
  ∃ (c : C) (v : Vector (S n)), WF_Matrix v ∧ v ≠ Zero ∧ A × v = c.* v.

Definition form_basis {n} (v : Vector n) (non_zero : nat) : Matrix n n :=
  fun x y ⇒ if (y =? non_zero)
             then (v x 0)
             else (@e_i n y x 0).

Lemma WF_form_basis : ∀ {n} (v : Vector n) (x : nat),
  WF_Matrix v → x < n → WF_Matrix (form_basis v x).

Lemma get_v_in_basis : ∀ {n} (v : Vector n) (x : nat),
  WF_Matrix v → get_vec x (form_basis v x) = v.

Lemma get_ei_in_basis : ∀ {n} (v : Vector n) (x y : nat),
  y < n → y ≠ x → get_vec y (form_basis v x) = e_i y.

Lemma form_basis_ver : ∀ {n} (v : Vector n) (x : nat),
  v ≠ Zero → WF_Matrix v → v x 0 ≠ C0 → x < n →
  linearly_independent (form_basis v x) ∧ get_vec x (form_basis v x) = v.

Lemma lin_indep_out_of_v : ∀ {n} (v : Vector n),
  WF_Matrix v → v ≠ Zero →
  ∃ S : Square n, WF_Matrix S ∧ linearly_independent S ∧ get_vec 0 S = v.

Definition proj {n} (u v : Vector n) : Vector n :=
  ((inner_product u v) / (inner_product u u)) .* u.

Definition proj_coef {n} (u v : Vector n) : C :=
  ((inner_product u v) / (inner_product u u)).

Lemma proj_inner_product : ∀ {n} (u v : Vector n),
  (norm u) ≠ 0%R → inner_product u (proj u v) = inner_product u v.

Definition gram_schmidt_on_v (n m : nat) (v : Vector n) (S : Matrix n m) :=
  v .+ (big_sum (fun i ⇒ (-C1) .* (proj (get_vec i S) v)) m).

Definition delta_T {n m} (T : Matrix n (S m)) (i : nat) : C :=
  match i =? m with
  | true ⇒ C1
  | _ ⇒ - (proj_coef (get_vec i T) (get_vec m T))
  end.

Definition gram_schmidt_on_T (n m : nat) (T : Matrix n (S m)) : Vector n :=
  big_sum (fun i ⇒ (delta_T T) i .* (get_vec i T)) (S m).

Lemma WF_gs_on_T : ∀ {n m} (T : Matrix n (S m)),
  WF_Matrix T → WF_Matrix (gram_schmidt_on_T n m T).

Lemma gram_schmidt_compare : ∀ {n m} (T : Matrix n (S m)),
  inner_product (get_vec m T) (get_vec m T) ≠ C0 →
  gram_schmidt_on_T n m T = gram_schmidt_on_v n m (get_vec m T) (reduce_col T m).

Lemma gram_schmidt_orthogonal : ∀ {n m} (v : Vector n) (S : Matrix n m) (i : nat),
  orthonormal S → i < m → inner_product (get_vec i S) (gram_schmidt_on_v n m v S) = C0.

Definition f_to_vec (n : nat) (f : nat → C) : Vector n :=
  fun i j ⇒ if (j =? 0) && (i <? n) then f i else C0.

Lemma WF_f_to_vec : ∀ (n : nat) (f : nat → C), WF_Matrix (f_to_vec n f).

Lemma Msum_to_Mmult : ∀ {n m} (T : Matrix n (S m)) (f : nat → C),
  big_sum (fun i ⇒ f i .* get_vec i T) (S m) = T × (f_to_vec (S m) f).

Lemma gram_schmidt_non_zero : ∀ {n m} (T : Matrix n (S m)),
  linearly_independent T → gram_schmidt_on_T n m T ≠ Zero.

Lemma inner_product_zero_normalize : ∀ {n} (u v : Vector n),
  inner_product u v = C0 → inner_product u (normalize v) = C0.

Lemma get_vec_reduce_append_miss : ∀ {n m} (T : Matrix n (S m)) (v : Vector n) (i : nat),
  i < m → get_vec i (col_append (reduce_col T m) v) = get_vec i T.

Lemma get_vec_reduce_append_hit : ∀ {n m} (T : Matrix n (S m)) (v : Vector n),
  WF_Matrix v → get_vec m (col_append (reduce_col T m) v) = v.

Lemma get_vec_reduce_append_over : ∀ {n m} (T : Matrix n (S m)) (v : Vector n) (i : nat),
  WF_Matrix T → i > m →
  get_vec i (col_append (reduce_col T m) v) = Zero.

Lemma extend_onb_ind_step_part1 : ∀ {n m} (T : Matrix n (S m)),
  WF_Matrix T → linearly_independent T → orthonormal (reduce_col T m) →
  orthonormal (col_append (reduce_col T m) (normalize (gram_schmidt_on_T n m T))).

Definition delta_T' {n m} (T : Matrix n m) (v : Vector n) (size : nat) : nat → C :=
  fun i ⇒ if (i <? size)
           then - (proj_coef (get_vec i T) v)
           else C0.

Lemma gs_on_T_cols_add : ∀ {n m1 m2} (T1 : Matrix n m1) (T2 : Matrix n m2) (v : Vector n),
  WF_Matrix v →
  smash (col_append T1 (gram_schmidt_on_T n m1 (col_append T1 v))) T2 =
  @col_add_many n ((S m1) + m2) m1 (f_to_vec (m1 + m2) (delta_T' T1 v m1))
                             (smash (col_append T1 v) T2).

Lemma smash_scale : ∀ {n m1 m2} (T1 : Matrix n m1) (T2 : Matrix n m2) (v : Vector n),
  smash (col_append T1 (normalize v)) T2 =
  col_scale (smash (col_append T1 v) T2) m1 (/ norm v).

Lemma extend_onb_ind_step_part2 : ∀ {n m1 m2} (T1 : Matrix n m1) (T2 : Matrix n m2)
                                                   (v : Vector n),
  WF_Matrix T1 → WF_Matrix T2 → WF_Matrix v → v ≠ Zero →
  linearly_independent (smash (col_append T1 v) T2) →
  linearly_independent (smash (col_append T1
                                    (normalize (gram_schmidt_on_T n m1 (col_append T1 v)))) T2).

Lemma extend_onb_ind_step : ∀ {n m1 m2} (T1 : Matrix n m1) (T2 : Matrix n m2) (v : Vector n),
  WF_Matrix T1 → WF_Matrix T2 → WF_Matrix v →
  linearly_independent (smash (col_append T1 v) T2) → orthonormal T1 →
  ∃ v1, WF_Matrix v1 ∧ orthonormal (col_append T1 v1) ∧
             linearly_independent (smash (col_append T1 v1) T2).

Lemma extend_onb : ∀ (n m2 m1 : nat) (T1 : Matrix n (S m1)) (T2 : Matrix n m2),
  WF_Matrix T1 → WF_Matrix T2 →
  linearly_independent (smash T1 T2) → orthonormal T1 →
  ∃ T2' : Matrix n m2, WF_Matrix T2' ∧ orthonormal (smash T1 T2').

Lemma get_vec_vec : ∀ {n} (v : Vector n),
  WF_Matrix v → get_vec 0 v = v.

Lemma orthonormal_normalize_v : ∀ {n} (v : Vector n),
  v ≠ Zero → WF_Matrix v → orthonormal (normalize v).

Theorem onb_out_of_v : ∀ {n} (v : Vector n),
  WF_Matrix v → v ≠ Zero →
  ∃ T : Square n, WF_Orthonormal T ∧ get_vec 0 T = normalize v.

Lemma unit_is_orthonormal : ∀ {n} (U : Square n),
  WF_Unitary U ↔ WF_Orthonormal U.

Lemma unit_out_of_v : ∀ {n} (v : Vector n) (x : nat),
  WF_Matrix v → v ≠ Zero →
  ∃ S : Square n, WF_Unitary S ∧ get_vec 0 S = normalize v.

Lemma det_by_unit : ∀ {n} (A B X : Square n),
  WF_Matrix A → WF_Matrix B →
  WF_Unitary X → (∀ i, A × (get_vec i X) = B × (get_vec i X)) → A = B.

Definition WF_Diagonal {n : nat} (A : Square n) : Prop :=
  WF_Matrix A ∧ ∀ i j, i ≠ j → A i j = C0.

Lemma diag_Zero : ∀ n : nat, WF_Diagonal (@Zero n n).

Lemma diag_I : ∀ n : nat, WF_Diagonal (I n).

Lemma diag_I1 : WF_Diagonal (I 1).

Lemma diag_scale : ∀ {n : nat} (r : C) (A : Square n),
  WF_Diagonal A → WF_Diagonal (r .* A).

Lemma diag_plus : ∀ {n} (A B : Square n),
  WF_Diagonal A → WF_Diagonal B → WF_Diagonal (A .+ B).

Lemma diag_mult : ∀ {n : nat} (A B : Square n),
  WF_Diagonal A → WF_Diagonal B → WF_Diagonal (A × B).

Lemma div_mod_eq : ∀ (a b m : nat),
  m ≠ 0 → (a / m = b / m) → (a mod m = b mod m) → a = b.

Lemma diag_kron : ∀ {n m : nat} (A : Square n) (B : Square m),
                  WF_Diagonal A → WF_Diagonal B → WF_Diagonal (A ⊗ B).

Lemma diag_transpose : ∀ {n : nat} (A : Square n),
                     WF_Diagonal A → WF_Diagonal A⊤.

Lemma diag_adjoint : ∀ {n : nat} (A : Square n),
      WF_Diagonal A → WF_Diagonal A†.

Lemma diag_kron_n : ∀ (n : nat) {m : nat} (A : Square m),
   WF_Diagonal A → WF_Diagonal (kron_n n A).

Lemma diag_big_kron : ∀ n (l : list (Square n)),
                        (∀ A, In A l → WF_Diagonal A) →
                         WF_Diagonal (⨂ l).

Lemma diag_Mmult_n : ∀ n {m} (A : Square m),
   WF_Diagonal A → WF_Diagonal (Mmult_n n A).

Definition WF_Diagonalizable {n : nat} (A : Square n) : Prop :=
  WF_Matrix A ∧ (∃ (X X' B: Square n),
    WF_Diagonal B ∧ WF_Matrix X ∧ WF_Matrix X' ∧ X × X' = I n ∧ B = X × A × X').

Lemma diag_imps_diagble : ∀ {n} (A : Square n),
  WF_Diagonal A → WF_Diagonalizable A.

Lemma diagble_Zero : ∀ n : nat, WF_Diagonalizable (@Zero n n).

Lemma diagble_I : ∀ n : nat, WF_Diagonalizable (I n).

Lemma diagble_I1 : WF_Diagonal (I 1).

Lemma diagble_scale : ∀ {n : nat} (r : C) (A : Square n),
  WF_Diagonalizable A → WF_Diagonalizable (r .* A).

Lemma diagble_switch : ∀ {n : nat} (A B X X' : Square n),
  WF_Matrix A → WF_Matrix B → WF_Matrix X → WF_Matrix X' →
  X × X' = I n → B = X × A × X' →
  A = X' × B × X.

Fixpoint Cprod (f : nat → C) (n : nat) : C :=
  match n with
  | 0 ⇒ C1
  | S n' ⇒ (Cprod f n' × f n')%C
  end.

Lemma Cprod_0_bounded : ∀ (f : nat → C) (n : nat),
  (∃ i, i < n ∧ f i = C0) → Cprod f n = C0.

Lemma Cprod_eq_bounded : ∀ (f g : nat → C) (n : nat),
  (∀ i : nat, i < n → f i = g i) → Cprod f n = Cprod g n.

Lemma Cprod_extend_r : ∀ (f : nat → C) (n : nat),
  (Cprod f n × f n)%C = Cprod f (S n).

Lemma Cprod_extend_l : ∀ (f : nat → C) (n : nat),
  (f 0 × (Cprod (fun x ⇒ f (S x)) n))%C = Cprod f (S n).

Lemma Cprod_product : ∀ (f g h : nat → C) (n : nat),
  (∀ i, i < n → h i = (f i × g i)%C) → ((Cprod f n) × (Cprod g n))%C = Cprod h n.

Definition upper_triangular {n} (A : Square n) : Prop :=
  ∀ i j, i > j → A i j = C0.

Lemma up_tri_Zero : ∀ n : nat, upper_triangular (@Zero n n).

Lemma up_tri_I : ∀ n : nat, upper_triangular (I n).

Lemma up_tri_I1 : upper_triangular (I 1).

Lemma up_tri_scale : ∀ {n : nat} (r : C) (A : Square n),
  upper_triangular A → upper_triangular (r .* A).

Lemma up_tri_plus : ∀ {n} (A B : Square n),
  upper_triangular A → upper_triangular B → upper_triangular (A .+ B).

Lemma up_tri_mult : ∀ {n : nat} (A B : Square n),
  upper_triangular A → upper_triangular B → upper_triangular (A × B).

Lemma up_tri_reduce_0 : ∀ {n : nat} (A : Square (S n)),
  upper_triangular A → upper_triangular (reduce A 0 0).

Lemma det_up_tri_diags : ∀ {n : nat} (A : Square n),
  upper_triangular A →
  Determinant A = Cprod (fun i ⇒ A i i) n.

Definition Eigenpair {n : nat} (U : Square n) (p : Vector n × C) : Prop :=
  U × (fst p) = (snd p) .* (fst p).

Lemma all_v_eigen_I : ∀ (n : nat) (v : Vector n),
   WF_Matrix v → Eigenpair (I n) (v, C1).

Lemma diags_have_basis_eigens : ∀ (n : nat) (U : Square n) (i : nat),
  (i < n) → WF_Diagonal U → Eigenpair U (e_i i, U i i).

Local Close Scope nat_scope.

Lemma eigen_scale : ∀ {n} (A : Square n) (v : Vector n) (c1 c2 : C),
  Eigenpair A (v, c1) → Eigenpair (c2 .* A) (v, c1 × c2).

Lemma eigen_scale_div : ∀ {n} (A : Square n) (v : Vector n) (c1 c2 : C),
  c2 ≠ C0 → Eigenpair (c2 .* A) (v, Cmult c2 c1) → Eigenpair A (v, c1).

Lemma eig_unit_invertible : ∀ {n} (v : Vector n) (c : C) (X X' B : Square n),
  WF_Matrix v → WF_Matrix X → WF_Matrix X' → X' × X = I n →
  Eigenpair B (X × v, c) → Eigenpair (X' × B × X) (v, c).

Lemma eig_unit_conv : ∀ {n} (v : Vector n) (c : C) (U B : Square n),
  WF_Matrix v → WF_Unitary U →
  Eigenpair B (U × v, c) → Eigenpair (U† × B × U) (v, c).

Lemma eig_unit_norm1 : ∀ {n} (U : Square n) (c : C),
  WF_Unitary U → (∃ v, WF_Matrix v ∧ v ≠ Zero ∧ Eigenpair U (v, c)) → (c × c^* = C1)%C.

Lemma unit_has_eigen : ∀ {n} (A : Square (S n)),
  WF_Unitary A →
  ∃ (c : C) (v : Vector (S n)), Eigenpair A (v, c) ∧ v ≠ Zero ∧ WF_Matrix v.

Lemma unitary_reduction_step1 : ∀ {n} (A : Square (S n)),
  WF_Unitary A →
  ∃ X, WF_Unitary X ∧
  (∃ c : C, get_vec 0 (X†×A×X) = c .* e_i 0).

Local Open Scope nat_scope.
Lemma unitary_reduction_step2 : ∀ {n} (A : Square (S n)),
  WF_Unitary A →
  (∃ c : C, get_vec 0 A = c .* e_i 0) →
  (∀ (i j : nat), (i = 0 ∨ j = 0) ∧ i ≠ j → A i j = C0).

Lemma unitary_reduction_step3 : ∀ {n} (A : Square (S n)),
  WF_Unitary A →
  (∀ (i j : nat), (i = 0 ∨ j = 0) ∧ i ≠ j → A i j = C0) →
  ∃ (A' : Square n), WF_Unitary A' ∧ pad1 A' (A 0 0) = A.

Lemma diagble_pad1 : ∀ {n} (A : Square n) (c : C),
  WF_Diagonalizable A → WF_Diagonalizable (pad1 A c).


Theorem unit_implies_diagble : ∀ {n} (A : Square n),
  WF_Unitary A → WF_Diagonalizable A.

Definition eq_eigs {n : nat} (U1 U2 : Square n) : Prop :=
  ∀ p, WF_Matrix (fst p) → (Eigenpair U1 p → Eigenpair U2 p).

Lemma eq_eigs_implies_eq_diag : ∀ {n} (D1 D2 : Square n),
  WF_Diagonal D1 → WF_Diagonal D2 → eq_eigs D1 D2 → D1 = D2.

Lemma diagble_eigenpairs_transfer : ∀ {n} (A B X X' : Square n),
  WF_Matrix A → WF_Diagonal B → WF_Matrix X → WF_Matrix X' →
  A = X' × B × X → X × X' = I n →
  (∀ x, x < n → Eigenpair A (X' × (e_i x), B x x)).

Lemma eq_eigs_implies_eq_diagble : ∀ {n} (D1 D2 : Square n),
  WF_Diagonalizable D1 → WF_Diagonalizable D2 → eq_eigs D1 D2 → D1 = D2.

Lemma eq_eigs_implies_eq_unit : ∀ {n} (U1 U2 : Square n),
  WF_Unitary U1 → WF_Unitary U2 → eq_eigs U1 U2 → U1 = U2.

Theorem eigs_eq_gate : ∀ {n} (U1 U2 : Square n),
  WF_Unitary U1 → WF_Unitary U2 → (U1 = U2 ↔ eq_eigs U1 U2).

Local Close Scope nat_scope.

Lemma EigenXp : Eigenpair σx (∣+⟩, C1).

Lemma EigenXm : Eigenpair σx (∣-⟩, -C1).

Lemma EigenYp : Eigenpair σy (∣R⟩, C1).

Lemma EigenYm : Eigenpair σy (∣L⟩, -C1).

Lemma EigenZp : Eigenpair σz (∣0⟩, C1).

Lemma EigenZm : Eigenpair σz (∣1⟩, -C1).

Lemma EigenXXB : Eigenpair (σx ⊗ σx) (∣Φ+⟩, C1).

#[export] Hint Resolve EigenXp EigenXm EigenYp EigenYm EigenZp EigenZm EigenXXB all_v_eigen_I : eig_db.