1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Tue May 21 11:33:03 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Thu May 23 12:20:28 2013 +0200
1.3 @@ -114,10 +114,10 @@
1.4 else modi (rold * (rinv + 1), rold, m, rinv + 1)
1.5 in modi (r, r, m, 1) end
1.6
1.7 - (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat
1.8 + (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
1.9 mod_div a b m = c
1.10 assumes True
1.11 - yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = *)
1.12 + yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
1.13 fun mod_div a b m = a * (mod_inv b m) mod m
1.14
1.15 (* required just for one approximation:
1.16 @@ -356,21 +356,21 @@
1.17 fun try_new_prime_up a b d M P g p =
1.18 if P > M then g else
1.19 let
1.20 - val p' = p next_prime_not_dvd d
1.21 - val g_p = centr_up (( (normalise (mod_up_gcd a b p') p')
1.22 - %* (d mod p'))
1.23 - mod_up p')
1.24 - p'
1.25 + val p = p next_prime_not_dvd d
1.26 + val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
1.27 + %* (d mod p))
1.28 + mod_up p)
1.29 + p
1.30 in
1.31 if deg_up g_p < deg_up g
1.32 then
1.33 - if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p'
1.34 + if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p
1.35 else
1.36 - if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p' else
1.37 + if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else
1.38 let
1.39 - val P = P * p'
1.40 - val g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P
1.41 - in try_new_prime_up a b d M P g p' end
1.42 + val P = P * p
1.43 + val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
1.44 + in try_new_prime_up a b d M P g p end
1.45 end
1.46
1.47 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
1.48 @@ -384,7 +384,7 @@
1.49 fun HENSEL_lifting_up a b d M p =
1.50 let
1.51 val p = p next_prime_not_dvd d
1.52 - val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p
1.53 + val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
1.54 in
1.55 if deg_up g_p = 0 then [1] else
1.56 let
2.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Tue May 21 11:33:03 2013 +0200
2.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Thu May 23 12:20:28 2013 +0200
2.3 @@ -69,13 +69,36 @@
2.4 value "mod_inv -3 7 = 2"
2.5 value "mod_inv -4 339 = 254"
2.6
2.7 -(* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat
2.8 +(* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
2.9 mod_div a b m = c
2.10 assumes True
2.11 - yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = *)
2.12 + yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
2.13 definition mod_div :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
2.14 "mod_div a b m = (nat a) * (mod_inv b m) mod m"
2.15
2.16 +(*if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";
2.17 +if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";
2.18 +if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";
2.19 +*)
2.20 +value "mod_div 21 3 5 = 2" value "(21::int) mod 5 = (3 * 2) mod 5"
2.21 +(* 21/3 = 7 mod 5 21 mod 5 = 6 mod 5
2.22 + = 2 1 1 *)
2.23 +value "mod_div 22 3 5 = 4" value "(22::int) mod 5 = (3 * 4) mod 5"
2.24 +(* 22/3 = ------- 22 mod 5 = 12 mod 5
2.25 + = 4 2 2 *)
2.26 +value "mod_div 23 3 5 = 1" value "(23::int) mod 5 = (3 * 1) mod 5"
2.27 +(* 23/3 = ------- 23 mod 5 = 3 mod 5
2.28 + = 1 3 3 *)
2.29 +value "mod_div 24 3 5 = 3" value "(24::int) mod 5 = (3 * 3) mod 5"
2.30 +(* 24/3 = ------- 24 mod 5 = 9 mod 5
2.31 + = 3 4 4 *)
2.32 +value "mod_div 25 3 5 = 0" value "(25::int) mod 5 = (3 * 0) mod 5"
2.33 +(* 25/3 = ------- 25 mod 5 = 0 mod 5
2.34 + = 0 0 0 *)
2.35 +value "mod_div 21 4 5 = 4" value "(21::int) mod 5 = (4 * 4) mod 5"
2.36 +value "mod_div 1 4 5 = 4" value "( 1::int) mod 5 = (4 * 4) mod 5"
2.37 +value "mod_div 0 4 5 = 0" value "( 0::int) mod 5 = (0 * 4) mod 5"
2.38 +
2.39 (* root1 is just local to approx_root *)
2.40 function root1 :: "int \<Rightarrow> nat \<Rightarrow> nat" where
2.41 "root1 a n = (if (int (n * n)) < a then root1 a (n + 1) else n)"
2.42 @@ -156,20 +179,41 @@
2.43 1 2 3 4 5 6
2.44 a: <= ? <= ? ? <=
2.45 b: <= ? <= <= <= <= *)
2.46 -termination make_primes (*not finished*)
2.47 +
2.48 +find_theorems "_ - _ < _ = _"
2.49 +find_theorems "?n < ?n + _"
2.50 + thm Rings.linordered_semidom_class.less_add_one --"?a < ?a + 1"
2.51 +find_theorems "_ + ?a < _ + ?a = _"
2.52 + thm Groups.ordered_ab_semigroup_add_imp_le_class.add_less_cancel_right
2.53 +find_theorems "_ - ?a < _ - ?a = _"
2.54 + thm Nat.less_diff_iff
2.55 +
2.56 +lemma termin_1: "n - Suc (length ps) < n - length ps" (*? \<not> True ?*) sorry
2.57 +
2.58 +termination make_primes (*by lexicographic_order +PROOF:2 GOLAS / size_change LOOPS*)
2.59 apply (relation "measures [\<lambda>(ps, last_p, n). n - (length ps),
2.60 \<lambda>(ps, last_p, n). n + 2 - last_p]")
2.61 apply auto
2.62 -sorry (*by lexicographic_order unfinished subgoals*)
2.63 +(*\<lambda>(ps, last_p, n). last_p - (length ps):
2.64 +goal (4 subgoals):
2.65 + 1. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc last_p - length ps < last_p - length ps \<Longrightarrow> Suc last_p - length ps = last_p - length ps
2.66 + 2. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc last_p - length ps < last_p - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p
2.67 + 3. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc (Suc last_p) - length ps < last_p - length ps \<Longrightarrow> Suc (Suc last_p) - length ps = last_p - length ps
2.68 + 4. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc (Suc last_p) - length ps < last_p - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p
2.69 +\<lambda>(ps, last_p, n). n - (length ps):
2.70 +goal (2 subgoals):
2.71 + 1. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> n - Suc (length ps) < n - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p
2.72 + 2. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> n - last_p < Suc (Suc n) - last_p*)
2.73 +sorry
2.74 declare make_primes.simps [simp del] -- "next_prime_not_dvd"
2.75
2.76 -value "make_primes [2, 3] 3 3 = [2, 3]"
2.77 -value "make_primes [2, 3] 3 4 = [2, 3, 5]"
2.78 -value "make_primes [2, 3] 3 5 = [2, 3, 5]"
2.79 -value "make_primes [2, 3] 3 6 = [2, 3, 5, 7]"
2.80 -value "make_primes [2, 3] 3 7 = [2, 3, 5, 7]"
2.81 -value "make_primes [2, 3] 3 8 = [2, 3, 5, 7, 11]"
2.82 -value "make_primes [2, 3] 3 9 = [2, 3, 5, 7, 11]"
2.83 +value "make_primes [2, 3] 3 3 = [2, 3]"
2.84 +value "make_primes [2, 3] 3 4 = [2, 3, 5]"
2.85 +value "make_primes [2, 3] 3 5 = [2, 3, 5]"
2.86 +value "make_primes [2, 3] 3 6 = [2, 3, 5, 7]"
2.87 +value "make_primes [2, 3] 3 7 = [2, 3, 5, 7]"
2.88 +value "make_primes [2, 3] 3 8 = [2, 3, 5, 7, 11]"
2.89 +value "make_primes [2, 3] 3 9 = [2, 3, 5, 7, 11]"
2.90 value "make_primes [2, 3] 3 10 = [2, 3, 5, 7, 11]"
2.91 value "make_primes [2, 3] 3 11 = [2, 3, 5, 7, 11]"
2.92 value "make_primes [2, 3] 3 12 = [2, 3, 5, 7, 11, 13]"
2.93 @@ -191,6 +235,8 @@
2.94 definition primes_upto :: "nat \<Rightarrow> nat list" where
2.95 "primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"
2.96
2.97 +value "primes_upto 0 = [2]"
2.98 +value "primes_upto 1 = [2]"
2.99 value "primes_upto 2 = [2]"
2.100 value "primes_upto 3 = [2, 3]"
2.101 value "primes_upto 4 = [2, 3, 5]"
2.102 @@ -202,6 +248,25 @@
2.103 value "primes_upto 10 = [2, 3, 5, 7, 11]"
2.104 value "primes_upto 11 = [2, 3, 5, 7, 11]"
2.105
2.106 +lemma primes_upto_0: "last (primes_upto n) > 0" (*see above*) sorry
2.107 +lemma primes_upto_1: "last (primes_upto (Suc n)) > n" (*see above*) sorry
2.108 +lemma primes_upto_2: "last (primes_upto n) >= n" (*see above*) sorry
2.109 +
2.110 +lemma termin_next_prime_not_dvd:
2.111 + shows "last (primes_upto (Suc n1)) * q - last (primes_upto (Suc n1))
2.112 + < last (primes_upto (Suc n1)) * q - n1"
2.113 +proof -
2.114 + from primes_upto_1 have "n1 < last (primes_upto (Suc n1))" by auto
2.115 + from this have "(int n1) - (int (last (primes_upto (Suc n1)) * q))
2.116 + < (int (last (primes_upto (Suc n1)))) - (int (last (primes_upto (Suc n1)) * q))" by auto
2.117 + from this show ?thesis sorry
2.118 +qed (*?FLORIAN : lifting nat -> int ?*)
2.119 +
2.120 +
2.121 +find_theorems "_ - _ < _ - _ = _" thm Nat.less_diff_iff
2.122 +find_theorems "_ * _ < _ * _ = _" thm Nat.Suc_mult_less_cancel1
2.123 +find_theorems "-_ < -_ = _" thm Groups.ordered_ab_group_add_class.neg_less_iff_less
2.124 +
2.125 (* max's' is analogous to Integer.gcds *)
2.126 definition maxs :: "nat list \<Rightarrow> nat" where "maxs ns = List.fold max ns (List.hd ns)"
2.127
2.128 @@ -210,7 +275,8 @@
2.129 (* find the next prime greater p not dividing the number n:
2.130 next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
2.131 n1 next_prime_not_dvd n2 = p
2.132 - assumes True
2.133 + assumes True assumes "0 < q"
2.134 +
2.135 yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
2.136 (\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
2.137 function next_prime_not_dvd :: "nat \<Rightarrow> nat \<Rightarrow> nat" (infixl "next'_prime'_not'_dvd" 70) where
2.138 @@ -223,19 +289,12 @@
2.139 then nxt
2.140 else nxt next_prime_not_dvd n2)"
2.141 by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare"
2.142 -termination (*next_prime_not_dvd ..Inner syntax error*) (*not finished*)
2.143 +termination (*next_prime_not_dvd: lexicographic_order +PROOF:?lifting nat-int? / size_change: Failed to apply initial proof method*)
2.144 apply (relation "measure (\<lambda>(n1, n2). n2 - n1)")
2.145 apply auto
2.146 -sorry (*by lexicographic_order unfinished subgoals*)
2.147 -(*Calls:
2.148 - a) (n1, n2) ~> (xa, n2)
2.149 -Measures:
2.150 - 1) \<lambda>p. size (snd p)
2.151 - 2) \<lambda>p. size (fst p)
2.152 - 3) size
2.153 -Result matrix:
2.154 - 1 2 3
2.155 -a: <= ? <= *)
2.156 + apply (rule termin_next_prime_not_dvd) (*?FLORIAN: IN Try_GCD... this is not necessary!?!+*)
2.157 +done
2.158 +
2.159 value "1 next_prime_not_dvd 15 = 2"
2.160 value "2 next_prime_not_dvd 15 = 7"
2.161 value "3 next_prime_not_dvd 15 = 7"
2.162 @@ -281,19 +340,19 @@
2.163 value "drop_lc0_up [0, 1, 2, 3, 4, 5] = [0, 1, 2, 3, 4, 5]"
2.164
2.165 (* degree *)
2.166 -(* THESE CREATE THE RESULT value "[-18,.." = value "[-1,.." = [1]
2.167 - ALTHOUGH THE TESTS BELOW SEEM TO WORK (*?FLORIAN*)
2.168 -(**)function deg_up :: "unipoly \<Rightarrow> nat" where
2.169 +(* THE VERSIONS BELOW CREATE THE RESULT value "[-18,.." --> value "[-1,.." = [1]
2.170 + ALTHOUGH THE TESTS BELOW SEEM TO WORK AND THE SAME CODE WORKS IN ML ?FLORIAN*)
2.171 +
2.172 +(*function deg_up :: "unipoly \<Rightarrow> nat" where
2.173 "deg_up p = ((op - 1) o length o drop_lc0_up) p"
2.174 -by auto
2.175 -termination sorry (* outcommented *)
2.176 +by auto termination sorry*)
2.177
2.178 -(**)fun deg_up :: "unipoly \<Rightarrow> nat" where
2.179 - "deg_up p = ((op - 1) o length o drop_lc0_up) p"
2.180 +(*fun deg_up :: "unipoly \<Rightarrow> nat" where
2.181 + "deg_up p = ((op - 1) o length o drop_lc0_up) p"*)
2.182
2.183 -(**)definition deg_up :: "unipoly \<Rightarrow> nat" where
2.184 - "deg_up p = ((op - 1) o length o drop_lc0_up) p"
2.185 -*)
2.186 +(*definition deg_up :: "unipoly \<Rightarrow> nat" where
2.187 + "deg_up p = ((op - 1) o length o drop_lc0_up) p"*)
2.188 +
2.189 function deg_up :: "unipoly \<Rightarrow> nat" where
2.190 "deg_up p =
2.191 (let len = List.length p - 1
2.192 @@ -301,21 +360,40 @@
2.193 if p ! len \<noteq> 0 then len else deg_up (nth_drop len p))"
2.194 by (*pat_completeness*) auto
2.195
2.196 -lemma termin_1:
2.197 - fixes p::"nat list"
2.198 - assumes "length p - Suc 0 = 0"
2.199 - shows "min (length p) (length p - Suc 0) < length p"
2.200 +lemma min_Suc:
2.201 + fixes a::nat
2.202 + assumes "0 < a"
2.203 + shows "min a (a - Suc 0) < a"
2.204 proof -
2.205 - from `length p - Suc 0 = 0` have "length p = Suc 0" sorry
2.206 - from `length p = Suc 0` have "min (length p) (length p - Suc 0) = 0" by auto
2.207 - from `length p = Suc 0` `min (length p) (length p - Suc 0) = 0` show ?thesis by auto
2.208 + from min_def have 1: "min a (a - Suc 0) = a - Suc 0" by auto
2.209 + from `0 < a` have 2: "a - Suc 0 < a" by auto
2.210 + from 1 2 show ?thesis by auto
2.211 qed
2.212 -thm termin_1 (*length ?p - Suc 0 = 0 \<Longrightarrow> min (length ?p) (length ?p - Suc 0) < length ?p*)
2.213 -termination deg_up (*not finished*)
2.214 +
2.215 +find_theorems "min _ _ = (_::nat)"
2.216 +find_theorems "min _ _ "
2.217 + thm Orderings.ord_class.min_def --"min ?a ?b = (if ?a \<le> ?b then ?a else ?b)"
2.218 +
2.219 +termination deg_up (*by lexicographic_order +PROOF:STRANGE GOAL+definition / by size_change: Failed to apply initial proof method*)
2.220 apply (relation "measure (\<lambda>(p). length p)")
2.221 apply auto
2.222 - (*apply (rule termin_1)*)
2.223 -sorry (*by lexicographic_order unfinished subgoals*)
2.224 +(*..RELATED?: https://lists.cam.ac.uk/mailman/htdig/cl-isabelle-users/2013-May/msg00075.html*)
2.225 + apply (rule min_Suc) (* 1. \<And>p. p ! (length p - Suc 0) = 0 \<Longrightarrow> 0 < length p*)
2.226 + apply auto (*?????: 1. [] ! 0 = 0 \<Longrightarrow> False *)
2.227 +sorry
2.228 +
2.229 +find_theorems "nth _ _"
2.230 + thm List.nth_mem (*?n < length ?xs \<Longrightarrow> ?xs ! ?n \<in> set ?xs*)
2.231 + (*(length p - Suc 0) < length ?xs *)
2.232 + thm List.set_conv_nth (*set ?xs = {?xs ! i |i. i < length ?xs}*)
2.233 +find_theorems "set _" (*found 378 theorem(s)*)
2.234 +find_theorems "length _" (*found 241 theorem(s)*)
2.235 +
2.236 +value "[1,2,3,4,5::int] ! 2"
2.237 +value "[1,2,3,4,5::int] ! 4"
2.238 +value "[1::int] ! 0"
2.239 +value "([]::int list) ! 0"
2.240 +
2.241 (*Calls:
2.242 a) p ~> nth_drop x p
2.243 Measures:
2.244 @@ -382,12 +460,12 @@
2.245 value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
2.246
2.247 (* lemmata for pattern compatibility in dvd_up *)
2.248 -lemma ex_fals0_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET" by simp
2.249 -lemma ex_fals0_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET" by simp
2.250 -lemma ex_fals0_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET" by simp
2.251 +lemma ex_falso_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET" by simp
2.252 +lemma ex_falso_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET" by simp
2.253 +lemma ex_falso_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET" by simp
2.254
2.255 function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where
2.256 -"[d] %|% [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) & (p mod d = 0))" |
2.257 +"[d] %|% [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) \<and> (p mod d = 0))" |
2.258 "ds %|% ps =
2.259 (let
2.260 ds = drop_lc0_up ds; ps = drop_lc0_up ps;
2.261 @@ -396,7 +474,7 @@
2.262 rest = drop_lc0_up (ps %-% (d000 %* quot))
2.263 in
2.264 if rest = [] then True else
2.265 - if quot \<noteq> 0 & List.length ds \<le> List.length rest then ds %|% rest else False)"
2.266 + if quot \<noteq> 0 \<and> List.length ds \<le> List.length rest then ds %|% rest else False)"
2.267 apply pat_completeness
2.268 apply simp
2.269 apply simp
2.270 @@ -417,18 +495,36 @@
2.271 apply simp (* > 1 sec IMPROVED BY declare simp del:
2.272 centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)
2.273 defer (*a "mixed" obligation*)
2.274 -apply (rule ex_fals0_1)
2.275 +apply (rule ex_falso_1)
2.276 apply simp
2.277 -apply (rule ex_fals0_2)
2.278 +apply (rule ex_falso_2)
2.279 apply simp
2.280 -apply (rule ex_fals0_3)
2.281 +apply (rule ex_falso_3)
2.282 apply simp
2.283 -apply (rule ex_fals0_1)
2.284 +apply (rule ex_falso_1)
2.285 apply simp
2.286 done (* > 1 sec IMPROVED BY declare simp del:
2.287 centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)
2.288 -termination (*dvd_up ..Inner syntax error*) sorry (*by lexicographic_order LOOPS*) (*not finished*)
2.289 -(*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
2.290 +termination (*dvd_up: by lexicographic_order LOOPS +PROOF:5 HUGE GOALS/ size_change LOOPS*)
2.291 +using [[linarith_split_limit = 999]]
2.292 +apply (relation "measure (\<lambda>(ds, ps). length (let
2.293 + ds = drop_lc0_up ds; ps = drop_lc0_up ps;
2.294 + d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
2.295 + quot = (lcoeff_up ps) div2 (lcoeff_up d000);
2.296 + rest = drop_lc0_up (ps %-% (d000 %* quot))
2.297 + in
2.298 + ds) - length ds)")
2.299 +apply auto
2.300 +(*apply (relation "measure (\<lambda>(ds, ps). length ps - length ds)"):
2.301 + 1. wf (measure (\<lambda>(ds, ps). length ps - length ds))
2.302 + 2. \<And>ps x xa xb xc xd. x = drop_lc0_up [] \<Longrightarrow> xa = drop_lc0_up ps \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), [], ps) \<in> measure (\<lambda>(ds, ps). length ps - length ds)
2.303 + 3. \<And>v vb vc ps x xa xb xc xd. x = drop_lc0_up (v # vb # vc) \<Longrightarrow> xa = drop_lc0_up ps \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), v # vb # vc, ps) \<in> measure (\<lambda>(ds, ps). length ps - length ds)
2.304 + 4. \<And>ds x xa xb xc xd. x = drop_lc0_up ds \<Longrightarrow> xa = drop_lc0_up [] \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), ds, []) \<in> measure (\<lambda>(ds, ps). length ps - length ds)
2.305 + 5. \<And>ds v vb vc x xa xb xc xd. x = drop_lc0_up ds \<Longrightarrow> xa = drop_lc0_up (v # vb # vc) \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), ds, v # vb # vc) \<in> measure (\<lambda>(ds, ps). length ps - length ds)
2.306 +apply auto
2.307 + 6 subgoals, even longer
2.308 +*)
2.309 +sorry
2.310
2.311 value "[4] %|% [6] = False"
2.312 value "[8] %|% [16, 0] = True"
2.313 @@ -505,7 +601,6 @@
2.314 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"
2.315
2.316 subsection {* modulo calculations for polynomials *}
2.317 -value "List.compound"
2.318
2.319 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)
2.320 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where
2.321 @@ -536,7 +631,7 @@
2.322 value "[5, 4, 7, 8, 1] mod_up 5 = [0, 4, 2, 3, 1]"
2.323 value "[5, 4,-7, 8,-1] mod_up 5 = [0, 4, 3, 3, 4]"
2.324
2.325 -(* euclidean algoritm in Z_p[x].
2.326 +(* euclidean algoritm in Z_p[x/m].
2.327 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly
2.328 mod_up_gcd p1 p2 m = g
2.329 assumes
2.330 @@ -550,28 +645,38 @@
2.331 quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
2.332 rest = drop_lc0_up ((p1m %-% (p2n %* (int quot))) mod_up m)
2.333 in
2.334 - if rest = []
2.335 - then p2
2.336 - else
2.337 + if rest = [] then p2 else
2.338 if List.length rest < List.length p2
2.339 then mod_up_gcd p2 rest m
2.340 else mod_up_gcd rest p2 m)"
2.341 -by auto
2.342 -termination mod_up_gcd (*not finished*)
2.343 -apply (relation "measures [\<lambda>(p1, p2, m). length p1,
2.344 - \<lambda>(p1, p2, m). 00000000000000]")
2.345 -apply auto
2.346 -sorry (*by lexicographic_order unfinished subgoals*)
2.347 +by auto
2.348 +termination mod_up_gcd (*by lexicographic_order +PROOF:3 HUGE SUGOALS / size_change LOOPS*)
2.349 +apply (relation "measures [\<lambda>(p1, p2, m). length p2 - length (let
2.350 + p1m = p1 mod_up m;
2.351 + p2m = drop_lc0_up (p2 mod_up m);
2.352 + p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
2.353 + quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
2.354 + rest = drop_lc0_up ((p1m %-% (p2n %* (int quot))) mod_up m)
2.355 + in rest),
2.356 + \<lambda>(p1, p2, m). length p1 - length (let
2.357 + p1m = p1 mod_up m;
2.358 + p2m = drop_lc0_up (p2 mod_up m);
2.359 + p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
2.360 + quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
2.361 + rest = drop_lc0_up ((p1m %-% (p2n %* (int quot))) mod_up m)
2.362 + in rest)]")
2.363 +(*apply auto ..LOOPS*)
2.364 +sorry
2.365 (*Calls:
2.366 a) (p1, p2, m) ~> (p2, xab, m)
2.367 b) (p1, p2, m) ~> (xab, p2, m)
2.368 Measures:
2.369 - 1) \<lambda>p. size (snd (snd p))
2.370 - 2) \<lambda>p. list_size (nat \<circ> abs) (fst (snd p))
2.371 - 3) \<lambda>p. length (fst (snd p))
2.372 - 4) \<lambda>p. size (snd p)
2.373 - 5) \<lambda>p. list_size (nat \<circ> abs) (fst p)
2.374 - 6) \<lambda>p. length (fst p)
2.375 + 1) \<lambda>p. size (snd (snd p)) m
2.376 + 2) \<lambda>p. list_size (nat \<circ> abs) (fst (snd p)) p2 ?(nat \<circ> abs) ...coeff?!?
2.377 + 3) \<lambda>p. length (fst (snd p)) p2
2.378 + 4) \<lambda>p. size (snd p) (p2, m)
2.379 + 5) \<lambda>p. list_size (nat \<circ> abs) (fst p) p1 ?(nat \<circ> abs) ...coeff?!?
2.380 + 6) \<lambda>p. length (fst p) p1
2.381 7) size
2.382 Result matrix:
2.383 1 2 3 4 5 6 7
2.384 @@ -582,6 +687,8 @@
2.385 value "mod_up_gcd [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] 7 = [2, 6, 0, 2, 6]"
2.386 value "mod_up_gcd [8, 28, 22, -11, -14, 1, 2] [2, 6, 0, 2, 6] 7 = [2, 6, 0, 2, 6]"
2.387 value "mod_up_gcd [20, 15, 8, 6] [8, -2, -2, 3] 2 = [0, 1]"
2.388 + value "[20, 15, 8, 6] %|% [0, 1] = False"
2.389 + value "[8, -2, -2, 3] %|% [0, 1] = False"
2.390
2.391 (* analogous to Integer.gcds
2.392 gcds :: int list \<Rightarrow> int
2.393 @@ -616,35 +723,44 @@
2.394 (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
2.395 try_new_prime_up p1 p2 d M P g p = new_g
2.396 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
2.397 - M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
2.398 + M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
2.399 p1 is primitive \<and> p2 is primitive
2.400 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
2.401
2.402 - argumentns "a b d M P g p" named according to [1] p.93 *)
2.403 + argumentns "a b d M P g p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
2.404 (*declare [[show_types]]*)
2.405 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
2.406 where
2.407 "try_new_prime_up a b d M P g p =
2.408 (if P > M then g else
2.409 - let p' = p next_prime_not_dvd d;
2.410 - g_p = centr_up ( ( (normalise (mod_up_gcd a b p') p')
2.411 - %* ((int d) mod (int p')))
2.412 - mod_up p')
2.413 - p'
2.414 + let p = p next_prime_not_dvd d;
2.415 + g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
2.416 + %* (int (d mod p)))
2.417 + mod_up p)
2.418 + p
2.419 in
2.420 if deg_up g_p < deg_up g
2.421 then
2.422 - if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p'
2.423 + if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p
2.424 else
2.425 - if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p' else
2.426 + if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p else
2.427 let
2.428 - P = P * p';
2.429 - g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P
2.430 + P = P * p;
2.431 + g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
2.432 in
2.433 - try_new_prime_up a b d M P g p')"
2.434 + try_new_prime_up a b d M P g p)"
2.435 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
2.436 -termination try_new_prime_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
2.437 -(*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
2.438 +termination try_new_prime_up (*by lexicographic_order LOOPS +PROOF:? / by size_change LOOPS*)
2.439 +(*apply (relation "measures
2.440 + [\<lambda>(a, b, d, M, P, g, p). ???,
2.441 + \<lambda>(a, b, d, M, P, g, p). ???,
2.442 + \<lambda>(a, b, d, M, P, g, p). ???]")
2.443 +Calls (manually):
2.444 + a) (a, b, d, M, P, g, p) ~> (a, b, d, M, p, g_p {length g_p <= length a b}, p' {> p})
2.445 + b) (a, b, d, M, P, g, p) ~> (a, b, d, M, P, g , p' {> p})
2.446 + c) (a, b, d, M, P, g, p) ~> (a, b, d, M, P {=P*p'}, g {length g <= length a b}, p' {> p})
2.447 +*)
2.448 +sorry
2.449
2.450 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
2.451 HENSEL_lifting_up p1 p2 d M p = gcd
2.452 @@ -653,12 +769,12 @@
2.453 p1 is primitive \<and> p2 is primitive
2.454 yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
2.455
2.456 - argumentns "a b d M p" named according to [1] p.93 *)
2.457 + argumentns "a b d M p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
2.458 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where
2.459 "HENSEL_lifting_up a b d M p =
2.460 (let
2.461 p = p next_prime_not_dvd d;
2.462 - g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p
2.463 + g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p (*see above*)
2.464 in
2.465 if deg_up g_p = 0 then [1] else
2.466 (let
2.467 @@ -666,8 +782,11 @@
2.468 in
2.469 if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"
2.470 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
2.471 -termination HENSEL_lifting_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
2.472 -(*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
2.473 +termination HENSEL_lifting_up (*by lexicographic_order LOOPS +PROOF:?next_prime_not_dvd / by size_change LOOPS*)
2.474 +sorry
2.475 +(*Calls (manually):
2.476 + a) (a, b, d, M, p) ~> (a, b, d, M, p {> p next_prime_not_dvd d})
2.477 +*)
2.478
2.479 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
2.480 gcd_up a b = c
2.481 @@ -684,7 +803,7 @@
2.482 termination by lexicographic_order (*works*)
2.483
2.484 value "gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
2.485 -(* gcd (-1 + x^2) (x + x^2) = (1 + x) *)
2.486 +(* gcd (-1 + x^2) (x + x^2) = (1 + x) ...*)
2.487 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
2.488
2.489 (*
3.1 --- a/test/Tools/isac/Knowledge/gcd_poly.sml Tue May 21 11:33:03 2013 +0200
3.2 +++ b/test/Tools/isac/Knowledge/gcd_poly.sml Thu May 23 12:20:28 2013 +0200
3.3 @@ -107,6 +107,12 @@
3.4 "----------- fun mod_div --------------------------------";
3.5 "----------- fun mod_div --------------------------------";
3.6 "----------- fun mod_div --------------------------------";
3.7 +if mod_div 21 3 5 = 2 then () else error "mod_div changed"; 21 mod 5 = (3 * 2) mod 5;
3.8 +if mod_div 22 3 5 = 4 then () else error "mod_div changed"; 22 mod 5 = (3 * 4) mod 5;
3.9 +if mod_div 23 3 5 = 1 then () else error "mod_div changed"; 23 mod 5 = (3 * 1) mod 5;
3.10 +if mod_div 24 3 5 = 3 then () else error "mod_div changed"; 24 mod 5 = (3 * 3) mod 5;
3.11 +if mod_div 25 3 5 = 0 then () else error "mod_div changed"; 25 mod 5 = (3 * 0) mod 5;
3.12 +
3.13 if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";
3.14 if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";
3.15 if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";
4.1 --- a/test/Tools/isac/Test_Some2.thy Tue May 21 11:33:03 2013 +0200
4.2 +++ b/test/Tools/isac/Test_Some2.thy Thu May 23 12:20:28 2013 +0200
4.3 @@ -17,15 +17,8 @@
4.4 [(1,[0,0,0]),(1,[1,1,0])]
4.5 *}
4.6 ML {*
4.7 -val (r1, m1, r2, m2) = (17, 9, 3, 4)
4.8 *}
4.9 ML {*
4.10 -(r1 mod m1);
4.11 -r2 - (r1 mod m1);
4.12 -mod_inv m1 m2;
4.13 -(r2 - (r1 mod m1)) * (mod_inv m1 m2);
4.14 -((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2;
4.15 -(((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) * m1;
4.16 *}
4.17 ML {*
4.18 *}