Theory Sqrt_Babylonian.Log_Impl

(*  Title:       Computing Square Roots using the Babylonian Method
    Author:      René Thiemann       <rene.thiemann@uibk.ac.at>
    Maintainer:  René Thiemann
    License:     LGPL
*)
section ‹A Fast Logarithm Algorithm›

theory Log_Impl
imports 
  Sqrt_Babylonian_Auxiliary
begin

text ‹We implement the discrete logarithm function in a manner similar to
  a repeated squaring exponentiation algorithm.›

text ‹In order to prove termination of the algorithm without intermediate checks 
  we need to ensure that we only use proper bases, 
  i.e., values of at least 2. This will be encoded into a separate type.›

typedef proper_base = "{x :: int. x ≥ 2}" by auto

setup_lifting type_definition_proper_base

lift_definition get_base :: "proper_base ⇒ int" is "λ x. x" .

lift_definition square_base :: "proper_base ⇒ proper_base" is "λ x. x * x" 
proof -
  fix i :: int
  assume i: "2 ≤ i"
  have "2 * 2 ≤ i * i" 
    by (rule mult_mono[OF i i], insert i, auto)
  thus "2 ≤ i * i" by auto
qed

lift_definition into_base :: "int ⇒ proper_base" is "λ x. if x ≥ 2 then x else 2" by auto

lemma square_base: "get_base (square_base b) = get_base b * get_base b" 
  by (transfer, auto)

lemma get_base_2: "get_base b ≥ 2"
  by (transfer, auto)

lemma b_less_square_base_b: "get_base b < get_base (square_base b)" 
  unfolding square_base using get_base_2[of b] by simp

lemma b_less_div_base_b: assumes xb: "¬ x < get_base b"
  shows "x div get_base b < x"
proof -
  from get_base_2[of b] have b: "get_base b ≥ 2" .
  with xb have x2: "x ≥ 2" by auto
  with b int_div_less_self[of x "(get_base b)"] 
  show ?thesis by auto
qed
    
text ‹We now state the main algorithm.›
    
function log_main :: "proper_base ⇒ int ⇒ nat × int" where
  "log_main b x = (if x < get_base b then (0,1) else
    case log_main (square_base b) x of 
      (z, bz) ⇒ 
    let l = 2 * z; bz1 = bz * get_base b
      in if x < bz1 then (l,bz) else (Suc l,bz1))" 
  by pat_completeness auto

termination by (relation "measure (λ (b,x). nat (1 + x - get_base b))",
  insert b_less_square_base_b, auto)   

lemma log_main: "x > 0 ⟹ log_main b x = (y,by) ⟹ by = (get_base b)^y ∧ (get_base b)^y ≤ x ∧ x < (get_base b)^(Suc y)" 
proof (induct b x arbitrary: y "by" rule: log_main.induct)
  case (1 b x y "by")
  note x = 1(2)
  note y = 1(3)
  note IH = 1(1)
  let ?b = "get_base b" 
  show ?case
  proof (cases "x < ?b")
    case True
    with x y show ?thesis by auto
  next
    case False
    obtain z bz where zz: "log_main (square_base b) x = (z,bz)" 
      by (cases "log_main (square_base b) x", auto)
    have id: "get_base (square_base b) ^ k = ?b ^ (2 * k)" for k unfolding square_base
      by (simp add: power_mult semiring_normalization_rules(29))
    from IH[OF False x zz, unfolded id] 
    have z: "?b ^ (2 * z) ≤ x" "x < ?b ^ (2 * Suc z)" and bz: "bz = get_base b ^ (2 * z)" by auto
    from y[unfolded log_main.simps[of b x] Let_def zz split] bz False
    have yy: "(if x < bz * ?b then (2 * z, bz) else (Suc (2 * z), bz * ?b)) =
      (y, by)" by auto
    show ?thesis
    proof (cases "x < bz * ?b")
      case True
      with yy have yz: "y = 2 * z" "by = bz" by auto
      from True z(1) bz show ?thesis unfolding yz by (auto simp: ac_simps)
    next
      case False
      with yy have yz: "y = Suc (2 * z)" "by = ?b * bz" by auto
      from False have "?b ^ Suc (2 * z) ≤ x" by (auto simp: bz ac_simps)
      with z(2) bz show ?thesis unfolding yz by auto
    qed
  qed
qed
    
text ‹We then derive the floor- and ceiling-log functions.›

definition log_floor :: "int ⇒ int ⇒ nat" where
  "log_floor b x = fst (log_main (into_base b) x)" 

definition log_ceiling :: "int ⇒ int ⇒ nat" where
  "log_ceiling b x = (case log_main (into_base b) x of
     (y,by) ⇒ if x = by then y else Suc y)" 

lemma log_floor_sound: assumes "b > 1" "x > 0" "log_floor b x = y"  
  shows "b^y ≤ x" "x < b^(Suc y)" 
proof -
  from assms(1,3) have id: "get_base (into_base b) = b" by transfer auto
  obtain yy bb where log: "log_main (into_base b) x = (yy,bb)" 
    by (cases "log_main (into_base b) x", auto)
  from log_main[OF assms(2) log] assms(3)[unfolded log_floor_def log] id
  show "b^y ≤ x" "x < b^(Suc y)" by auto
qed

lemma log_ceiling_sound: assumes "b > 1" "x > 0" "log_ceiling b x = y"  
  shows "x ≤ b^y" "y ≠ 0 ⟹ b^(y - 1) < x" 
proof -
  from assms(1,3) have id: "get_base (into_base b) = b" by transfer auto
  obtain yy bb where log: "log_main (into_base b) x = (yy,bb)" 
    by (cases "log_main (into_base b) x", auto)
  from log_main[OF assms(2) log, unfolded id] assms(3)[unfolded log_ceiling_def log split]
  have bnd: "b ^ yy ≤ x" "x < b ^ Suc yy" and
    y: "y = (if x = b ^ yy then yy else Suc yy)" by auto
  have "x ≤ b^y ∧ (y ≠ 0 ⟶ b^(y - 1) < x)"
  proof (cases "x = b ^ yy")
    case True
    with y bnd assms(1) show ?thesis by (cases yy, auto)
  next
    case False
    with y bnd show ?thesis by auto
  qed
  thus "x ≤ b^y" "y ≠ 0 ⟹ b^(y - 1) < x" by auto
qed

text ‹Finally, we connect it to the @{const log} function working on real numbers.›

lemma log_floor[simp]: assumes b: "b > 1" and x: "x > 0"
  shows "log_floor b x = ⌊log b x⌋"
proof -
  obtain y where y: "log_floor b x = y" by auto
  note main = log_floor_sound[OF assms y]
  from b x have *: "1 < real_of_int b" "0 < real_of_int (b ^ y)" "0 < real_of_int x" 
    and **: "1 < real_of_int b" "0 < real_of_int x" "0 < real_of_int (b ^ Suc y)" 
    by auto
  show ?thesis unfolding y
  proof (rule sym, rule floor_unique)
    show "real_of_int (int y) ≤ log (real_of_int b) (real_of_int x)" 
      using main(1)[folded log_le_cancel_iff[OF *, unfolded of_int_le_iff]]
      using log_pow_cancel[of b y] b by auto
    show "log (real_of_int b) (real_of_int x) < real_of_int (int y) + 1" 
      using main(2)[folded log_less_cancel_iff[OF **, unfolded of_int_less_iff]]
      using log_pow_cancel[of b "Suc y"] b by auto
  qed
qed
    
lemma log_ceiling[simp]: assumes b: "b > 1" and x: "x > 0"
  shows "log_ceiling b x = ⌈log b x⌉"
proof -
  obtain y where y: "log_ceiling b x = y" by auto
  note main = log_ceiling_sound[OF assms y]
  from b x have *: "1 < real_of_int b" "0 < real_of_int (b ^ (y - 1))" "0 < real_of_int x" 
    and **: "1 < real_of_int b" "0 < real_of_int x" "0 < real_of_int (b ^ y)" 
    by auto
  show ?thesis unfolding y
  proof (rule sym, rule ceiling_unique)
    show "log (real_of_int b) (real_of_int x) ≤ real_of_int (int y)" 
      using main(1)[folded log_le_cancel_iff[OF **, unfolded of_int_le_iff]]
      using log_pow_cancel[of b y] b by auto
    from x have x: "x ≥ 1" by auto
    show "real_of_int (int y) - 1 < log (real_of_int b) (real_of_int x)" 
    proof (cases "y = 0")
      case False
      thus ?thesis 
        using main(2)[folded log_less_cancel_iff[OF *, unfolded of_int_less_iff]]
        using log_pow_cancel[of b "y - 1"] b x by auto
    next
      case True
      have "real_of_int (int y) - 1 = log b (1/b)" using True b 
        by (subst log_divide, auto)
      also have "… < log b 1"
        by (subst log_less_cancel_iff, insert b, auto) 
      also have "… ≤ log b x" 
        by (subst log_le_cancel_iff, insert b x, auto) 
      finally show "real_of_int (int y) - 1 < log (real_of_int b) (real_of_int x)" .
    qed
  qed
qed

end