Exponentiation by squaring, fast modular exponentiation
authorManuel Eberl <eberlm@in.tum.de>
Mon Feb 04 15:39:37 2019 +0100 (4 months ago)
changeset 69790154cf64e403e
parent 69789 2c3e5e58d93f
child 69791 195aeee8b30a
Exponentiation by squaring, fast modular exponentiation
CONTRIBUTORS
NEWS
src/HOL/Library/Library.thy
src/HOL/Library/Power_By_Squaring.thy
src/HOL/Number_Theory/Mod_Exp.thy
src/HOL/Number_Theory/Number_Theory.thy
     1.1 --- a/CONTRIBUTORS	Mon Feb 04 16:01:44 2019 +0100
     1.2 +++ b/CONTRIBUTORS	Mon Feb 04 15:39:37 2019 +0100
     1.3 @@ -7,6 +7,10 @@
     1.4  --------------------------------------
     1.5  
     1.6  * February 2019: Manuel Eberl
     1.7 +  Exponentiation by squaring, used to implement "power" in monoid_mult and
     1.8 +  fast modular exponentiation.
     1.9 +
    1.10 +* February 2019: Manuel Eberl
    1.11    Carmichael's function, primitive roots in residue rings, more properties
    1.12    of the order in residue rings.
    1.13  
     2.1 --- a/NEWS	Mon Feb 04 16:01:44 2019 +0100
     2.2 +++ b/NEWS	Mon Feb 04 15:39:37 2019 +0100
     2.3 @@ -87,6 +87,8 @@
     2.4  
     2.5  *** HOL ***
     2.6  
     2.7 +* exponentiation by squaring in HOL-Library; used for computing powers in monoid_mult and modular exponentiation in HOL-Number_Theory
     2.8 +
     2.9  * more material on residue rings in HOL-Number_Theory:
    2.10  Carmichael's function, primitive roots, more properties for "ord"
    2.11  
     3.1 --- a/src/HOL/Library/Library.thy	Mon Feb 04 16:01:44 2019 +0100
     3.2 +++ b/src/HOL/Library/Library.thy	Mon Feb 04 15:39:37 2019 +0100
     3.3 @@ -63,6 +63,7 @@
     3.4    Perm
     3.5    Permutation
     3.6    Permutations
     3.7 +  Power_By_Squaring
     3.8    Preorder
     3.9    Product_Plus
    3.10    Quadratic_Discriminant
     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/src/HOL/Library/Power_By_Squaring.thy	Mon Feb 04 15:39:37 2019 +0100
     4.3 @@ -0,0 +1,68 @@
     4.4 +(*
     4.5 +  File:     Power_By_Squaring.thy
     4.6 +  Author:   Manuel Eberl, TU München
     4.7 +  
     4.8 +  Fast computing of funpow (applying some functon n times) for weakly associative binary
     4.9 +  functions using exponentiation by squaring. Yields efficient exponentiation algorithms on
    4.10 +  monoid_mult and for modular exponentiation "b ^ e mod m" (and thus also for "cong")
    4.11 +*)
    4.12 +section \<open>Exponentiation by Squaring\<close>
    4.13 +theory Power_By_Squaring
    4.14 +  imports Main
    4.15 +begin
    4.16 +
    4.17 +context
    4.18 +  fixes f :: "'a \<Rightarrow> 'a \<Rightarrow> 'a"
    4.19 +begin
    4.20 +
    4.21 +function efficient_funpow :: "'a \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> 'a" where
    4.22 +  "efficient_funpow y x 0 = y"
    4.23 +| "efficient_funpow y x (Suc 0) = f x y"
    4.24 +| "n \<noteq> 0 \<Longrightarrow> even n \<Longrightarrow> efficient_funpow y x n = efficient_funpow y (f x x) (n div 2)"
    4.25 +| "n \<noteq> 1 \<Longrightarrow> odd n \<Longrightarrow> efficient_funpow y x n = efficient_funpow (f x y) (f x x) (n div 2)"
    4.26 +  by force+
    4.27 +termination by (relation "measure (snd \<circ> snd)") (auto elim: oddE)
    4.28 +
    4.29 +lemma efficient_funpow_code [code]:
    4.30 +  "efficient_funpow y x n =
    4.31 +     (if n = 0 then y
    4.32 +      else if n = 1 then f x y
    4.33 +      else if even n then efficient_funpow y (f x x) (n div 2)
    4.34 +      else efficient_funpow (f x y) (f x x) (n div 2))"
    4.35 +  by (induction y x n rule: efficient_funpow.induct) auto
    4.36 +
    4.37 +end
    4.38 +
    4.39 +lemma efficient_funpow_correct:
    4.40 +  assumes f_assoc: "\<And>x z. f x (f x z) = f (f x x) z"
    4.41 +  shows "efficient_funpow f y x n = (f x ^^ n) y"
    4.42 +proof -
    4.43 +  have [simp]: "f ^^ 2 = (\<lambda>x. f (f x))" for f :: "'a \<Rightarrow> 'a"
    4.44 +    by (simp add: eval_nat_numeral o_def)
    4.45 +  show ?thesis
    4.46 +    by (induction y x n rule: efficient_funpow.induct[of _ f])
    4.47 +       (auto elim!: evenE oddE simp: funpow_mult [symmetric] funpow_Suc_right f_assoc
    4.48 +             simp del: funpow.simps(2))
    4.49 +qed
    4.50 +
    4.51 +(*
    4.52 +  TODO: This could be used as a code_unfold rule or something like that but the
    4.53 +  implications are not quite clear. Would this be a good default implementation
    4.54 +  for powers?
    4.55 +*)
    4.56 +context monoid_mult
    4.57 +begin
    4.58 +
    4.59 +lemma power_by_squaring: "efficient_funpow (*) (1 :: 'a) = (^)"
    4.60 +proof (intro ext)
    4.61 +  fix x :: 'a and n
    4.62 +  have "efficient_funpow (*) 1 x n = ((*) x ^^ n) 1"
    4.63 +    by (subst efficient_funpow_correct) (simp_all add: mult.assoc)
    4.64 +  also have "\<dots> = x ^ n"
    4.65 +    by (induction n) simp_all
    4.66 +  finally show "efficient_funpow (*) 1 x n = x ^ n" .
    4.67 +qed
    4.68 +
    4.69 +end
    4.70 +
    4.71 +end
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/src/HOL/Number_Theory/Mod_Exp.thy	Mon Feb 04 15:39:37 2019 +0100
     5.3 @@ -0,0 +1,117 @@
     5.4 +(*
     5.5 +  File:    HOL/Number_Theory/Mod_Exp
     5.6 +  Author:  Manuel Eberl, TU München
     5.7 +
     5.8 +  Fast implementation of modular exponentiation and "cong" using exponentiation by squaring.
     5.9 +  Includes code setup for nat and int.
    5.10 +*)
    5.11 +section \<open>Fast modular exponentiation\<close>
    5.12 +theory Mod_Exp
    5.13 +  imports Cong "HOL-Library.Power_By_Squaring"
    5.14 +begin
    5.15 +
    5.16 +context euclidean_semiring_cancel
    5.17 +begin
    5.18 +
    5.19 +definition mod_exp_aux :: "'a \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> 'a"
    5.20 +  where "mod_exp_aux m = efficient_funpow (\<lambda>x y. x * y mod m)"
    5.21 +
    5.22 +lemma mod_exp_aux_code [code]:
    5.23 +  "mod_exp_aux m y x n =
    5.24 +     (if n = 0 then y
    5.25 +      else if n = 1 then (x * y) mod m
    5.26 +      else if even n then mod_exp_aux m y ((x * x) mod m) (n div 2)
    5.27 +      else mod_exp_aux m ((x * y) mod m) ((x * x) mod m) (n div 2))"
    5.28 +  unfolding mod_exp_aux_def by (rule efficient_funpow_code)
    5.29 +
    5.30 +lemma mod_exp_aux_correct:
    5.31 +  "mod_exp_aux m y x n mod m = (x ^ n * y) mod m"
    5.32 +proof -
    5.33 +  have "mod_exp_aux m y x n = efficient_funpow (\<lambda>x y. x * y mod m) y x n"
    5.34 +    by (simp add: mod_exp_aux_def)
    5.35 +  also have "\<dots> = ((\<lambda>y. x * y mod m) ^^ n) y"
    5.36 +    by (rule efficient_funpow_correct) (simp add: mod_mult_left_eq mod_mult_right_eq mult_ac)
    5.37 +  also have "((\<lambda>y. x * y mod m) ^^ n) y mod m = (x ^ n * y) mod m"
    5.38 +  proof (induction n)
    5.39 +    case (Suc n)
    5.40 +    hence "x * ((\<lambda>y. x * y mod m) ^^ n) y mod m = x * x ^ n * y mod m"
    5.41 +      by (metis mod_mult_right_eq mult.assoc)
    5.42 +    thus ?case by auto
    5.43 +  qed auto
    5.44 +  finally show ?thesis .
    5.45 +qed
    5.46 +
    5.47 +definition mod_exp :: "'a \<Rightarrow> nat \<Rightarrow> 'a \<Rightarrow> 'a"
    5.48 +  where "mod_exp b e m = (b ^ e) mod m"
    5.49 +
    5.50 +lemma mod_exp_code [code]: "mod_exp b e m = mod_exp_aux m 1 b e mod m"
    5.51 +  by (simp add: mod_exp_def mod_exp_aux_correct)
    5.52 +
    5.53 +end
    5.54 +
    5.55 +(*
    5.56 +  TODO: Setup here only for nat and int. Could be done for any
    5.57 +  euclidean_semiring_cancel. Should it?
    5.58 +*)
    5.59 +lemmas [code_abbrev] = mod_exp_def[where ?'a = nat] mod_exp_def[where ?'a = int]
    5.60 +
    5.61 +lemma cong_power_nat_code [code_unfold]:
    5.62 +  "[b ^ e = (x ::nat)] (mod m) \<longleftrightarrow> mod_exp b e m = x mod m"
    5.63 +  by (simp add: mod_exp_def cong_def)
    5.64 +
    5.65 +lemma cong_power_int_code [code_unfold]:
    5.66 +  "[b ^ e = (x ::int)] (mod m) \<longleftrightarrow> mod_exp b e m = x mod m"
    5.67 +  by (simp add: mod_exp_def cong_def)
    5.68 +
    5.69 +
    5.70 +text \<open>
    5.71 +  The following rules allow the simplifier to evaluate @{const mod_exp} efficiently.
    5.72 +\<close>
    5.73 +lemma eval_mod_exp_aux [simp]:
    5.74 +  "mod_exp_aux m y x 0 = y"
    5.75 +  "mod_exp_aux m y x (Suc 0) = (x * y) mod m"
    5.76 +  "mod_exp_aux m y x (numeral (num.Bit0 n)) =
    5.77 +     mod_exp_aux m y (x\<^sup>2 mod m) (numeral n)"
    5.78 +  "mod_exp_aux m y x (numeral (num.Bit1 n)) =
    5.79 +     mod_exp_aux m ((x * y) mod m) (x\<^sup>2 mod m) (numeral n)"
    5.80 +proof -
    5.81 +  define n' where "n' = (numeral n :: nat)"
    5.82 +  have [simp]: "n' \<noteq> 0" by (auto simp: n'_def)
    5.83 +  
    5.84 +  show "mod_exp_aux m y x 0 = y" and "mod_exp_aux m y x (Suc 0) = (x * y) mod m"
    5.85 +    by (simp_all add: mod_exp_aux_def)
    5.86 +
    5.87 +  have "numeral (num.Bit0 n) = (2 * n')"
    5.88 +    by (subst numeral.numeral_Bit0) (simp del: arith_simps add: n'_def)
    5.89 +  also have "mod_exp_aux m y x \<dots> = mod_exp_aux m y (x^2 mod m) n'"
    5.90 +    by (subst mod_exp_aux_code) (simp_all add: power2_eq_square)
    5.91 +  finally show "mod_exp_aux m y x (numeral (num.Bit0 n)) =
    5.92 +                  mod_exp_aux m y (x\<^sup>2 mod m) (numeral n)"
    5.93 +    by (simp add: n'_def)
    5.94 +
    5.95 +  have "numeral (num.Bit1 n) = Suc (2 * n')"
    5.96 +    by (subst numeral.numeral_Bit1) (simp del: arith_simps add: n'_def)
    5.97 +  also have "mod_exp_aux m y x \<dots> = mod_exp_aux m ((x * y) mod m) (x^2 mod m) n'"
    5.98 +    by (subst mod_exp_aux_code) (simp_all add: power2_eq_square)
    5.99 +  finally show "mod_exp_aux m y x (numeral (num.Bit1 n)) =
   5.100 +                  mod_exp_aux m ((x * y) mod m) (x\<^sup>2 mod m) (numeral n)"
   5.101 +    by (simp add: n'_def)
   5.102 +qed
   5.103 +
   5.104 +lemma eval_mod_exp [simp]:
   5.105 +  "mod_exp b' 0 m' = 1 mod m'"
   5.106 +  "mod_exp b' 1 m' = b' mod m'"
   5.107 +  "mod_exp b' (Suc 0) m' = b' mod m'"
   5.108 +  "mod_exp b' e' 0 = b' ^ e'"  
   5.109 +  "mod_exp b' e' 1 = 0"
   5.110 +  "mod_exp b' e' (Suc 0) = 0"
   5.111 +  "mod_exp 0 1 m' = 0"
   5.112 +  "mod_exp 0 (Suc 0) m' = 0"
   5.113 +  "mod_exp 0 (numeral e) m' = 0"
   5.114 +  "mod_exp 1 e' m' = 1 mod m'"
   5.115 +  "mod_exp (Suc 0) e' m' = 1 mod m'"
   5.116 +  "mod_exp (numeral b) (numeral e) (numeral m) =
   5.117 +     mod_exp_aux (numeral m) 1 (numeral b) (numeral e) mod numeral m"
   5.118 +  by (simp_all add: mod_exp_def mod_exp_aux_correct)
   5.119 +
   5.120 +end
     6.1 --- a/src/HOL/Number_Theory/Number_Theory.thy	Mon Feb 04 16:01:44 2019 +0100
     6.2 +++ b/src/HOL/Number_Theory/Number_Theory.thy	Mon Feb 04 15:39:37 2019 +0100
     6.3 @@ -6,6 +6,7 @@
     6.4    Fib
     6.5    Residues
     6.6    Eratosthenes
     6.7 +  Mod_Exp
     6.8    Quadratic_Reciprocity
     6.9    Pocklington
    6.10    Prime_Powers