--- a/CONTRIBUTORS Mon Feb 04 16:01:44 2019 +0100
+++ b/CONTRIBUTORS Mon Feb 04 15:39:37 2019 +0100
@@ -7,6 +7,10 @@
--------------------------------------
* February 2019: Manuel Eberl
+ Exponentiation by squaring, used to implement "power" in monoid_mult and
+ fast modular exponentiation.
+
+* February 2019: Manuel Eberl
Carmichael's function, primitive roots in residue rings, more properties
of the order in residue rings.
--- a/NEWS Mon Feb 04 16:01:44 2019 +0100
+++ b/NEWS Mon Feb 04 15:39:37 2019 +0100
@@ -87,6 +87,8 @@
*** HOL ***
+* exponentiation by squaring in HOL-Library; used for computing powers in monoid_mult and modular exponentiation in HOL-Number_Theory
+
* more material on residue rings in HOL-Number_Theory:
Carmichael's function, primitive roots, more properties for "ord"
--- a/src/HOL/Library/Library.thy Mon Feb 04 16:01:44 2019 +0100
+++ b/src/HOL/Library/Library.thy Mon Feb 04 15:39:37 2019 +0100
@@ -63,6 +63,7 @@
Perm
Permutation
Permutations
+ Power_By_Squaring
Preorder
Product_Plus
Quadratic_Discriminant
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/Power_By_Squaring.thy Mon Feb 04 15:39:37 2019 +0100
@@ -0,0 +1,68 @@
+(*
+ File: Power_By_Squaring.thy
+ Author: Manuel Eberl, TU München
+
+ Fast computing of funpow (applying some functon n times) for weakly associative binary
+ functions using exponentiation by squaring. Yields efficient exponentiation algorithms on
+ monoid_mult and for modular exponentiation "b ^ e mod m" (and thus also for "cong")
+*)
+section \<open>Exponentiation by Squaring\<close>
+theory Power_By_Squaring
+ imports Main
+begin
+
+context
+ fixes f :: "'a \<Rightarrow> 'a \<Rightarrow> 'a"
+begin
+
+function efficient_funpow :: "'a \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> 'a" where
+ "efficient_funpow y x 0 = y"
+| "efficient_funpow y x (Suc 0) = f x y"
+| "n \<noteq> 0 \<Longrightarrow> even n \<Longrightarrow> efficient_funpow y x n = efficient_funpow y (f x x) (n div 2)"
+| "n \<noteq> 1 \<Longrightarrow> odd n \<Longrightarrow> efficient_funpow y x n = efficient_funpow (f x y) (f x x) (n div 2)"
+ by force+
+termination by (relation "measure (snd \<circ> snd)") (auto elim: oddE)
+
+lemma efficient_funpow_code [code]:
+ "efficient_funpow y x n =
+ (if n = 0 then y
+ else if n = 1 then f x y
+ else if even n then efficient_funpow y (f x x) (n div 2)
+ else efficient_funpow (f x y) (f x x) (n div 2))"
+ by (induction y x n rule: efficient_funpow.induct) auto
+
+end
+
+lemma efficient_funpow_correct:
+ assumes f_assoc: "\<And>x z. f x (f x z) = f (f x x) z"
+ shows "efficient_funpow f y x n = (f x ^^ n) y"
+proof -
+ have [simp]: "f ^^ 2 = (\<lambda>x. f (f x))" for f :: "'a \<Rightarrow> 'a"
+ by (simp add: eval_nat_numeral o_def)
+ show ?thesis
+ by (induction y x n rule: efficient_funpow.induct[of _ f])
+ (auto elim!: evenE oddE simp: funpow_mult [symmetric] funpow_Suc_right f_assoc
+ simp del: funpow.simps(2))
+qed
+
+(*
+ TODO: This could be used as a code_unfold rule or something like that but the
+ implications are not quite clear. Would this be a good default implementation
+ for powers?
+*)
+context monoid_mult
+begin
+
+lemma power_by_squaring: "efficient_funpow (*) (1 :: 'a) = (^)"
+proof (intro ext)
+ fix x :: 'a and n
+ have "efficient_funpow (*) 1 x n = ((*) x ^^ n) 1"
+ by (subst efficient_funpow_correct) (simp_all add: mult.assoc)
+ also have "\<dots> = x ^ n"
+ by (induction n) simp_all
+ finally show "efficient_funpow (*) 1 x n = x ^ n" .
+qed
+
+end
+
+end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Number_Theory/Mod_Exp.thy Mon Feb 04 15:39:37 2019 +0100
@@ -0,0 +1,117 @@
+(*
+ File: HOL/Number_Theory/Mod_Exp
+ Author: Manuel Eberl, TU München
+
+ Fast implementation of modular exponentiation and "cong" using exponentiation by squaring.
+ Includes code setup for nat and int.
+*)
+section \<open>Fast modular exponentiation\<close>
+theory Mod_Exp
+ imports Cong "HOL-Library.Power_By_Squaring"
+begin
+
+context euclidean_semiring_cancel
+begin
+
+definition mod_exp_aux :: "'a \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> 'a"
+ where "mod_exp_aux m = efficient_funpow (\<lambda>x y. x * y mod m)"
+
+lemma mod_exp_aux_code [code]:
+ "mod_exp_aux m y x n =
+ (if n = 0 then y
+ else if n = 1 then (x * y) mod m
+ else if even n then mod_exp_aux m y ((x * x) mod m) (n div 2)
+ else mod_exp_aux m ((x * y) mod m) ((x * x) mod m) (n div 2))"
+ unfolding mod_exp_aux_def by (rule efficient_funpow_code)
+
+lemma mod_exp_aux_correct:
+ "mod_exp_aux m y x n mod m = (x ^ n * y) mod m"
+proof -
+ have "mod_exp_aux m y x n = efficient_funpow (\<lambda>x y. x * y mod m) y x n"
+ by (simp add: mod_exp_aux_def)
+ also have "\<dots> = ((\<lambda>y. x * y mod m) ^^ n) y"
+ by (rule efficient_funpow_correct) (simp add: mod_mult_left_eq mod_mult_right_eq mult_ac)
+ also have "((\<lambda>y. x * y mod m) ^^ n) y mod m = (x ^ n * y) mod m"
+ proof (induction n)
+ case (Suc n)
+ hence "x * ((\<lambda>y. x * y mod m) ^^ n) y mod m = x * x ^ n * y mod m"
+ by (metis mod_mult_right_eq mult.assoc)
+ thus ?case by auto
+ qed auto
+ finally show ?thesis .
+qed
+
+definition mod_exp :: "'a \<Rightarrow> nat \<Rightarrow> 'a \<Rightarrow> 'a"
+ where "mod_exp b e m = (b ^ e) mod m"
+
+lemma mod_exp_code [code]: "mod_exp b e m = mod_exp_aux m 1 b e mod m"
+ by (simp add: mod_exp_def mod_exp_aux_correct)
+
+end
+
+(*
+ TODO: Setup here only for nat and int. Could be done for any
+ euclidean_semiring_cancel. Should it?
+*)
+lemmas [code_abbrev] = mod_exp_def[where ?'a = nat] mod_exp_def[where ?'a = int]
+
+lemma cong_power_nat_code [code_unfold]:
+ "[b ^ e = (x ::nat)] (mod m) \<longleftrightarrow> mod_exp b e m = x mod m"
+ by (simp add: mod_exp_def cong_def)
+
+lemma cong_power_int_code [code_unfold]:
+ "[b ^ e = (x ::int)] (mod m) \<longleftrightarrow> mod_exp b e m = x mod m"
+ by (simp add: mod_exp_def cong_def)
+
+
+text \<open>
+ The following rules allow the simplifier to evaluate @{const mod_exp} efficiently.
+\<close>
+lemma eval_mod_exp_aux [simp]:
+ "mod_exp_aux m y x 0 = y"
+ "mod_exp_aux m y x (Suc 0) = (x * y) mod m"
+ "mod_exp_aux m y x (numeral (num.Bit0 n)) =
+ mod_exp_aux m y (x\<^sup>2 mod m) (numeral n)"
+ "mod_exp_aux m y x (numeral (num.Bit1 n)) =
+ mod_exp_aux m ((x * y) mod m) (x\<^sup>2 mod m) (numeral n)"
+proof -
+ define n' where "n' = (numeral n :: nat)"
+ have [simp]: "n' \<noteq> 0" by (auto simp: n'_def)
+
+ show "mod_exp_aux m y x 0 = y" and "mod_exp_aux m y x (Suc 0) = (x * y) mod m"
+ by (simp_all add: mod_exp_aux_def)
+
+ have "numeral (num.Bit0 n) = (2 * n')"
+ by (subst numeral.numeral_Bit0) (simp del: arith_simps add: n'_def)
+ also have "mod_exp_aux m y x \<dots> = mod_exp_aux m y (x^2 mod m) n'"
+ by (subst mod_exp_aux_code) (simp_all add: power2_eq_square)
+ finally show "mod_exp_aux m y x (numeral (num.Bit0 n)) =
+ mod_exp_aux m y (x\<^sup>2 mod m) (numeral n)"
+ by (simp add: n'_def)
+
+ have "numeral (num.Bit1 n) = Suc (2 * n')"
+ by (subst numeral.numeral_Bit1) (simp del: arith_simps add: n'_def)
+ also have "mod_exp_aux m y x \<dots> = mod_exp_aux m ((x * y) mod m) (x^2 mod m) n'"
+ by (subst mod_exp_aux_code) (simp_all add: power2_eq_square)
+ finally show "mod_exp_aux m y x (numeral (num.Bit1 n)) =
+ mod_exp_aux m ((x * y) mod m) (x\<^sup>2 mod m) (numeral n)"
+ by (simp add: n'_def)
+qed
+
+lemma eval_mod_exp [simp]:
+ "mod_exp b' 0 m' = 1 mod m'"
+ "mod_exp b' 1 m' = b' mod m'"
+ "mod_exp b' (Suc 0) m' = b' mod m'"
+ "mod_exp b' e' 0 = b' ^ e'"
+ "mod_exp b' e' 1 = 0"
+ "mod_exp b' e' (Suc 0) = 0"
+ "mod_exp 0 1 m' = 0"
+ "mod_exp 0 (Suc 0) m' = 0"
+ "mod_exp 0 (numeral e) m' = 0"
+ "mod_exp 1 e' m' = 1 mod m'"
+ "mod_exp (Suc 0) e' m' = 1 mod m'"
+ "mod_exp (numeral b) (numeral e) (numeral m) =
+ mod_exp_aux (numeral m) 1 (numeral b) (numeral e) mod numeral m"
+ by (simp_all add: mod_exp_def mod_exp_aux_correct)
+
+end
--- a/src/HOL/Number_Theory/Number_Theory.thy Mon Feb 04 16:01:44 2019 +0100
+++ b/src/HOL/Number_Theory/Number_Theory.thy Mon Feb 04 15:39:37 2019 +0100
@@ -6,6 +6,7 @@
Fib
Residues
Eratosthenes
+ Mod_Exp
Quadratic_Reciprocity
Pocklington
Prime_Powers