Theory Log_Impl
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