Theory Exponentiation

section ‹Exponentiation is Diophaninte›

subsection ‹Expressing Exponentiation in terms of the alpha function›

theory Exponentiation
  imports "HOL-Library.Discrete"
begin

locale Exp_Matrices
  begin

subsubsection ‹2x2 matrices and operations›
datatype mat2 = mat (mat_11 : int) (mat_12 : int) (mat_21 : int) (mat_22 : int)
datatype vec2 = vec (vec_1: int) (vec_2: int)

fun mat_plus:: "mat2 ⇒ mat2 ⇒ mat2" where 
  "mat_plus A B = mat (mat_11 A + mat_11 B) (mat_12 A + mat_12 B) 
                      (mat_21 A + mat_21 B) (mat_22 A + mat_22 B)"

fun mat_mul:: "mat2 ⇒ mat2 ⇒ mat2" where
  "mat_mul A B = mat (mat_11 A * mat_11 B + mat_12 A * mat_21 B) 
                     (mat_11 A * mat_12 B + mat_12 A * mat_22 B) 
                     (mat_21 A * mat_11 B + mat_22 A * mat_21 B) 
                     (mat_21 A * mat_12 B + mat_22 A * mat_22 B)"

fun mat_pow:: "nat ⇒ mat2 ⇒ mat2" where
  "mat_pow 0 _ = mat 1 0 0 1" | 
  "mat_pow n A = mat_mul A (mat_pow (n - 1) A)"

lemma mat_pow_2[simp]: "mat_pow 2 A = mat_mul A A"
  by (simp add: numeral_2_eq_2)

fun mat_det::"mat2 ⇒ int" where
  "mat_det M = mat_11 M * mat_22 M - mat_12 M * mat_21 M"

fun mat_scalar_mult::"int ⇒ mat2 ⇒ mat2" where
  "mat_scalar_mult a M = mat (a * mat_11 M) (a * mat_12 M) (a * mat_21 M) (a * mat_22 M)"

fun mat_minus:: "mat2 ⇒ mat2 ⇒ mat2" where
  "mat_minus A B = mat (mat_11 A - mat_11 B) (mat_12 A - mat_12 B) 
                       (mat_21 A - mat_21 B) (mat_22 A - mat_22 B)"

fun mat_vec_mult:: "mat2 ⇒ vec2 ⇒ vec2" where
  "mat_vec_mult M v = vec (mat_11 M * vec_1 v + mat_12 M * vec_2 v) 
                          (mat_21 M * vec_1 v + mat_21 M * vec_2 v)"

definition ID :: "mat2" where "ID = mat 1 0 0 1"
declare mat_det.simps[simp del]

subsubsection ‹Properties of 2x2 matrices›
lemma mat_neutral_element: "mat_mul ID N = N" by (auto simp: ID_def)

lemma mat_associativity: "mat_mul (mat_mul D B) C = mat_mul D (mat_mul B C)"
  apply auto by algebra+ 

lemma mat_exp_law: "mat_mul (mat_pow n M) (mat_pow m M) = mat_pow (n+m) M" 
  apply (induction n, auto) by (metis mat2.sel(1,2) mat_associativity mat_mul.simps)+

lemma mat_exp_law_mult: "mat_pow (n*m) M = mat_pow n (mat_pow m M)" (is "?P n")
  apply (induction n, auto) using mat_exp_law by (metis mat_mul.simps)

lemma det_mult: "mat_det (mat_mul M1 M2) = (mat_det M1) * (mat_det M2)"
  by (auto simp: mat_det.simps algebra_simps)

subsubsection ‹Special second-order recurrent sequences›
text ‹Equation 3.2›
fun α:: "nat ⇒ nat ⇒ int" where
           "α b 0 = 0" | 
           "α b (Suc 0) = 1" |
  alpha_n: "α b (Suc (Suc n)) = (int b) * (α b (Suc n)) - (α b n)"

text ‹Equation 3.3›
lemma alpha_strictly_increasing:
  shows "int b ≥ 2 ⟹ α b n < α b (Suc n) ∧ 0 < α b (Suc n)"
proof (induct n)
  case 0
  show ?case by simp
next
  case (Suc n)
  have pos: "0 < α b (Suc n)" by (simp add:Suc)
  have "α b (Suc n) ≤ (int b) * (α b (Suc n)) - α b (Suc n)" using pos Suc by simp
  also have "... < α b (Suc (Suc n))" by (simp add: Suc)
  finally show ?case using pos Suc by simp
qed

lemma alpha_strictly_increasing_general:
  fixes b n m::nat
  assumes "b > 2 ∧ m > n"
  shows "α b m > α b n"
proof -
  from alpha_strictly_increasing assms have S2: "α b n < α b m"
    by (smt less_imp_of_nat_less lift_Suc_mono_less of_nat_0_less_iff pos2)
  show ?thesis using S2 by simp
qed


text ‹Equation 3.4›
lemma alpha_superlinear: "b > 2 ⟹ int n ≤ α b n"
  apply (induction n, auto) 
  by (smt Suc_1 alpha_strictly_increasing less_imp_of_nat_less of_nat_1 of_nat_Suc)

text ‹A simple consequence that's often useful; could also be generalized to alpha using 
      alpha linear›

lemma alpha_nonnegative:
  shows "b > 2 ⟹ α b n ≥ 0"
  using of_nat_0_le_iff alpha_superlinear order_trans by blast

text ‹Equation 3.5›
lemma alpha_linear: "α 2 n = n"
proof(induct n rule: nat_less_induct)
  case (1 n)
  have s0: "n=0 ⟹ α 2 n = n" by simp
  have s1: "n=1 ⟹ α 2 n = n" by simp
  note hyp = ‹∀m < n. α 2 m = m›
  from hyp have s2: "n>1 ⟹ α 2 (n-1) = n-1 ∧ α 2 (n-2) = n-2" by simp
  have s3: "n>1 ⟹ α 2 (Suc (Suc (n-2))) = 2*α 2 (Suc (n-2)) - α 2 (n-2)" by simp
  have s4: "n>1 ⟹ Suc (Suc (n-2)) = n" by simp
  have s5: "n>1 ⟹ Suc (n-2) = n-1" by simp
  from s3 s4 s5 have s6: "n>1 ⟹ α 2 n = 2*α 2 (n-1) - α 2 (n-2)" by simp
  from s2 s6 have s7: "n>1 ⟹ α 2 n = 2*(n-1) - (n-2)" by simp
  from s7 have s8: "n>1 ⟹ α 2 n = n" by simp
  from s0 s1 s8 show ?case by linarith
qed

text ‹Equation 3.6 (modified)›
lemma alpha_exponential_1: "b > 0 ⟹ int b ^ n ≤ α (b + 1) (n+1)"
proof(induction n)
case 0
  thus ?case by(simp)
next
  case (Suc n)
  hence "((int b)*(int b)^n) ≤ (int b)*(α (b+1) (n+1))" by simp
  hence r2: "((int b)^(Suc n)) ≤ (int (b+1))*(α (b+1) (n+1)) - (α (b+1) (n+1))" 
    by (simp add: algebra_simps)
  have "(int b+1) *(α (b+1) (n+1)) - (α (b+1) (n+1)) ≤ (int b+1)*(α (b+1) (n+1)) - α (b+1) n" 
    using alpha_strictly_increasing Suc by (smt Suc_eq_plus1 of_nat_0_less_iff of_nat_Suc)
  thus ?case using r2 by auto
qed

lemma alpha_exponential_2: "int b>2 ⟹ α b (n+1) ≤ (int b)^(n)"
proof(induction n)
  case 0
  thus ?case by simp
next
  case (Suc n)
  hence s1: "α b (n+2) ≤ (int b)^(n+1) - α b n" by simp
  have "(int b)^(n+1) - (α b n) ≤ (int b)^(n+1)" 
    using alpha_strictly_increasing Suc by (smt α.simps(1) alpha_superlinear of_nat_1 of_nat_add 
                                            of_nat_le_0_iff of_nat_less_iff one_add_one)
  thus ?case using s1 by simp
qed

subsubsection ‹First order relation›
text ‹Equation 3.7 - Definition of A›
fun A :: "nat ⇒ nat ⇒ mat2" where
  "A b 0 = mat 1 0 0 1" | 
  A_n: "A b n = mat (α b (n + 1)) (-(α b n)) (α b n) (-(α b (n - 1)))"

text ‹Equation 3.9 - Definition of B›
fun B :: "nat ⇒ mat2" where
  "B b = mat (int b) (-1) 1 0"

declare A.simps[simp del]
declare B.simps[simp del]

text ‹Equation 3.8›
lemma A_rec: "b>2 ⟹ A b (Suc n) = mat_mul (A b n) (B b)" 
  by (induction n, auto simp: A.simps B.simps)

text ‹Equation 3.10›
lemma A_pow: "b>2 ⟹ A b n = mat_pow n (B b)"
  apply (induction n, auto simp: A.simps B.simps)
    subgoal by (smt A.elims Suc_eq_plus1 α.simps α.simps(2) mat2.sel)
    subgoal for n apply (cases "n=0", auto) 
      using A.simps(2)[of b "n-1"] gr0_conv_Suc mult.commute by auto
    subgoal by (metis A.simps(2) Suc_eq_plus1 α.simps(2) mat2.sel(1) mat_pow.elims)
    subgoal by (metis A.simps(2) α.simps(1) add.inverse_neutral mat2.sel(2) mat_pow.elims)
    done


subsubsection ‹Characteristic equation›
text ‹Equation 3.11›
lemma A_det: "b>2 ⟹ mat_det (A b n) = 1"
  apply (auto simp: A_pow, induction n, simp add: mat_det.simps)
  using det_mult apply (auto simp del: mat_mul.simps) by (simp add: B.simps mat_det.simps)

text ‹Equation 3.12›
lemma alpha_det1:
  assumes "b>2"
  shows "(α b (Suc n))^2 - (int b) * α b (Suc n) * α b n + (α b n)^2 = 1"
proof(cases "n = 0")
  case True
  thus ?thesis by auto
next
  case False
  hence "A b n = mat (α b (n + 1)) (-(α b n)) (α b n) (-(α b (n - 1)))" using A.elims neq0_conv by blast
  hence "mat_det (A b n)  = (α b n)^2 - (α b (Suc n)) * α b (n-1)" 
    apply (auto simp: mat_det.simps) by (simp add: power2_eq_square)
  moreover hence "... =  (α b (Suc n))^2 - b * (α b (Suc n)) * α b n + (α b n)^2"
    using False alpha_n[of b "n-1"] apply(auto simp add: algebra_simps)
    by (metis Suc_1 distrib_left mult.commute mult.left_commute power_Suc power_one_right)
  ultimately show ?thesis using A_det assms by auto
qed

text ‹Equation 3.12›
lemma alpha_det2:
  assumes "b>2" "n>0"
  shows "(α b (n-1))^2 - (int b) * (α b (n-1) * (α b n)) + (α b n)^2 = 1"
  using alpha_det1 assms by (smt One_nat_def Suc_diff_Suc diff_zero mult.commute mult.left_commute)

text ‹Equations 3.14 to 3.17›
lemma alpha_char_eq:
  fixes x y b:: nat
  shows  "(y < x ∧ x * x + y * y = 1 + b * x * y) ⟹ (∃m. int y = α b m ∧ int x = α b (Suc m))"
proof (induct y arbitrary: x rule:nat_less_induct)
  case (1 n)

  note pre = ‹n < x ∧ (x * x + n * n = 1 + b * x * n)›

  have h0: "int (x * x + n * n) = int (x * x) + int (n * n)" by simp
  from pre h0 have pre1: "int x * int x + int(n * n) = int 1 + int(b * x * n)" by simp

  have i0: "int (n * n) = int n * int n" by simp
  have i1: "int (b * x * n) = int b * int x * int n" by simp
  from pre1 i0 i1 have pre2: "int x * int x + int n * int n = 1 + int b * int x * int n" by simp

  from pre2 have j0: "int n * int n - 1 = int b * int x * int n - int x * int x" by simp
  have j1:"… = int x * (int b * int n - int x)" by (simp add: right_diff_distrib)
  from j0 j1 have pre3:"int n * int n - 1 = int x * (int b * int n - int x)" by simp

  have k0: "int n * int n - 1 < int n * int n" by simp
  from pre3 k0 have k1:"int n * int n >  int x * (int b * int n - int x)" by simp
  from pre have k2: "int n ≤ int x" by simp
  from k2 have k3: "int x * int n ≥ int n * int n" by (simp add: mult_mono)
  from k1 k3 have k4: "int x * int n > int x * (int b * int n - int x)" by linarith
  from pre k4 have k5: "int n > int b * int n - int x" by simp

  from pre have l0:"n = 0 ⟹ x = 1" by simp
  from l0 have l1: "n = 0 ⟹ x = Suc 0" by simp
  from l1 have l2: "n = 0 ⟹ int n = α b 0 ∧ int x = α b (Suc 0)" by simp
  from l2 have l3: "n = 0 ⟹ ∃m. int n = α b m ∧ int x = α b (Suc m)" by blast

  have m0: "n > 0 ⟹ int n * int n - 1 ≥ 0" by simp
  from pre3 m0 have m1: "n > 0 ⟹ int x * (int b * int n - int x) ≥ 0" by simp
  from m1 have m2: "n > 0 ⟹ int b * int n - int x ≥ 0" using zero_le_mult_iff by force

  from j0 have n0: "int x * int x - int b * int x * int n + int n * int n = 1" by simp
  have n1: "(int b * int n - int x) * (int b * int n - int x) = int b * int n * (int b * int n - int x) - int x * (int b * int n - int x)" by (simp add: left_diff_distrib)
  from n1 have n2: "int n * int n - int b * int n * (int b * int n - int x) + (int b * int n - int x) * (int b * int n - int x) = int n * int n - int x * (int b * int n - int x)" by simp
  from n0 n2 j1 have n3: "int n * int n - int b * int n * (int b * int n - int x) + (int b * int n - int x) * (int b * int n - int x) = 1" by linarith
  from n3 have n4: "int n * int n + (int b * int n - int x) * (int b * int n - int x) = 1 + int b * int n * (int b * int n - int x)" by simp
  have n5: "int b * int n = int (b * n)" by simp
  from n5 m2 have n6: "n > 0 ⟹ int b * int n - int x = int (b * n - x)" by linarith
  from n4 n6 have n7: "n > 0 ⟹ int (n * n + (b * n - x) * (b * n - x)) = int (1 + b * n * (b * n - x))" by simp
  from n7 have n8: "n > 0 ⟹ n * n + (b * n - x) * (b * n - x) = 1 + b * n * (b * n - x)" using of_nat_eq_iff by blast

  note hyp = ‹∀m<n. ∀x. m < x ∧ x * x + m * m = 1 + b * x * m ⟶
          (∃ma. int m = α b ma ∧ int x = α b (Suc ma))›

  from k5 n6 n8 have o0: "n > 0 ⟹ (b * n - x) < n ∧ n * n + (b * n - x) * (b * n - x) = 1 + b * n * (b * n - x)" by simp
  from o0 hyp have o1: "n > 0 ⟹ (∃ma. int (b * n - x) = α b ma ∧ int n = α b (Suc ma))" by simp

  from o1 l3 n6 show ?case by force
qed

lemma alpha_char_eq2:
  assumes "(x*x + y*y = 1 + b * x * y)" "b>2"
  shows  "(∃n. int x = α b n)"
proof -
  have "x ≠ y"
  proof(rule ccontr, auto)
    assume "x=y"
    hence "2*x*x = 1+b*x*x" using assms by simp
    hence "2*x*x ≥ 1+2*x*x" using assms by (metis add_le_mono le_less mult_le_mono1)
    thus False by auto
  qed
  thus ?thesis
  proof(cases "x<y")
    case True
    hence "∃n. int x = α b n ∧ int y = α b (Suc n)" using alpha_char_eq assms
      by (simp add: add.commute power2_eq_square)
    thus ?thesis by auto
  next
    case False
    hence "∃j. int y = α b j ∧ int x = α b (Suc j)" using alpha_char_eq assms ‹x ≠ y› by auto
    thus ?thesis by blast
  qed
qed


subsubsection ‹Divisibility properties›
text ‹The following lemmas are needed in the proof of equation 3.25›
lemma representation:
  fixes k m :: nat
  assumes "k > 0" "n = m mod k" "l = (m-n)div k"
  shows "m = n+k*l ∧ 0≤n ∧ n≤k-1" by (metis Suc_pred' assms le_add2 le_add_same_cancel2 
      less_Suc_eq_le minus_mod_eq_div_mult minus_mod_eq_mult_div mod_div_mult_eq
      mod_less_divisor neq0_conv nonzero_mult_div_cancel_left)

lemma div_3251:
  fixes b k m:: nat
  assumes "b>2" and "k>0"
  defines "n ≡ m mod k"
  defines "l ≡ (m-n) div k"
  shows "A b m = mat_mul (A b n) (mat_pow l (A b k))"
proof -
  from assms(2) l_def n_def representation have m: "m = n+k*l ∧ 0≤n ∧ n≤k-1" by simp
  from A_pow assms(1) have Abm2: "A b m = mat_pow m (B b)" by simp
  from m have Bm: "mat_pow m (B b) =mat_pow (n+k*l) (B b)" by simp
  from mat_exp_law have as1: "mat_pow (n+k*l) (B b) 
                            = mat_mul (mat_pow n (B b)) (mat_pow (k*l) (B b))" by simp
  from mat_exp_law_mult have as2: "mat_pow (k*l) (B b) = mat_pow l (mat_pow k (B b))" 
    by (metis mult.commute)
  from A_pow assms have Abn: "mat_pow n (B b) = A b n" by simp
  from A_pow assms(1) have Ablk: "mat_pow l (mat_pow k (B b)) = mat_pow l (A b k)" by simp
  from Ablk Abm2 Abn Bm as1 as2 show Abm: "A b m = mat_mul (A b n) (mat_pow l (A b k))" by simp
qed

lemma div_3252:
  fixes a b c d m :: int and l :: nat
  defines "M ≡ mat a b c d"
  assumes "mat_21 M mod m = 0"
  shows "(mat_21 (mat_pow l M)) mod m = 0 " (is "?P l")
proof(induction l)
  show "?P 0" by simp
next
  fix l assume IH: "?P l"
  define Ml where "Ml = mat_pow l M"
  have S1: "mat_pow (Suc(l)) M = mat_mul M (mat_pow l M)" by simp
  have S2: "mat_21 (mat_mul M Ml) = mat_21 M * mat_11 Ml + mat_22 M * mat_21 Ml" 
    by (rule_tac mat_mul.induct mat_plus.induct, auto)
  have S3: "mat_21 (mat_pow (Suc(l)) M) = mat_21 M * mat_11 Ml + mat_22 M * mat_21 Ml" 
    using S1 S2 Ml_def by simp
  from assms(2) have S4: "(mat_21 M * mat_11 Ml) mod m = 0" by auto
  from IH Ml_def have S5: " mat_22 M * mat_21 Ml mod m = 0" by auto
  from S4 S5 have S6: "(mat_21 M * mat_11 Ml + mat_22 M * mat_21 Ml) mod m = 0" by auto
  from S3 S6 show "?P (Suc(l))" by simp
qed

lemma div_3253:
  fixes a b c d m:: int and l :: nat
  defines "M ≡ mat a b c d"
  assumes "mat_21 M mod m = 0"
  shows "((mat_11 (mat_pow l M)) - a^l) mod m = 0" (is "?P l")
proof(induction l)
  show "?P 0" by simp
next
  fix l assume IH: "?P l"
  define Ml where "Ml = mat_pow l M"
  from Ml_def have S1: "mat_pow (Suc(l)) M = mat_mul M Ml" by simp
  have S2: "mat_11 (mat_mul M Ml) = mat_11 M * mat_11 Ml + mat_12 M * mat_21 Ml" 
    by (rule_tac mat_mul.induct mat_plus.induct, auto)
  hence S3: "mat_11 (mat_pow (Suc(l)) M) = mat_11 M * mat_11 Ml + mat_12 M * mat_21 Ml" 
    using S1 by simp
  from M_def Ml_def assms(2) div_3252 have S4: "mat_21 Ml mod m = 0" by auto
  from IH Ml_def have S5: "(mat_11 Ml - a ^ l) mod m = 0" by auto
  from IH M_def have S6: "(mat_11 M -a) mod m = 0" by simp
  from S4 have S7: "(mat_12 M * mat_21 Ml) mod m = 0" by auto
  from S5 S6 have S8: "(mat_11 M * mat_11 Ml- a^(Suc(l))) mod m = 0" 
  by (metis M_def mat2.sel(1) mod_0 mod_mult_right_eq mult_zero_right power_Suc right_diff_distrib)
  have S9: "(mat_11 M * mat_11 Ml - a^(Suc(l)) + mat_12 M * mat_21 Ml ) mod m = 0" 
    using S7 S8 by auto
  from S9 have S10: "(mat_11 M * mat_11 Ml + mat_12 M * mat_21 Ml - a^(Suc(l))) mod m = 0" by smt
  from S3 S10 show "?P (Suc(l))" by auto
qed

text ‹Equation 3.25›
lemma divisibility_lemma1:
    fixes b k m:: nat
  assumes "b>2" and "k>0"
  defines "n ≡ m mod k"
  defines "l ≡ (m-n) div k"
  shows "α b m mod α b k = α b n * (α b (k+1)) ^ l mod α b k"
proof -
  from assms(2) l_def n_def representation have m: "m = n+k*l ∧ 0≤n ∧ n≤k-1" by simp
  consider (eq0) "n = 0" | (neq0) "n > 0" by auto
  thus ?thesis
  proof cases
    case eq0
    have Abm_gen: "A b m = mat_mul (A b n) (mat_pow l (A b k))" 
      using assms div_3251 l_def n_def by blast
    have Abk: "mat_pow l (A b k) = mat_pow l (mat (α b (k+1)) (-α b k) (α b k) (-α b (k-1)))" 
      using assms(2) neq0_conv by (metis A.elims)
    from eq0 have Abm: "A b m = mat_pow l (mat (α b (k+1)) (-α b k) (α b k) (-α b (k-1)))" 
      using A_pow ‹b>2› apply (auto simp: A.simps B.simps)
      by (metis Abk Suc_eq_plus1 add.left_neutral m mat_exp_law_mult mult.commute)
    have Abm1: "mat_21 (A b m) = α b m" by (metis A.elims α.simps(1) mat2.sel(3))
    have Abm2: "mat_21 (mat_pow l (mat (α b (k+1)) (-α b k) (α b k) (-α b (k-1)))) mod (α b k) = 0" 
      using Abm div_3252 by simp
    from Abm Abm1 Abm2 have MR0: "α b m mod α b k = 0" by simp
    from MR0 eq0 show ?thesis by simp
  next case neq0
    from assms have Abm_gen: "A b m = mat_mul (A b n) (mat_pow l (A b k))" 
      using div_3251 l_def n_def by blast
    from assms(2) neq0_conv have Abk: "mat_pow l (A b k) 
               = mat_pow l (mat (α b (k+1)) (-α b k) (α b k) (-α b (k-1)))" by (metis A.elims)
  from n_def neq0 have N0: "n>0" by simp
  define M where "M = mat (α b (n + 1)) (-(α b n)) (α b n) (-(α b (n - 1)))"
  define N where "N = mat_pow l (mat (α b (k+1)) (-α b k) (α b k) (-α b (k-1)))"
  from Suc_pred' neq0 have Abn: "A b n = mat (α b (n + 1)) (-(α b n)) (α b n) (-(α b (n - 1)))" 
    by (metis A.elims neq0_conv)
  from Abm_gen Abn Abk M_def N_def have Abm: "A b m = mat_mul M N" by simp
  (* substitutions done! next: calculations *)
  from Abm have S1: "mat_21 (mat_mul M N) = mat_21 M * mat_11 N + mat_22 M * mat_21 N" 
    by (rule_tac mat_mul.induct mat_plus.induct, auto)
  have S2: "mat_21 (A b m) = α b m" by (metis A.elims α.simps(1) mat2.sel(3))
  from S1 S2 Abm have S3: "α b m = mat_21 M * mat_11 N + mat_22 M * mat_21 N" by simp
  from S3 have S4: "(α b m - (mat_21 M * mat_11 N + mat_22 M * mat_21 N)) mod (α b k) = 0" by simp
  from M_def have S5: "mat_21 M = α b n" by simp
  from div_3253 N_def have S6: "(mat_11 N - (α b (k+1)) ^ l) mod (α b k) = 0" by simp
  from N_def Abm div_3252 have S7: "mat_21 N mod (α b k) = 0" by simp
  from S4 S7 have S8: "(α b m - mat_21 M * mat_11 N) mod (α b k) = 0" by algebra
  from S5 S6 have S9: "(mat_21 M * mat_11 N - (α b n) * (α b (k+1)) ^ l) mod (α b k) = 0"
    by (metis mod_0 mod_mult_left_eq mult.commute mult_zero_left right_diff_distrib')
  from S8 S9 show ?thesis
  proof -
    have "(mat_21 M * mat_11 N - α b m) mod α b k = 0"
      using S8 by presburger
    hence "∀i. (α b m - (mat_21 M * mat_11 N - i)) mod α b k = i mod α b k"
      by (metis (no_types) add.commute diff_0_right diff_diff_eq2 mod_diff_right_eq)
    thus ?thesis
      by (metis (no_types) S9 diff_0_right mod_diff_right_eq)
  qed
    qed
  qed

text ‹Prerequisite lemma for 3.27›
lemma div_coprime:
  assumes "b>2" "n ≥ 0"
    shows "coprime (α b k) (α b (k+1))" (is "?P")
proof(rule ccontr)
  assume as: "¬ ?P"
  define n where "n = gcd (α b k) (α b (k+1))"
  from n_def have S1: "n > 1"
    using alpha_det1 as assms(1) coprime_iff_gcd_eq_1 gcd_pos_int right_diff_distrib' 
    by (smt add.commute plus_1_eq_Suc)
  have S2: "(α b (Suc k))^2 - (int b) * α b (Suc k) * (α b k) + (α b k)^2 = 1" 
    using alpha_det1 assms by auto
  from n_def have D1: " n dvd (α b (k+1))^2" by (simp add: numeral_2_eq_2)
  from n_def have D2: " n dvd (- (int b) * α b (k+1) * (α b k))" by simp
  from n_def have D3: "n dvd (α b k)^2" by (simp add: gcd_dvdI1)
  have S3: "n dvd ((α b (Suc k))^2 - (int b) * α b (Suc k) * (α b k) + (α b k)^2)" 
    using D1 D2 D3 by simp
  from S2 S3 have S4: "n dvd 1" by simp
  from S4 n_def as is_unit_gcd show "False" by blast
qed

text ‹Equation 3.27›
lemma divisibility_lemma2:
  fixes b k m:: nat
  assumes "b>2" and "k>0"
  defines "n ≡ m mod k"
  defines "l ≡ (m-n) div k"
  assumes "α b k dvd α b m"
    shows "α b k dvd α b n"
proof -
  from assms(2) l_def n_def representation have m: "0≤n ∧ n≤k-1" by simp
  from divisibility_lemma1 assms(1) assms(2) l_def n_def have S1: 
    "(α b m) mod (α b k)  = (α b n) * (α b (k+1)) ^ l mod (α b k)" by blast
  from S1 assms(5) have S2: "(α b k) dvd ((α b n) * (α b (k+1)) ^ l)" by auto
  show ?thesis
    using S1 div_coprime S2 assms(1) apply auto
    using coprime_dvd_mult_left_iff coprime_power_right_iff by blast  
qed

text ‹Equation 3.23 - main result of this section›
theorem divisibility_alpha:
  assumes "b>2" and "k>0"
    shows "α b k dvd α b m ⟷ k dvd m" (is "?P ⟷ ?Q")
proof
  assume Q: "?Q"
  define n where "n = m mod k"
  have N: "n=0" by (simp add: Q n_def)
  from N have Abn: "α b n = 0" by simp
  from Abn divisibility_lemma1 assms(1) assms(2) mult_eq_0_iff n_def show "?P"
    by (metis dvd_0_right dvd_imp_mod_0 mod_0_imp_dvd) 
next 
  assume P: "?P"
  define n where "n = m mod k"
  define l where "l = (m-n) div k"
  define B where "B = (mat (int b) (-1) 1 0)"
  have S1: "(α b n) mod (α b k) = 0" 
    using divisibility_lemma2 assms(1) assms(2) n_def P by simp
  from n_def assms(2) have m: "n < k" using mod_less_divisor by blast
  from alpha_strictly_increasing m assms(1) have S2: "α b n < α b k"
    by (smt less_imp_of_nat_less lift_Suc_mono_less of_nat_0_less_iff pos2)
  from S1 S2 have S3: "n=0"
    by (smt alpha_superlinear assms(1) mod_pos_pos_trivial neq0_conv of_nat_0_less_iff)
  from S3 n_def show "?Q" by auto
qed

subsubsection ‹Divisibility properties (continued)›

text ‹Equation 3.28 - main result of this section›

lemma divisibility_equations:
  assumes 0: "m = k*l" and "b>2" "m>0"
  shows  "A b m = mat_pow l (mat_minus (mat_scalar_mult (α b k) (B b)) 
                                        (mat_scalar_mult (α b (k-1)) ID))"
    apply (auto simp del: mat_pow.simps mat_mul.simps mat_minus.simps mat_scalar_mult.simps
                simp add: A_pow mult.commute[of k l] assms mat_exp_law_mult)
    using A_pow[of b k] ‹m>0›
    apply (auto simp: A.simps ‹m>0› ID_def B.simps) 
    using A.simps(2) alpha_n One_nat_def Suc_eq_plus1 Suc_pred assms ‹m>0› assms
        mult.commute nat_0_less_mult_iff
    by (smt mat_exp_law_mult)

lemma divisibility_cong:
  fixes e f :: int
  fixes l :: nat
  fixes M :: mat2
  assumes "mat_22 M = 0" "mat_21 M = 1"
  shows "(mat_21 (mat_pow l (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID)))) mod e^2 = (-1)^(l-1)*l*e*f^(l-1)*(mat_21 M) mod e^2
        ∧ mat_22  (mat_pow l (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))) mod e^2 = (-1)^l *f^l mod e^2" 
(is "?P l ∧ ?Q l" )
proof(induction l)
  case 0
  then show ?case by simp
next
  case (Suc l)
  have S2: "mat_pow (Suc(l)) (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID)) =
    mat_mul (mat_pow l (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))) (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))"
    using mat_exp_law[of l _ 1] mat2.sel by (auto, metis)+
  define a1 where "a1 = mat_11 (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))"
  define b1 where "b1 = mat_12 (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))"
  define c1 where "c1 = mat_21 (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))"
  define d1 where "d1 = mat_22 (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))"
  define a where "a = mat_11 M"
  define b where "b = mat_12 M"
  define c where "c = mat_21 M"
  define d where "d = mat_22 M"
  define g where "g = mat_21 (mat_pow l (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID)))"
  define h where "h = mat_22 (mat_pow l (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID)))"
  from S2 g_def a1_def h_def c1_def have S3: "mat_21 (mat_pow (Suc(l)) (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID))) = g*a1 + h*c1"
    by simp
  from S2 g_def b1_def h_def d1_def have S4: "mat_22 (mat_pow (Suc(l)) (mat_minus (mat_scalar_mult e M) (mat_scalar_mult f ID)))
    = g*b1+h*d1" by simp
  have S5: "mat_11 (mat_scalar_mult e M) = e*a" by (simp add: a_def)
  have S6: "mat_12 (mat_scalar_mult e M) = e*b" by (simp add: b_def)
  have S7: "mat_21 (mat_scalar_mult e M) = e*c" by (simp add: c_def)
  have S8: "mat_22 (mat_scalar_mult e M) = e*d" by (simp add: d_def)
  from a1_def S5  have S9: "a1 = e*a-f"by (simp add: Exp_Matrices.ID_def)
  from b1_def S6  have S10: "b1 = e*b" by (simp add: Exp_Matrices.ID_def)
  from c1_def S7  have S11: "c1 = e*c" by (simp add: Exp_Matrices.ID_def)
  from S11 assms(2) c_def have S115: "c1 = e" by simp
  from d1_def S8  have S12: "d1 = e*d - f" by (simp add: Exp_Matrices.ID_def)
  from S12 assms(1) d_def have S125: "d1 = - f" by simp
  from assms(2) c_def Suc g_def c_def have S13: "g mod e^2 = (-1)^(l-1)*l*e*f^(l-1)*c mod e^2" by blast
  from assms(2) c_def S13 have S135: "g mod e^2 = (-1)^(l-1)*l*e*f^(l-1) mod e^2" by simp
  from Suc h_def have S14: "h mod e^2 = (-1)^l *f^l mod e^2" by simp
  from S10 S135 have S27: "g*b1 mod e^2 = (-1)^(l-1)*l*e*f^(l-1)*e*b mod e^2" by (metis mod_mult_left_eq mult.assoc)
  from S27 have S28: "g*b1 mod e^2 = 0 mod e^2" by (simp add: power2_eq_square)
  from S125 S14 mod_mult_cong have S29: "h*d1 mod e^2 = (-1)^l *f^l*(- f) mod e^2" by blast
  from S29 have S30: "h*d1 mod e^2 = (-1)^(l+1) *f^l*f mod e^2" by simp
  from S30 have S31: "h*d1 mod e^2 = (-1)^(l+1) *f^(l+1) mod e^2" by (metis mult.assoc power_add power_one_right)
  from S31 have F2: "?Q (Suc(l))" by (metis S28 S4 Suc_eq_plus1 add.left_neutral mod_add_cong)
  from S9 S13 have S15: "g*a1 mod e^2 = ((-1)^(l-1)*l*e*f^(l-1)*c*(e*a-f))mod e^2" by (metis mod_mult_left_eq)
  have S16: "((-1)^(l-1)*l*e*f^(l-1)*c*(e*a-f)) = ((-1)^(l-1)*l*e^2*f^(l-1)*c*a) - f*(-1)^(l-1)*l*e*f^(l-1)*c" by algebra
  have S17: "((-1)^(l-1)*l*e^2*f^(l-1)*c*a) mod e^2 = 0 mod e^2" by simp
  from S17 have S18: "(((-1)^(l-1)*l*e^2*f^(l-1)*c*a) - f*(-1)^(l-1)*l*e*f^(l-1)*c) mod e^2 =
    - f*(-1)^(l-1)*l*e*f^(l-1)*c mod e^2"
    proof -
      have f1: "∀i ia. (ia::int) - (0 - i) = ia + i"
        by auto
      have "∀i ia. ((0::int) - ia) * i = 0 - ia * i"
        by auto
      then show ?thesis using f1
      proof -
        have f1: "⋀i. (0::int) - i = - i"
          by presburger
        then have "⋀i. (i - - ((- 1) ^ (l - 1) * int l * e⇧2 * f ^ (l - 1) * c * a)) mod e⇧2 = i mod e⇧2"
          by (metis (no_types) S17 ‹∀i ia. ia - (0 - i) = ia + i› add.right_neutral mod_add_right_eq)
        then have "⋀i. ((- 1) ^ (l - 1) * int l * e⇧2 * f ^ (l - 1) * c * a - i) mod e⇧2 = - i mod e⇧2"
          using f1 by (metis ‹∀i ia. ia - (0 - i) = ia + i› uminus_add_conv_diff)
        then show ?thesis
          using f1 ‹∀i ia. (0 - ia) * i = 0 - ia * i› by presburger
      qed
  qed
  from S15 S16 S18 have S19: "g*a1 mod e^2 = - f*(-1)^(l-1)*l*e*f^(l-1)*c mod e^2" by presburger
  from S11 S14 have S20: "h*c1 mod e^2 = (-1)^l *f^l*e*c mod e^2" by (metis mod_mult_left_eq mult.assoc)
  from S19 S20 have S21: "(g*a1 + h*c1) mod e^2 = (- f*(-1)^(l-1)*l*e*f^(l-1)*c + (-1)^l *f^l*e*c) mod e^2" using mod_add_cong by blast
  from assms(2) c_def have S22: "(- f*(-1)^(l-1)*l*e*f^(l-1)*c + (-1)^l *f^l*e*c) mod e^2=(- f*(-1)^(l-1)*l*e*f^(l-1) + (-1)^l *f^l*e) mod e^2" by simp
  have S23: "(- f*(-1)^(l-1)*l*e*f^(l-1) + (-1)^l *f^l*e) mod e^2 = (f*(-1)^(l)*l*e*f^(l-1) + (-1)^l *f^l*e) mod e^2"
    by (smt One_nat_def Suc_pred mult.commute mult_cancel_left2 mult_minus_left neq0_conv of_nat_eq_0_iff power.simps(2))
  have S24: "(f*(-1)^(l)*l*e*f^(l-1) + (-1)^l *f^l*e) mod e^2 = ((-1)^(l)*l*e*f^l + (-1)^l *f^l*e) mod e^2"
    by (smt One_nat_def Suc_pred mult.assoc mult.commute mult_eq_0_iff neq0_conv of_nat_eq_0_iff power.simps(2))
  have S25: "((-1)^(l)*l*e*f^l + (-1)^l *f^l*e) mod e^2 = ((-1)^(l)*(l+1)*e*f^l) mod e^2"
  proof -
    have f1: "∀i ia. (ia::int) * i = i * ia"
      by simp
    then have f2: "∀i ia. (ia::int) * i - - i = i * (ia - - 1)"
      by (metis (no_types) mult.right_neutral mult_minus_left right_diff_distrib')
    have "∀n. int n - - 1 = int (n + 1)"
      by simp
    then have "e * (f ^ l * (int l * (- 1) ^ l - - ((- 1) ^ l))) mod e⇧2 = e * (f ^ l * ((- 1) ^ l * int (l + 1))) mod e⇧2"
      using f2 by presburger
    then have "((- 1) ^ l * int l * e * f ^ l - - ((- 1) ^ l) * f ^ l * e) mod e⇧2 = (- 1) ^ l * int (l + 1) * e * f ^ l mod e⇧2"
      using f1
    proof -
      have f1: "⋀i ia ib. (i::int) * (ia * ib) = ia * (i * ib)"
        by simp
      then have "⋀i ia ib. (i::int) * (ia * ib) - - (i * ib) = (ia - - 1) * (i * ib)"
        by (metis (no_types) ‹∀i ia. ia * i = i * ia› f2)
      then show ?thesis
        using f1 by (metis (no_types) ‹∀i ia. ia * i = i * ia› ‹e * (f ^ l * (int l * (- 1) ^ l - - ((- 1) ^ l))) mod e⇧2 = e * (f ^ l * ((- 1) ^ l * int (l + 1))) mod e⇧2› f2 mult_minus_right)
    qed
    then show ?thesis
      by simp
  qed
  from S21 S22 S23 S24 S25 have S26: "(g*a1 + h*c1) mod e^2 = ((-1)^(l)*(l+1)*e*f^l) mod e^2" by presburger
  from S3 S26 have F1: "?P (Suc(l))" by (metis Suc_eq_plus1 assms(2) diff_Suc_1 mult.right_neutral)
  from F1 F2 show ?case by simp
qed

lemma divisibility_congruence:
  assumes "m = k*l" and "b>2" "m>0"
  shows "α b m mod (α b k)^2 = ((-1)^(l-1)*l*(α b k)*(α b (k-1))^(l-1)) mod (α b k)^2"
proof -
  have S0: "α b m = mat_21 (A b m)" by (metis A.elims assms(3) mat2.sel(3) neq0_conv)
  from assms S0 divisibility_equations have S1: "α b m =
    mat_21 ( mat_pow l (mat_minus (mat_scalar_mult (α b k) (B b)) 
                                        (mat_scalar_mult (α b (k-1)) ID)))" by auto
  have S2: "mat_21 (B b) = 1" using Binomial.binomial_ring by (simp add: Exp_Matrices.B.simps)
  have S3: "mat_22 (B b) = 0"  by (simp add: Exp_Matrices.B.simps)
  from S1 S2 S3 divisibility_cong show ?thesis by (metis mult.right_neutral)
qed

text ‹ Main result section 3.5 ›
theorem divisibility_alpha2:
  assumes "b>2" "m>0"
  shows "(α b k)^2 dvd (α b m) ⟷ k*(α b k) dvd m" (is "?P ⟷ ?Q")
proof
  assume Q: "?Q"
  then show "?P"
  proof(cases "k dvd m")
    case True
    then obtain l where mkl: "m = k * l" by blast
    from Q assms mkl have S0: "l mod α b k = 0" by simp
    from S0 have S1: "l*(α b k) mod (α b k)^2 = 0" by (simp add: power2_eq_square)
    from S1 have S2: "((-1)^(l-1)*l*(α b k)*(α b (k-1))^(l-1)) mod (α b k)^2 = 0"
      proof -
        have "∀i. α b k * (int l * i) mod (α b k)⇧2 = 0"
          by (metis (no_types) S1 mod_0 mod_mult_left_eq mult.assoc mult.left_commute mult_zero_left)
        then show ?thesis
          by (simp add: mult.assoc mult.left_commute)
    qed
    from assms divisibility_congruence mkl have S3: 
      "α b m mod (α b k)^2 = ((-1)^(l-1)*l*(α b k)*(α b (k-1))^(l-1)) mod (α b k)^2" by simp
    from S2 S3 have S4: "α b m mod (α b k)^2 = 0" by linarith
    then show ?thesis by auto
  next
    case False
    then show ?thesis using Q dvd_mult_left int_dvd_int_iff by blast
  qed
next 
  assume P: "?P"
  show "?Q" 
  proof(cases "k dvd m")
    case True
      then obtain l where mkl: "m = k * l" by blast
      from assms mkl divisibility_congruence have S0: 
          "α b m mod (α b k)^2 = ((-1)^(l-1)*l*(α b k)*(α b (k-1))^(l-1)) mod (α b k)^2" by simp
      from S0 P have S1: "(α b k)^2 dvd ((-1)^(l-1)*l*(α b k)*(α b (k-1))^(l-1))" by auto
      from S1 have S2: "(α b k)^2 dvd l*(α b k)*(α b (k-1))^(l-1)"
        by (metis (no_types, opaque_lifting) Groups.mult_ac(1) dvd_trans dvd_triv_right left_minus_one_mult_self)
      from S2 have S3: "(α b k) dvd l*(α b (k-1))^(l-1)"
        by (metis (full_types) Exp_Matrices.alpha_superlinear assms(1) assms(2) mkl 
            mult.assoc mult.commute mult_0 not_less_zero of_nat_le_0_iff power2_eq_square zdvd_mult_cancel)
      from div_coprime Suc_eq_plus1 Suc_pred' assms(1) assms(2) mkl less_imp_le_nat nat_0_less_mult_iff
      have S4: "coprime (α b k) (α b (k-1))" by (metis coprime_commute)
      hence "coprime (α b k) ((α b (k-1))^(l-1))" using coprime_power_right_iff by blast
      hence "(α b k) dvd l" using S3 using coprime_dvd_mult_left_iff by blast
      then show ?thesis by (simp add: mkl)
  next
    case False
    then show ?thesis 
      apply(cases "0<k")
      subgoal using divisibility_alpha[of b k m] assms using dvd_mult_left P by auto
      subgoal using Exp_Matrices.alpha_strictly_increasing_general assms(1) P by fastforce
    done
  qed
qed

subsubsection ‹Congruence properties›
text ‹In this section we will need the inverse matrices of A and B›
fun A_inv :: "nat ⇒ nat ⇒ mat2" where
  "A_inv b n = mat (-α b (n-1)) (α b n) (-α b n) (α b (n+1))"

fun B_inv :: "nat ⇒ mat2" where
  "B_inv b = mat 0 1 (-1) b"

lemma A_inv_aux: "b>2 ⟹ n>0 ⟹ α b n * α b n - α b (Suc n) * α b (n - Suc 0) = 1" 
  apply (induction n, auto) subgoal for n using alpha_det1[of b n] apply auto by algebra done

lemma A_inverse[simp]: "b>2 ⟹ n>0 ⟹ mat_mul (A_inv b n) (A b n) = ID"
  using mat2.expand[of "mat_mul (A_inv b n) (A b n)" ID] apply rule
  using ID_def A.simps(2)[of _ "n-1"] ID_def apply (auto) 
    subgoal using mat2.sel(1)[of 1 0 0 1] apply (auto) 
      using A_inv_aux[of b n] by (auto simp: mult.commute) 
    subgoal by (metis mat2.sel(2)) 
    subgoal by (metis mat2.sel(3))
    subgoal using mat2.sel(4)[of 1 0 0 1] apply (auto) 
          using A_inv_aux[of b n] by (auto simp: mult.commute) 
  done

lemma B_inverse[simp]: "mat_mul (B b) (B_inv b) = ID" using B.simps ID_def by auto

declare A_inv.simps B_inv.simps[simp del]


text ‹Equation 3.33›
lemma congruence:
  assumes "b1 mod q = b2 mod q"
  shows "α b1 n mod q = α b2 n mod q"
proof (induct n rule:nat_less_induct)
  case (1 n)
  note hyps = ‹∀m<n. α b1 m mod q = α b2 m mod q›
  have n0:"(α b1 0) mod q = (α b2 0) mod q" by simp
  have n1:"(α b1 1) mod q = (α b2 1) mod q" by simp
  from hyps have s1: "n>1⟹α b1 (n-1) mod q = α b2 (n-1) mod q" by auto
  from hyps have s2: "n>1⟹α b1 (n-2) mod q = α b2 (n-2) mod q" by auto
  have s3: "n>1 ⟹ α b1 (Suc (Suc n)) = (int b1) * (α b1 (Suc n)) - (α b1 n)" by simp
  from s3 have s4: "n>1 ⟹ (α b1 n = (int b1*(α b1 (n-1)) - α b1 (n-2)))" 
    by (smt Suc_1 Suc_diff_Suc diff_Suc_1 alpha_n lessE)
  have sw: "n>1 ⟹ α b2 (Suc (Suc n)) = (int b2) * (α b2 (Suc n)) - (α b2 n)" by simp
  from sw have sx: "n>1 ⟹ (α b2 n = (int b2*(α b2 (n-1)) - α b2 (n-2)))" 
    by (smt Suc_1 Suc_diff_Suc diff_Suc_1 alpha_n lessE)
  from n0 n1 s1 s2 s3 s4 assms(1) mod_mult_cong have s5: "n>1 
        ⟹ b1*(α b1 (n-1)) mod q =  b2*(α b2 (n-1)) mod q " by (smt mod_mult_eq of_nat_mod)
   from hyps have sq: "n>1 ⟹ α b1 (n-2) mod q = α b2 (n-2) mod q " by simp
   from s5 sq have sd: "n>1 ⟹-( (α b1 (n-2))) mod q = -((α b2 (n-2))) mod q " 
     by (metis mod_minus_eq)
   from sd s5 mod_add_cong have s6: "n>1 ⟹ ( b1*(α b1 (n-1)) - α b1 (n-2)) mod q  
                                      =  (b2*(α b2 (n-1)) - α b2 (n-2)) mod q" by force
  from s4 have sa: "n>1 ⟹( b1*(α b1 (n-1)) - α b1 (n-2)) mod q = (α b1 n) mod q" by simp
  from sx have sb: "n>1 ⟹( b2*(α b2 (n-1)) - α b2 (n-2)) mod q = (α b2 n) mod q" by simp
  from sb sa s6 sx have s7: "n>1 ⟹ (α b1 n) mod q = (b2*(α b2 (n-1)) - α b2 (n-2)) mod q" by simp
  from s7 sx s6 have s9: "α b1 n mod q = α b2 n mod q" 
    by (metis One_nat_def α.simps(1) α.simps(2) less_Suc0 nat_neq_iff)
  from s9 n0 n1 show ?case by simp
qed

text ‹Equation 3.34›
lemma congruence2:
  fixes b1 :: nat
  assumes "b>=2"
  shows "(α b n) mod (b - 2) = n mod (b - 2)"
proof-
  from alpha_linear have S1: "α (nat 2) n = n" by simp
  define q where "q = b - (nat 2)"
  from q_def assms le_mod_geq have S4: "b mod q = 2 mod q" by auto
  from assms S4 congruence have SN: "(α b n) mod q = (α 2 n) mod q" by blast
  from S1 SN q_def zmod_int show ?thesis by simp
qed

lemma congruence_jpos:
  fixes b m j l :: nat
  assumes "b>2" and "2*l*m+j>0"
  defines "n ≡ 2*l*m+j"
  shows "A b n = mat_mul (mat_pow l (mat_pow 2 (A b m))) (A b j)"
proof-
  from A_pow assms(1) have Abm2: "A b n = mat_pow n (B b)" by simp
  from Abm2 n_def have Bn: "mat_pow n (B b) =mat_pow (2*l*m+j) (B b)" by simp
  from mat_exp_law have as1: "mat_pow (2*l*m+j) (