1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Tue May 21 10:38:16 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Tue May 21 11:33:03 2013 +0200
1.3 @@ -33,25 +33,21 @@
1.4 value "gcd (10::int) (3::int) = 1"
1.5 value "gcd (6::int) (24::int) = 6"
1.6
1.7 -(* why has this been removed since 2002 ? *)
1.8 -definition last_elem :: "'a list \<Rightarrow> 'a" where
1.9 - "last_elem ls = (hd o rev) ls"
1.10 -
1.11 (* drop tail elements equal 0 *)
1.12 primrec drop_hd_zeros :: "int list \<Rightarrow> int list" where
1.13 - "drop_hd_zeros (p # ps) = (if p = 0 then drop_hd_zeros ps else (p # ps))"
1.14 +"drop_hd_zeros (p # ps) = (if p = 0 then drop_hd_zeros ps else (p # ps))"
1.15 definition drop_tl_zeros :: "int list \<Rightarrow> int list" where
1.16 - "drop_tl_zeros p = (rev o drop_hd_zeros o rev) p"
1.17 +"drop_tl_zeros p = (rev o drop_hd_zeros o rev) p"
1.18
1.19 subsection {* modulo calculations for integers *}
1.20 (* modi is just local for mod_inv *)
1.21 function modi :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where
1.22 - "modi r rold m rinv =
1.23 - (if m \<le> rinv
1.24 - then 0 else
1.25 - if r mod (int m) = 1
1.26 - then rinv
1.27 - else modi (rold * ((int rinv) + 1)) rold m (rinv + 1))"
1.28 +"modi r rold m rinv =
1.29 +(if m \<le> rinv
1.30 + then 0 else
1.31 + if r mod (int m) = 1
1.32 + then rinv
1.33 + else modi (rold * ((int rinv) + 1)) rold m (rinv + 1))"
1.34 by auto
1.35 termination by (relation "measure (\<lambda>(r, rold, m, rinv). m - rinv)") auto
1.36
1.37 @@ -59,8 +55,7 @@
1.38 mod_inv r m = s
1.39 assumes True
1.40 yields r * s mod m = 1 *)
1.41 -definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where
1.42 - "mod_inv r m = modi r r m 1"
1.43 +definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where "mod_inv r m = modi r r m 1"
1.44
1.45 value "modi 5 5 7 1 = 3"
1.46 value "modi 3 3 7 1 = 5"
1.47 @@ -79,11 +74,11 @@
1.48 assumes True
1.49 yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = *)
1.50 definition mod_div :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
1.51 - "mod_div a b m = (nat a) * (mod_inv b m) mod m"
1.52 +"mod_div a b m = (nat a) * (mod_inv b m) mod m"
1.53
1.54 (* root1 is just local to approx_root *)
1.55 function root1 :: "int \<Rightarrow> nat \<Rightarrow> nat" where
1.56 - "root1 a n = (if (int (n * n)) < a then root1 a (n + 1) else n)"
1.57 +"root1 a n = (if (int (n * n)) < a then root1 a (n + 1) else n)"
1.58 by auto
1.59 termination by (relation "measure (\<lambda>(a, n). nat (a - (int (n * n))))") auto
1.60
1.61 @@ -92,17 +87,16 @@
1.62 approx_root a = r
1.63 assumes True
1.64 yields r * r \<ge> a *)
1.65 -definition approx_root :: "int \<Rightarrow> nat" where
1.66 - "approx_root a = root1 a 1"
1.67 +definition approx_root :: "int \<Rightarrow> nat" where "approx_root a = root1 a 1"
1.68
1.69 (* chinese_remainder :: int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int
1.70 chinese_remainder (r1, m1) (r2, m2) = r
1.71 assumes True
1.72 yields r = r1 mod m1 \<and> r = r2 mod m2 *)
1.73 definition chinese_remainder :: "int \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
1.74 - "chinese_remainder r1 m1 r2 m2 =
1.75 - ((nat (r1 mod (int m1))) +
1.76 - (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv (int m1) m2))) mod (int m2))) * m1)"
1.77 +"chinese_remainder r1 m1 r2 m2 =
1.78 + ((nat (r1 mod (int m1))) +
1.79 + (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv (int m1) m2))) mod (int m2))) * m1)"
1.80
1.81 value "chinese_remainder 17 9 3 4 = 35"
1.82 value "chinese_remainder 7 2 6 11 = 17"
1.83 @@ -116,13 +110,13 @@
1.84 (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and> (\<forall>p. p < n \<and> Primes.prime p \<longrightarrow> List.member ps p)
1.85 yields Primes.prime n *)
1.86 fun is_prime :: "nat list \<Rightarrow> nat \<Rightarrow> bool" where
1.87 - "is_prime ps n =
1.88 - (if List.length ps > 0
1.89 - then
1.90 - if (n mod (List.hd ps)) = 0
1.91 - then False
1.92 - else is_prime (List.tl ps) n
1.93 - else True)"
1.94 +"is_prime ps n =
1.95 +(if List.length ps > 0
1.96 +then
1.97 + if (n mod (List.hd ps)) = 0
1.98 + then False
1.99 + else is_prime (List.tl ps) n
1.100 +else True)"
1.101 declare is_prime.simps [simp del] -- "make_primes, next_prime_not_dvd"
1.102
1.103 value "is_prime [2, 3] 2 = False" --"... precondition!"
1.104 @@ -142,11 +136,11 @@
1.105 yields n \<le> maxs pps \<and> (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p) \<and>
1.106 (\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
1.107 function make_primes :: "nat list \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat list" where
1.108 - "make_primes ps last_p n =
1.109 - (if n <= last_elem ps then ps else
1.110 - if is_prime ps (last_p + 2)
1.111 - then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
1.112 - else make_primes ps (last_p + 2) n)"
1.113 +"make_primes ps last_p n =
1.114 +(if n <= last ps then ps else
1.115 + if is_prime ps (last_p + 2)
1.116 + then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
1.117 + else make_primes ps (last_p + 2) n)"
1.118 by pat_completeness auto --"simp del: is_prime.simps <-- declare"
1.119 (*Calls:
1.120 a) (ps, last_p, n) ~> (ps @ [last_p + 2], last_p + 2, n)
1.121 @@ -195,7 +189,7 @@
1.122 n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
1.123 (\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
1.124 definition primes_upto :: "nat \<Rightarrow> nat list" where
1.125 - "primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"
1.126 +"primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"
1.127
1.128 value "primes_upto 2 = [2]"
1.129 value "primes_upto 3 = [2, 3]"
1.130 @@ -209,8 +203,8 @@
1.131 value "primes_upto 11 = [2, 3, 5, 7, 11]"
1.132
1.133 (* max's' is analogous to Integer.gcds *)
1.134 -definition maxs :: "nat list \<Rightarrow> nat" where
1.135 - "maxs ns = List.fold max ns (List.hd ns)"
1.136 +definition maxs :: "nat list \<Rightarrow> nat" where "maxs ns = List.fold max ns (List.hd ns)"
1.137 +
1.138 value "maxs [5, 3, 7, 1, 2, 4] = 7"
1.139
1.140 (* find the next prime greater p not dividing the number n:
1.141 @@ -220,14 +214,14 @@
1.142 yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
1.143 (\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
1.144 function next_prime_not_dvd :: "nat \<Rightarrow> nat \<Rightarrow> nat" (infixl "next'_prime'_not'_dvd" 70) where
1.145 - "n1 next_prime_not_dvd n2 =
1.146 - (let
1.147 - ps = primes_upto (n1 + 1);
1.148 - nxt = last ps
1.149 - in
1.150 - if n2 mod nxt \<noteq> 0
1.151 - then nxt
1.152 - else nxt next_prime_not_dvd n2)"
1.153 +"n1 next_prime_not_dvd n2 =
1.154 +(let
1.155 + ps = primes_upto (n1 + 1);
1.156 + nxt = last ps
1.157 +in
1.158 + if n2 mod nxt \<noteq> 0
1.159 + then nxt
1.160 + else nxt next_prime_not_dvd n2)"
1.161 by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare"
1.162 termination (*next_prime_not_dvd ..Inner syntax error*) (*not finished*)
1.163 apply (relation "measure (\<lambda>(n1, n2). n2 - n1)")
1.164 @@ -254,15 +248,14 @@
1.165
1.166 (* not in List.thy, copy from library.ML *)
1.167 fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list" where
1.168 - "nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"
1.169 +"nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"
1.170 value "nth_drop 0 [] = []"
1.171 value "nth_drop 0 [1, 2, 3::int] = [2, 3]"
1.172 value "nth_drop 2 [1, 2, 3::int] = [1, 2]"
1.173 value "nth_drop 77 [1, 2, 3::int] = [1, 2, 3]"
1.174
1.175 (* leading coefficient *)
1.176 -definition lcoeff_up :: "unipoly \<Rightarrow> int" where
1.177 - "lcoeff_up p = (last_elem o drop_tl_zeros) p"
1.178 +definition lcoeff_up :: "unipoly \<Rightarrow> int" where "lcoeff_up p = (last o drop_tl_zeros) p"
1.179
1.180 value "lcoeff_up [3, 4, 5, 6] = 6"
1.181 value "lcoeff_up [3, 4, 5, 6, 0] = 6"
1.182 @@ -276,13 +269,13 @@
1.183 "drop_lc0_up p = drop_tl_zeros p"
1.184 *)
1.185 fun drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where
1.186 - "drop_lc0_up [] = []"
1.187 -| "drop_lc0_up p =
1.188 - (let l = List.length p - 1
1.189 - in
1.190 - if p ! l = 0
1.191 - then drop_lc0_up (nth_drop l p)
1.192 - else p)"
1.193 +"drop_lc0_up [] = []" |
1.194 +"drop_lc0_up p =
1.195 + (let l = List.length p - 1
1.196 + in
1.197 + if p ! l = 0
1.198 + then drop_lc0_up (nth_drop l p)
1.199 + else p)"
1.200
1.201 value "drop_lc0_up [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5]"
1.202 value "drop_lc0_up [0, 1, 2, 3, 4, 5] = [0, 1, 2, 3, 4, 5]"
1.203 @@ -302,10 +295,10 @@
1.204 "deg_up p = ((op - 1) o length o drop_lc0_up) p"
1.205 *)
1.206 function deg_up :: "unipoly \<Rightarrow> nat" where
1.207 - "deg_up p =
1.208 - (let len = List.length p - 1
1.209 - in
1.210 - if p ! len \<noteq> 0 then len else deg_up (nth_drop len p))"
1.211 +"deg_up p =
1.212 +(let len = List.length p - 1
1.213 +in
1.214 + if p ! len \<noteq> 0 then len else deg_up (nth_drop len p))"
1.215 by (*pat_completeness*) auto
1.216
1.217 lemma termin_1:
1.218 @@ -337,22 +330,22 @@
1.219
1.220 (* norm is just local to normalise *)
1.221 fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where
1.222 - "norm p nrm m lcp i =
1.223 - (if i = 0
1.224 - then [int (mod_div (p ! i) lcp m)] @ nrm
1.225 - else norm p ([int (mod_div (p ! i) lcp m)] @ nrm) m lcp (i - 1))"
1.226 +"norm p nrm m lcp i =
1.227 +(if i = 0
1.228 + then [int (mod_div (p ! i) lcp m)] @ nrm
1.229 + else norm p ([int (mod_div (p ! i) lcp m)] @ nrm) m lcp (i - 1))"
1.230 (* normalise a unipoly such that lcoeff_up mod m = 1.
1.231 normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
1.232 normalise [p_0, .., p_k] m = [q_0, .., q_k]
1.233 assumes
1.234 yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i \<and> (p_k * t) mod m = 1 *)
1.235 fun normalise :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
1.236 - "normalise p m =
1.237 - (let
1.238 - p = drop_lc0_up p;
1.239 - lcp = lcoeff_up p
1.240 - in
1.241 - if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))"
1.242 +"normalise p m =
1.243 +(let
1.244 + p = drop_lc0_up p;
1.245 + lcp = lcoeff_up p
1.246 +in
1.247 + if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))"
1.248 declare normalise.simps [simp del] --"HENSEL_lifting_up"
1.249
1.250 value "normalise [-18, -15, -20, 12, 20, -13, 2] 5 = [1, 0, 0, 1, 0, 1, 1]"
1.251 @@ -360,54 +353,50 @@
1.252
1.253 subsection {* addition, multiplication, division *}
1.254
1.255 -(* scalar multiplication, %* in GCD_Poly.thy *)
1.256 -definition mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "mult'_ups" 70) where
1.257 - "p mult_ups m = List.map (op * m) p"
1.258 +(* scalar multiplication *)
1.259 +definition mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
1.260 +"p %* m = List.map (op * m) p"
1.261
1.262 -value "[5, 4, 7, 8, 1] mult_ups 5 = [25, 20, 35, 40, 5]"
1.263 -value "[5, 4, -7, 8, -1] mult_ups 5 = [25, 20, -35, 40, -5]"
1.264 +value "[5, 4, 7, 8, 1] %* 5 = [25, 20, 35, 40, 5]"
1.265 +value "[5, 4, -7, 8, -1] %* 5 = [25, 20, -35, 40, -5]"
1.266
1.267 -(* scalar divison, %/ in GCD_Poly.thy *)
1.268 -definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where
1.269 - "swapf f a b = f b a"
1.270 -definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" 70) where
1.271 - "p div_ups m = map (swapf op div2 m) p"
1.272 +(* scalar divison *)
1.273 +definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
1.274 +definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" (* %/ error FLORIAN*) 70) where
1.275 +"p div_ups m = map (swapf op div2 m) p"
1.276
1.277 value "[4, 3, 2, 5, 6] div_ups 3 = [1, 1, 0, 1, 2]"
1.278 value "[4, 3, 2, 0] div_ups 3 = [1, 1, 0, 0]"
1.279
1.280 (* not in List.thy, copy from library.ML *)
1.281 fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where
1.282 - "map2 _ [] [] = []"
1.283 -| "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys"
1.284 -| "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
1.285 +"map2 _ [] [] = []" |
1.286 +"map2 f (x # xs) (y # ys) = f x y # map2 f xs ys" |
1.287 +"map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
1.288
1.289 -(* minus of polys, %-% in GCD_Poly.thy *)
1.290 -definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "minus'_up" 70) where
1.291 - "p1 minus_up p2 = map2 (op -) p1 p2"
1.292 +(* minus of polys *)
1.293 +definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "%-%" 70) where
1.294 +"p1 %-% p2 = map2 (op -) p1 p2"
1.295
1.296 -value "[8, -7, 0, 1] minus_up [-2, 2, 3, 0] = [10, -9, -3, 1]"
1.297 -value "[8, 7, 6, 5, 4] minus_up [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
1.298 +value "[8, -7, 0, 1] %-% [-2, 2, 3, 0] = [10, -9, -3, 1]"
1.299 +value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
1.300
1.301 (* lemmata for pattern compatibility in dvd_up *)
1.302 -lemma ex_fals0_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET"
1.303 -by simp
1.304 -lemma ex_fals0_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET"
1.305 -by simp
1.306 -lemma ex_fals0_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET"
1.307 -by simp
1.308 +lemma ex_fals0_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET" by simp
1.309 +lemma ex_fals0_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET" by simp
1.310 +lemma ex_fals0_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET" by simp
1.311
1.312 -function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "dvd'_up" 70) where
1.313 - "[d] dvd_up [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) & (p mod d = 0))"
1.314 -| "ds dvd_up ps =
1.315 - (let
1.316 - ds = drop_lc0_up ds; ps = drop_lc0_up ps;
1.317 - d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
1.318 - quot = (lcoeff_up ps) div2 (lcoeff_up d000);
1.319 - rest = drop_lc0_up (ps minus_up (d000 mult_ups quot))
1.320 - in
1.321 - if rest = [] then True else
1.322 - if quot \<noteq> 0 & List.length ds \<le> List.length rest then ds dvd_up rest else False)"
1.323 +function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where
1.324 +"[d] %|% [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) & (p mod d = 0))" |
1.325 +"ds %|% ps =
1.326 + (let
1.327 + ds = drop_lc0_up ds; ps = drop_lc0_up ps;
1.328 + d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
1.329 + quot = (lcoeff_up ps) div2 (lcoeff_up d000);
1.330 + rest = drop_lc0_up (ps %-% (d000 %* quot))
1.331 + in
1.332 + if rest = [] then True else
1.333 + if quot \<noteq> 0 & List.length ds \<le> List.length rest then ds %|% rest else False)"
1.334 apply pat_completeness
1.335 apply simp
1.336 apply simp
1.337 @@ -441,16 +430,16 @@
1.338 termination (*dvd_up ..Inner syntax error*) sorry (*by lexicographic_order LOOPS*) (*not finished*)
1.339 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
1.340
1.341 -value "[4] dvd_up [6] = False"
1.342 -value "[8] dvd_up [16, 0] = True"
1.343 -value "[3, 2] dvd_up [0, 0, 9, 12, 4] = True"
1.344 -value "[8, 0] dvd_up [16] = True"
1.345 +value "[4] %|% [6] = False"
1.346 +value "[8] %|% [16, 0] = True"
1.347 +value "[3, 2] %|% [0, 0, 9, 12, 4] = True"
1.348 +value "[8, 0] %|% [16] = True"
1.349
1.350 subsection {* normalisation and Landau-Mignotte bound *}
1.351
1.352 (* centr is just local to centr_up *)
1.353 definition centr :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where
1.354 - "centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"
1.355 +"centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"
1.356
1.357 (* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
1.358 normalise [p_0, .., p_k] m = [q_0, .., q_k]
1.359 @@ -458,11 +447,11 @@
1.360 yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
1.361 (where |^ x ^| means round x up to the next greater integer) *)
1.362 definition centr_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
1.363 - "centr_up p m =
1.364 - (let
1.365 - mi = (int m) div2 2;
1.366 - mid = if (int m) mod mi = 0 then mi else mi + 1
1.367 - in map (centr m mid) p)"
1.368 +"centr_up p m =
1.369 +(let
1.370 + mi = (int m) div2 2;
1.371 + mid = if (int m) mod mi = 0 then mi else mi + 1
1.372 +in map (centr m mid) p)"
1.373
1.374 value "centr_up [7, 3, 5, 8, 1, 3] 10 = [-3, 3, 5, -2, 1, 3]"
1.375 value "centr_up [1, 2, 3, 4, 5] 2 = [1, 0, 1, 2, 3]"
1.376 @@ -476,12 +465,12 @@
1.377 value "centr_up [1, 2, 3, 4, 5] 10 = [1, 2, 3, 4, 5]"
1.378 value "centr_up [1, 2, 3, 4, 5] 11 = [1, 2, 3, 4, 5]"
1.379
1.380 - (* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
1.381 +(* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
1.382 sum_lmb [p_0, .., p_k] e = s
1.383 assumes exp >= 0
1.384 yields. p_0^e + p_1^e + ... + p_k^e *)
1.385 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where
1.386 - "sum_lmb p exp = List.fold ((op +) o (swapf power exp)) p 0"
1.387 +"sum_lmb p exp = List.fold ((op +) o (swapf power exp)) p 0"
1.388
1.389 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"
1.390 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"
1.391 @@ -496,14 +485,10 @@
1.392 yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
1.393 min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)
1.394 definition LANDAU_MIGNOTTE_bound :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat" where
1.395 - "LANDAU_MIGNOTTE_bound p1 p2 =
1.396 - ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) *
1.397 - (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))
1.398 - (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"
1.399 -(* |------------| -|------------|
1.400 - |--------------------------|
1.401 - |--------------------------------|
1.402 - |-------------------------------------------------------| *)
1.403 +"LANDAU_MIGNOTTE_bound p1 p2 =
1.404 + ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) *
1.405 + (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))
1.406 + (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"
1.407
1.408 value "LANDAU_MIGNOTTE_bound [1] [4, 5] = 1"
1.409 value "LANDAU_MIGNOTTE_bound [1, 2] [4, 5] = 2"
1.410 @@ -520,22 +505,23 @@
1.411 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"
1.412
1.413 subsection {* modulo calculations for polynomials *}
1.414 +value "List.compound"
1.415
1.416 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)
1.417 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where
1.418 - "([] pair []) = []"
1.419 -| "([] pair ys) = []" (*raise ListPair.UnequalLengths*)
1.420 -| "(xs pair []) = []" (*raise ListPair.UnequalLengths*)
1.421 -| "((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
1.422 +"([] pair []) = []" |
1.423 +"([] pair ys) = []" | (*raise ListPair.UnequalLengths*)
1.424 +"(xs pair []) = []" | (*raise ListPair.UnequalLengths*)
1.425 +"((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
1.426 fun chinese_rem :: "nat \<times> nat \<Rightarrow> int \<times> int \<Rightarrow> int" where
1.427 - "chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"
1.428 +"chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"
1.429
1.430 (* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
1.431 chinese_remainder_up (m1, m2) (p1, p2) = p
1.432 assume m1, m2 relatively prime
1.433 yields p1 = p mod m1 \<and> p2 = p mod m2 *)
1.434 fun chinese_remainder_up :: "nat \<times> nat \<Rightarrow> unipoly \<times> unipoly \<Rightarrow> unipoly" where
1.435 - "chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"
1.436 +"chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"
1.437
1.438 value "chinese_remainder_up (5, 7) ([2, 2, 4, 3], [3, 2, 3, 5]) = [17, 2, 24, 33]"
1.439
1.440 @@ -545,36 +531,31 @@
1.441 yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
1.442 definition mod' :: "nat \<Rightarrow> int \<Rightarrow> int" where "mod' m i = i mod (int m)"
1.443 definition mod_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" (infixl "mod'_up" 70) where
1.444 - "p mod_up m = map (mod' m) p"
1.445 +"p mod_up m = map (mod' m) p"
1.446
1.447 value "[5, 4, 7, 8, 1] mod_up 5 = [0, 4, 2, 3, 1]"
1.448 value "[5, 4,-7, 8,-1] mod_up 5 = [0, 4, 3, 3, 4]"
1.449
1.450 -value "curry"
1.451 -value "curry"
1.452 -value "5 mod 3::int"
1.453 -
1.454 -
1.455 (* euclidean algoritm in Z_p[x].
1.456 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly
1.457 mod_up_gcd p1 p2 m = g
1.458 assumes
1.459 yields gcd (p1 mod m) (p2 mod m) = g *)
1.460 function mod_up_gcd :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
1.461 - "mod_up_gcd p1 p2 m =
1.462 - (let
1.463 - p1m = p1 mod_up m;
1.464 - p2m = drop_lc0_up (p2 mod_up m);
1.465 - p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
1.466 - quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
1.467 - rest = drop_lc0_up ((p1m minus_up (p2n mult_ups (int quot))) mod_up m)
1.468 - in
1.469 - if rest = []
1.470 - then p2
1.471 - else
1.472 - if List.length rest < List.length p2
1.473 - then mod_up_gcd p2 rest m
1.474 - else mod_up_gcd rest p2 m)"
1.475 +"mod_up_gcd p1 p2 m =
1.476 +(let
1.477 + p1m = p1 mod_up m;
1.478 + p2m = drop_lc0_up (p2 mod_up m);
1.479 + p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
1.480 + quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
1.481 + rest = drop_lc0_up ((p1m %-% (p2n %* (int quot))) mod_up m)
1.482 +in
1.483 + if rest = []
1.484 + then p2
1.485 + else
1.486 + if List.length rest < List.length p2
1.487 + then mod_up_gcd p2 rest m
1.488 + else mod_up_gcd rest p2 m)"
1.489 by auto
1.490 termination mod_up_gcd (*not finished*)
1.491 apply (relation "measures [\<lambda>(p1, p2, m). length p1,
1.492 @@ -608,8 +589,7 @@
1.493 assumes True
1.494 yields THE d. ((\<forall>n. List.member ns n \<longrightarrow> d dvd n) \<and>
1.495 (\<forall>d'. (\<forall>n. List.member ns n \<and> d' dvd n) \<longrightarrow> d'modp \<le> d)) *)
1.496 -fun gcds :: "int list \<Rightarrow> int" where
1.497 - "gcds ns = List.fold gcd ns (List.hd ns)"
1.498 +fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)"
1.499
1.500 value "gcds [6, 9, 12] = 3"
1.501 value "gcds [6, -9, 12] = 3"
1.502 @@ -621,11 +601,11 @@
1.503 assumes True
1.504 yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
1.505 fun primitive_up :: "unipoly \<Rightarrow> unipoly" where
1.506 - "primitive_up [c] = (if c = 0 then [0] else [1])"
1.507 -| "primitive_up p =
1.508 - (let d = gcds p
1.509 - in
1.510 - if d = 1 then p else p div_ups d)"
1.511 +"primitive_up [c] = (if c = 0 then [0] else [1])" |
1.512 +"primitive_up p =
1.513 + (let d = gcds p
1.514 + in
1.515 + if d = 1 then p else p div_ups d)"
1.516
1.517 value "primitive_up [12, 16, 32, 44] = [3, 4, 8, 11]"
1.518 value "primitive_up [4, 5, 12] = [4, 5, 12]"
1.519 @@ -643,26 +623,25 @@
1.520 argumentns "a b d M P g p" named according to [1] p.93 *)
1.521 (*declare [[show_types]]*)
1.522 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
1.523 - where
1.524 - "try_new_prime_up a b d M P g p =
1.525 - (if P > M then g else
1.526 - let p' = p next_prime_not_dvd d;
1.527 - g_p = centr_up ( ( (normalise (mod_up_gcd a b p') p')
1.528 - mult_ups ((int d) mod (int p')))
1.529 - mod_up p')
1.530 - p'
1.531 - in
1.532 - if deg_up g_p < deg_up g
1.533 - then
1.534 - if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p'
1.535 - else
1.536 - if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p' else
1.537 - let
1.538 - P = P * p';
1.539 - g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P
1.540 - in
1.541 - try_new_prime_up a b d M P g p'
1.542 - )"
1.543 +where
1.544 +"try_new_prime_up a b d M P g p =
1.545 +(if P > M then g else
1.546 + let p' = p next_prime_not_dvd d;
1.547 + g_p = centr_up ( ( (normalise (mod_up_gcd a b p') p')
1.548 + %* ((int d) mod (int p')))
1.549 + mod_up p')
1.550 + p'
1.551 + in
1.552 + if deg_up g_p < deg_up g
1.553 + then
1.554 + if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p'
1.555 + else
1.556 + if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p' else
1.557 + let
1.558 + P = P * p';
1.559 + g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P
1.560 + in
1.561 + try_new_prime_up a b d M P g p')"
1.562 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
1.563 termination try_new_prime_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
1.564 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
1.565 @@ -676,16 +655,16 @@
1.566
1.567 argumentns "a b d M p" named according to [1] p.93 *)
1.568 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where
1.569 - "HENSEL_lifting_up a b d M p =
1.570 - (let
1.571 - p = p next_prime_not_dvd d;
1.572 - g_p = centr_up (((normalise (mod_up_gcd a b p) p) mult_ups (int (d mod p))) mod_up p) p
1.573 - in
1.574 - if deg_up g_p = 0 then [1] else
1.575 - (let
1.576 - g = primitive_up (try_new_prime_up a b d M p g_p p)
1.577 - in
1.578 - if (g dvd_up a) \<and> (g dvd_up b) then g else HENSEL_lifting_up a b d M p))"
1.579 +"HENSEL_lifting_up a b d M p =
1.580 +(let
1.581 + p = p next_prime_not_dvd d;
1.582 + g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p
1.583 +in
1.584 + if deg_up g_p = 0 then [1] else
1.585 + (let
1.586 + g = primitive_up (try_new_prime_up a b d M p g_p p)
1.587 + in
1.588 + if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"
1.589 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
1.590 termination HENSEL_lifting_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
1.591 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
1.592 @@ -696,11 +675,11 @@
1.593 a is primitive \<and> b is primitive
1.594 yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
1.595 function gcd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" where
1.596 - "gcd_up a b =
1.597 - (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>
1.598 - in
1.599 - if a = b then a else
1.600 - HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)"
1.601 +"gcd_up a b =
1.602 +(let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>
1.603 +in
1.604 + if a = b then a else
1.605 + HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)"
1.606 by pat_completeness auto --"simp del: lcoeff_up.simps ?+ others?"
1.607 termination by lexicographic_order (*works*)
1.608
1.609 @@ -708,13 +687,6 @@
1.610 (* gcd (-1 + x^2) (x + x^2) = (1 + x) *)
1.611 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
1.612
1.613 -primrec nth :: "'a list => nat => 'a" (infixl "!!" 100) where
1.614 -nth_Cons: "(x # xs) !! n = (case n of 0 \<Rightarrow> x | Suc k \<Rightarrow> xs !! k)"
1.615 -value "[1,2,3,4,5] ! 2"
1.616 -
1.617 -definition mult_ups2 :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
1.618 - "p %* m = List.map (op * m) p"
1.619 -
1.620 (*
1.621 print_configs
1.622 declare [[simp_trace_depth_limit = 99]]