1 (* GCD for polynomials by the function package following GCD_Poly_ML *)
3 imports "HOL-Computational_Algebra.Polynomial"
4 "HOL-Computational_Algebra.Primes"
8 This code has been translated from GCD_Poly.thy by Diana Meindl,
9 who follows Franz Winkler, Polynomial algorithyms in computer algebra, Springer 1996.
10 Winkler's original identifiers are in test/./gcd_poly_winkler.sml;
11 test/../gcd_poly.sml documents the changes towards GCD_Poly.thy;
12 the style of GCD_Poly.thy has been adapted to the function package.
15 section {* Isabelle's predefined polynomials *}
16 \<comment> \<open>TODO\<close>
18 section {* gcd for univariate polynomials *}
20 type_synonym unipoly = "int list" (*TODO: compare Polynomial.thy*)
21 value "[0, 1, 2, 3, 4, 5] :: unipoly"
23 subsection {* auxiliary functions *}
25 5 div 2 = 2; ~5 div 2 = ~3; BUT WEE NEED ~5 div2 2 = ~2; *)
26 definition div2 :: "int \<Rightarrow> int \<Rightarrow> int" (infixl "div2" 70) where
27 "a div2 b = (if a div b < 0 then (\<bar>a\<bar> div \<bar>b\<bar>) * -1 else a div b)"
30 value "-5 div2 2 = -2"
31 value "-5 div2 -2 = 2"
32 value " 5 div2 -2 = -2"
34 value "gcd (15::int) (6::int) = 3"
35 value "gcd (10::int) (3::int) = 1"
36 value "gcd (6::int) (24::int) = 6"
38 (* drop tail elements equal 0 *)
39 primrec drop_hd_zeros :: "int list \<Rightarrow> int list" where
40 "drop_hd_zeros [] = []" |
41 "drop_hd_zeros (p # ps) = (if p = 0 then drop_hd_zeros ps else (p # ps))"
43 (* drop leading coefficients equal 0 *)
44 definition drop_tl_zeros :: "int list \<Rightarrow> int list" where
45 "drop_tl_zeros = rev o drop_hd_zeros o rev"
46 value "drop_tl_zeros [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5]"
47 value "drop_tl_zeros [0, 1, 2, 3, 4, 5] = [0, 1, 2, 3, 4, 5]"
49 subsection {* modulo calculations for integers *}
50 (* modi is just local for mod_inv *)
51 function modi :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where
53 (if m \<le> rinv then 0 else
56 else modi (rold * ((int rinv) + 1)) rold m (rinv + 1))"
58 termination by (relation "measure (\<lambda>(r, rold, m, rinv). m - rinv)") auto
60 (* mod_inv :: int \<Rightarrow> nat \<Rightarrow> nat
63 yields r * s mod m = 1 *)
64 definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where "mod_inv r m = modi r r m 1"
66 value "modi 5 5 7 1 = 3"
67 value "modi 3 3 7 1 = 5"
68 value "modi 4 4 339 1 = 85"
70 value "mod_inv 5 7 = 3"
71 value "mod_inv 3 7 = 5"
72 value "mod_inv 4 339 = 85"
74 value "mod_inv (-5) 7 = 4"
75 value "mod_inv (-3) 7 = 2"
76 value "mod_inv (-4) 339 = 254"
78 (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
81 yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
82 definition mod_div :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
83 "mod_div a b m = ((nat a) * (mod_inv b m)) mod m"
85 definition "ASSERT_mod_div1 \<longleftrightarrow> mod_div 21 4 5 = 4" ML {* @{assert} @{code ASSERT_mod_div1} *}
86 definition "ASSERT_mod_div2 \<longleftrightarrow> mod_div 1 4 5 = 4" ML {* @{assert} @{code ASSERT_mod_div2} *}
87 definition "ASSERT_mod_div3 \<longleftrightarrow> mod_div 0 4 5 = 0" ML {* @{assert} @{code ASSERT_mod_div3} *}
89 value "mod_div 21 3 5 = 2" value "(21::int) mod 5 = (3 * 2) mod 5"
90 (* 21/3 = 7 mod 5 21 mod 5 = 6 mod 5
92 value "mod_div 22 3 5 = 4" value "(22::int) mod 5 = (3 * 4) mod 5"
93 (* 22/3 = ------- 22 mod 5 = 12 mod 5
95 value "mod_div 23 3 5 = 1" value "(23::int) mod 5 = (3 * 1) mod 5"
96 (* 23/3 = ------- 23 mod 5 = 3 mod 5
98 value "mod_div 24 3 5 = 3" value "(24::int) mod 5 = (3 * 3) mod 5"
99 (* 24/3 = ------- 24 mod 5 = 9 mod 5
101 value "mod_div 25 3 5 = 0" value "(25::int) mod 5 = (3 * 0) mod 5"
102 (* 25/3 = ------- 25 mod 5 = 0 mod 5
104 value "mod_div 21 4 5 = 4" value "(21::int) mod 5 = (4 * 4) mod 5"
105 value "mod_div 1 4 5 = 4" value "( 1::int) mod 5 = (4 * 4) mod 5"
106 value "mod_div 0 4 5 = 0" value "( 0::int) mod 5 = (0 * 4) mod 5"
108 (* root1 is just local to approx_root *)
109 function root1 :: "int \<Rightarrow> nat \<Rightarrow> nat" where
110 "root1 a n = (if (int (n * n)) < a then root1 a (n + 1) else n)"
112 termination sorry (*by (relation "measure (\<lambda>(a, n). nat (a - (int (n * n))))") auto*)
114 (* required just for one approximation:
115 approx_root :: nat \<Rightarrow> nat
118 yields r * r \<ge> a *)
119 definition approx_root :: "int \<Rightarrow> nat" where "approx_root a = root1 a 1"
121 (* chinese_remainder :: int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int
122 chinese_remainder (r1, m1) (r2, m2) = r
124 yields r = r1 mod m1 \<and> r = r2 mod m2 *)
125 definition chinese_remainder :: "int \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
126 "chinese_remainder r1 m1 r2 m2 =
127 ((nat (r1 mod (int m1))) +
128 (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv (int m1) m2))) mod (int m2))) *
131 value "chinese_remainder 17 9 3 4 = 35"
132 value "chinese_remainder 7 2 6 11 = 17"
134 subsection {* creation of lists of primes for efficiency *}
136 (* is_prime :: int list \<Rightarrow> int \<Rightarrow> bool
138 assumes max ps < n \<and> n \<le> (max ps)^2 \<and>
139 (* FIXME: really ^^^^^^^^^^^^^^^? *)
140 (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and> (\<forall>p. p < n \<and> Primes.prime p \<longrightarrow> List.member ps p)
141 yields Primes.prime n *)
142 fun is_prime :: "nat list \<Rightarrow> nat \<Rightarrow> bool" where
144 (if List.length ps > 0
146 if (n mod (List.hd ps)) = 0
148 else is_prime (List.tl ps) n
150 declare is_prime.simps [simp del] \<comment> \<open>make_primes, next_prime_not_dvd\<close>
152 value "is_prime [2, 3] 2 = False" \<comment> \<open>... precondition!\<close>
153 value "is_prime [2, 3] 3 = False" \<comment> \<open>... precondition!\<close>
154 value "is_prime [2, 3] 4 = False"
155 value "is_prime [2, 3] 5 = True"
156 value "is_prime [2, 3, 5] 5 = False" \<comment> \<open>... precondition!\<close>
157 value "is_prime [2, 3] 6 = False"
158 value "is_prime [2, 3] 7 = True"
159 value "is_prime [2, 3] 25 = True" \<comment> \<open>... because 5 not in list\<close>
161 (* make_primes is just local to primes_upto only:
162 make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list
163 make_primes ps last_p n = pps
164 assumes last_p = maxs ps \<and> (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
165 (\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)
166 yields n \<le> maxs pps \<and> (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p) \<and>
167 (\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
168 function make_primes :: "nat \<Rightarrow> nat \<Rightarrow> nat list \<Rightarrow> nat list" where
169 "make_primes last_p n ps =
170 (if n \<le> last ps then ps
171 else make_primes (last_p + 2) n
172 (if is_prime ps (last_p + 2) then ps @ [last_p + 2] else ps))"
173 by pat_completeness auto \<comment> \<open>simp del: is_prime.simps <-- declare\<close>
174 termination make_primes (*by lexicographic_order +PROOF primes? / size_change LOOPS*)
175 sorry \<comment> \<open>FIXME proof needs semantic properties of primes themselves\<close>
177 declare make_primes.simps [simp del] \<comment> \<open>next_prime_not_dvd\<close>
179 value "make_primes 3 3 [2, 3] = [2, 3]"
180 value "make_primes 3 4 [2, 3] = [2, 3, 5]"
181 value "make_primes 3 5 [2, 3] = [2, 3, 5]"
182 value "make_primes 3 6 [2, 3] = [2, 3, 5, 7]"
183 value "make_primes 3 7 [2, 3] = [2, 3, 5, 7]"
184 value "make_primes 3 8 [2, 3] = [2, 3, 5, 7, 11]"
185 value "make_primes 3 9 [2, 3] = [2, 3, 5, 7, 11]"
186 value "make_primes 3 10 [2, 3] = [2, 3, 5, 7, 11]"
187 value "make_primes 3 11 [2, 3] = [2, 3, 5, 7, 11]"
188 value "make_primes 3 12 [2, 3] = [2, 3, 5, 7, 11, 13]"
189 value "make_primes 3 13 [2, 3] = [2, 3, 5, 7, 11, 13]"
190 value "make_primes 3 14 [2, 3] = [2, 3, 5, 7, 11, 13, 17]"
191 value "make_primes 3 15 [2, 3] = [2, 3, 5, 7, 11, 13, 17]"
192 value "make_primes 3 16 [2, 3] = [2, 3, 5, 7, 11, 13, 17]"
193 value "make_primes 3 17 [2, 3] = [2, 3, 5, 7, 11, 13, 17]"
194 value "make_primes 3 18 [2, 3] = [2, 3, 5, 7, 11, 13, 17, 19]"
195 value "make_primes 3 19 [2, 3] = [2, 3, 5, 7, 11, 13, 17, 19]"
196 value "make_primes 7 4 [2, 3, 5, 7] = [2, 3, 5, 7]"
198 (* primes_upto :: nat \<Rightarrow> nat list
201 yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
202 n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
203 (\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
204 definition primes_upto :: "nat \<Rightarrow> nat list" where
205 "primes_upto n = (if n < 3 then [2] else make_primes 3 n [2, 3])"
207 value "primes_upto 0 = [2]"
208 value "primes_upto 1 = [2]"
209 value "primes_upto 2 = [2]"
210 value "primes_upto 3 = [2, 3]"
211 value "primes_upto 4 = [2, 3, 5]"
212 value "primes_upto 5 = [2, 3, 5]"
213 value "primes_upto 6 = [2, 3, 5, 7]"
214 value "primes_upto 7 = [2, 3, 5, 7]"
215 value "primes_upto 8 = [2, 3, 5, 7, 11]"
216 value "primes_upto 9 = [2, 3, 5, 7, 11]"
217 value "primes_upto 10 = [2, 3, 5, 7, 11]"
218 value "primes_upto 11 = [2, 3, 5, 7, 11]"
220 lemma primes_upto_0: "last (primes_upto n) > 0" (*see above*) sorry
221 lemma primes_upto_1: "last (primes_upto (Suc n)) > n" (*see above*)
222 apply (simp add: primes_upto_def)
223 apply (induct rule: make_primes.induct)
224 apply auto (*... same problems as with "make_primes" *)
226 lemma primes_upto_2: "last (primes_upto n) >= n" (*see above*) sorry
228 (* max's' is analogous to Integer.gcds; used ONLY in specifications, not in functions *)
229 definition maxs :: "nat list \<Rightarrow> nat" where "maxs ns = List.fold max ns (List.hd ns)"
230 value "maxs [5, 3, 7, 1, 2, 4] = 7"
232 (* find the next prime greater p not dividing the number n:
233 next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
234 n1 next_prime_not_dvd n2 = p
235 assumes True assumes "0 < q"
237 yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
238 (\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
239 function next_prime_not_dvd :: "nat \<Rightarrow> nat \<Rightarrow> nat" (infixl "next'_prime'_not'_dvd" 70) where
240 "n1 next_prime_not_dvd n2 =
242 ps = primes_upto (n1 + 1);
245 if n2 mod nxt \<noteq> 0
247 else nxt next_prime_not_dvd n2)"
248 by auto \<comment> \<open>simp del: is_prime.simps, make_primes.simps, primes_upto.simps < -- declare*\<close>
249 termination sorry (*next_prime_not_dvd: lexicographic_order +PROOF primes? / size_change: Failed*)
251 value "1 next_prime_not_dvd 15 = 2"
252 value "2 next_prime_not_dvd 15 = 7"
253 value "3 next_prime_not_dvd 15 = 7"
254 value "4 next_prime_not_dvd 15 = 7"
255 value "5 next_prime_not_dvd 15 = 7"
256 value "6 next_prime_not_dvd 15 = 7"
257 value "7 next_prime_not_dvd 15 =11"
259 subsection {* basic notions for univariate polynomials *}
261 (* not in List.thy, copy from library.ML *)
262 fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list" where
263 "nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"
264 value "nth_drop 0 [] = []"
265 value "nth_drop 0 [1, 2, 3::int] = [2, 3]"
266 value "nth_drop 2 [1, 2, 3::int] = [1, 2]"
267 value "nth_drop 77 [1, 2, 3::int] = [1, 2, 3]"
269 (* leading coefficient *)
270 definition lcoeff_up :: "unipoly \<Rightarrow> int" where "lcoeff_up p = (last o drop_tl_zeros) p"
272 value "lcoeff_up [3, 4, 5, 6] = 6"
273 value "lcoeff_up [3, 4, 5, 6, 0] = 6"
276 definition deg_up :: "unipoly \<Rightarrow> nat" where
277 "deg_up p = ((\<lambda>k. k - 1) o length o drop_tl_zeros) p"
278 (* FH wrong: (op - 1) o *)
280 value "degree (Coeff [1::int, 2, 3])"
281 value "deg_up [1, 2, 3]"
282 value "deg_up [3, 4, 5, 6] = 3"
283 value "deg_up [3, 4, 5, 6, 0] = 3"
284 value "deg_up [1, 0, 3, 0, 0] = 2"
286 (* norm is just local to normalise *)
287 fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where
288 "norm p nrm m lcp i =
290 then [int (mod_div (p ! i) lcp m)] @ nrm
291 else norm p ([int (mod_div (p ! i) lcp m)] @ nrm) m lcp (i - 1))"
292 (* normalise a unipoly such that lcoeff_up mod m = 1.
293 normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
294 normalise [p_0, .., p_k] m = [q_0, .., q_k]
296 yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i \<and> (p_k * t) mod m = 1 *)
297 fun normalise :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
303 if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))"
304 declare normalise.simps [simp del] \<comment> \<open>HENSEL_lifting_up\<close>
306 value "normalise [-18, -15, -20, 12, 20, -13, 2] 5 = [1, 0, 0, 1, 0, 1, 1]"
307 value "normalise [9, 6, 3] 10 = [3, 2, 1]"
309 subsection {* addition, multiplication, division *}
311 (* scalar multiplication *)
312 definition mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
313 "p %* m = List.map (op * m) p"
315 value "[5, 4, 7, 8, 1] %* 5 = [25, 20, 35, 40, 5]"
316 value "[5, 4, -7, 8, -1] %* 5 = [25, 20, -35, 40, -5]"
319 (*definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
320 CODEGEN CAUSES ERROR:
321 ML error (line 913 of "/home/neuper/devel/isac/codegen/gcd_univariate.sml"):
322 Type error in function application. Function: div2 : inta -> inta -> inta
323 Argument: (swapf, m) : (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd
324 Reason: Can't unify inta to (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd (Incompatible types)
325 THUS TYPE CONSTRAINED...
327 definition swapf1 :: "(int \<Rightarrow> int \<Rightarrow> int) \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where "swapf1 f a b = f b a"
328 definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%'/" 70) where
329 "p %/ m = map (swapf1 op div2 m) p"
331 value "[4, 3, 2, 5, 6] %/ 3 = [1, 1, 0, 1, 2]"
332 value "[4, 3, 2, 0] %/ 3 = [1, 1, 0, 0]"
334 (* not in List.thy, copy from library.ML *)
335 fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where
336 "map2 _ [] [] = []" |
337 "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys" |
338 "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
341 definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "%-%" 70) where
342 "p1 %-% p2 = map2 (op -) p1 p2"
344 value "[8, -7, 0, 1] %-% [-2, 2, 3, 0] = [10, -9, -3, 1]"
345 value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
347 function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where
348 "[d] %|% [p] \<longleftrightarrow> (\<bar>d\<bar> \<le> \<bar>p\<bar>) \<and> (p mod d = 0)" |
349 "ds %|% ps \<longleftrightarrow> (*a hint by FH*)
351 ds = drop_tl_zeros ds; ps = drop_tl_zeros ps;
352 d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
353 quot = (lcoeff_up ps) div2 (lcoeff_up d000);
354 rest = drop_tl_zeros (ps %-% (d000 %* quot))
356 rest = [] \<or> (quot \<noteq> 0 \<and> List.length ds \<le> List.length rest \<and> ds %|% rest))"
357 apply pat_completeness
359 done (* > 1 sec IMPROVED BY FLORIAN'S drop_tl_zeros AND declare simp del:
360 centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)
361 termination (*dvd_up: by lexicographic_order LOOPS +PROOF primes? / size_change LOOPS*)
362 using [[linarith_split_limit = 999]]
363 apply (relation "measure (\<lambda>(_, ps). length ps)") (*a hint by FH*)
367 value "[4] %|% [6] = False"
368 value "[8] %|% [16, 0] = True"
369 value "[3, 2] %|% [0, 0, 9, 12, 4] = True"
370 value "[8, 0] %|% [16] = True"
372 subsection {* normalisation and Landau-Mignotte bound *}
374 (* centr is just local to centr_up *)
375 definition centr :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where
376 "centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"
378 (* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
379 normalise [p_0, .., p_k] m = [q_0, .., q_k]
381 yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
382 (where |^ x ^| means round x up to the next greater integer) *)
383 definition centr_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
387 mid = if (int m) mod mi = 0 then mi else mi + 1
388 in map (centr m mid) p)"
390 value "centr_up [7, 3, 5, 8, 1, 3] 10 = [-3, 3, 5, -2, 1, 3]"
391 value "centr_up [1, 2, 3, 4, 5] 2 = [1, 0, 1, 2, 3]"
392 value "centr_up [1, 2, 3, 4, 5] 3 = [1, -1, 0, 1, 2]"
393 value "centr_up [1, 2, 3, 4, 5] 4 = [1, 2, -1, 0, 1]"
394 value "centr_up [1, 2, 3, 4, 5] 5 = [1, 2, 3, -1, 0]"
395 value "centr_up [1, 2, 3, 4, 5] 6 = [1, 2, 3, -2, -1]"
396 value "centr_up [1, 2, 3, 4, 5] 7 = [1, 2, 3, 4, -2]"
397 value "centr_up [1, 2, 3, 4, 5] 8 = [1, 2, 3, 4, -3]"
398 value "centr_up [1, 2, 3, 4, 5] 9 = [1, 2, 3, 4, 5]"
399 value "centr_up [1, 2, 3, 4, 5] 10 = [1, 2, 3, 4, 5]"
400 value "centr_up [1, 2, 3, 4, 5] 11 = [1, 2, 3, 4, 5]"
402 (*definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
403 CODEGEN CAUSES ERROR:
404 ML error (line 913 of "/home/neuper/devel/isac/codegen/gcd_univariate.sml"):
405 Type error in function application. Function: div2 : inta -> inta -> inta
406 Argument: (swapf, m) : (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd
407 Reason: Can't unify inta to (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd (Incompatible types)
408 THUS TYPE CONSTRAINED...
410 definition swapf2 :: "(int \<Rightarrow> nat \<Rightarrow> int) \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> int" where "swapf2 f a b = f b a"
412 (* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
413 sum_lmb [p_0, .., p_k] e = s
415 yields. p_0^e + p_1^e + ... + p_k^e *)
416 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where
417 "sum_lmb p e = List.fold ((op +) o (swapf2 power e)) p 0"
419 (*value "power" at outcommented Isabelle2015->17*)
421 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"
422 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"
423 value "sum_lmb [-1, 2, -3, 4, -5] 3 = -81"
424 value "sum_lmb [-1, 2, -3, 4, -5] 4 = 979"
425 value "sum_lmb [-1, 2, -3, 4, -5] 5 = -2313"
426 value "sum_lmb [-1, 2, -3, 4, -5] 6 = 20515"
428 (* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly => unipoly \<Rightarrow> int
429 LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
431 yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
432 min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)
433 definition LANDAU_MIGNOTTE_bound :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat" where
434 "LANDAU_MIGNOTTE_bound p1 p2 =
435 ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) *
436 (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))
437 (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"
439 value "LANDAU_MIGNOTTE_bound [1] [4, 5] = 1"
440 value "LANDAU_MIGNOTTE_bound [1, 2] [4, 5] = 2"
441 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5] = 2"
442 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4] = 1"
443 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5] = 2"
444 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5, 6] = 12"
446 value "LANDAU_MIGNOTTE_bound [-1] [4, 5] = 1"
447 value "LANDAU_MIGNOTTE_bound [-1, 2] [4, 5] = 2"
448 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5] = 2"
449 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4] = 1"
450 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5] = 2"
451 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"
453 subsection {* modulo calculations for polynomials *}
455 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)
456 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where
457 "([] pair []) = []" |
458 "([] pair ys) = []" | (*raise ListPair.UnequalLengths*)
459 "(xs pair []) = []" | (*raise ListPair.UnequalLengths*)
460 "((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
461 fun chinese_rem :: "nat \<times> nat \<Rightarrow> int \<times> int \<Rightarrow> int" where
462 "chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"
464 (* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
465 chinese_remainder_up (m1, m2) (p1, p2) = p
466 assume m1, m2 relatively prime
467 yields p1 = p mod m1 \<and> p2 = p mod m2 *)
468 fun chinese_remainder_up :: "nat \<times> nat \<Rightarrow> unipoly \<times> unipoly \<Rightarrow> unipoly" where
469 "chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"
471 value "chinese_remainder_up (5, 7) ([2, 2, 4, 3], [3, 2, 3, 5]) = [17, 2, 24, 33]"
473 (* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
474 mod_up [p1, p2, ..., pk] m = up
476 yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
477 definition mod' :: "nat \<Rightarrow> int \<Rightarrow> int" where "mod' m i = i mod (int m)"
478 definition mod_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" (infixl "mod'_up" 70) where
479 "p mod_up m = map (mod' m) p"
481 value "[5, 4, 7, 8, 1] mod_up 5 = [0, 4, 2, 3, 1]"
482 value "[5, 4,-7, 8,-1] mod_up 5 = [0, 4, 3, 3, 4]"
484 (* euclidean algoritm in Z_p[x/m].
485 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly
486 mod_up_gcd p1 p2 m = g
488 yields gcd (p1 mod m) (p2 mod m) = g *)
489 function mod_up_gcd :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
490 "mod_up_gcd p1 p2 m =
493 p2m = drop_tl_zeros (p2 mod_up m);
494 p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
495 quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
496 rest = drop_tl_zeros ((p1m %-% (p2n %* (int quot))) mod_up m)
498 if rest = [] then p2 else
499 if List.length rest < List.length p2
500 then mod_up_gcd p2 rest m
501 else mod_up_gcd rest p2 m)"
503 termination mod_up_gcd (*by lexicographic_order +PROOF primes? / size_change LOOPS*)
505 declare mod_up_gcd.simps [simp del] \<comment> \<open>HENSEL_lifting_up\<close>
507 value "mod_up_gcd [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] 7 = [2, 6, 0, 2, 6]"
508 value "mod_up_gcd [8, 28, 22, -11, -14, 1, 2] [2, 6, 0, 2, 6] 7 = [2, 6, 0, 2, 6]"
509 value "mod_up_gcd [20, 15, 8, 6] [8, -2, -2, 3] 2 = [0, 1]"
510 value "[20, 15, 8, 6] %|% [0, 1] = False"
511 value "[8, -2, -2, 3] %|% [0, 1] = False"
513 (* analogous to Integer.gcds
514 gcds :: int list \<Rightarrow> int
517 yields THE d. ((\<forall>n. List.member ns n \<longrightarrow> d dvd n) \<and>
518 (\<forall>d'. (\<forall>n. List.member ns n \<and> d' dvd n) \<longrightarrow> d'modp \<le> d)) *)
519 fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)" (*FH Gcd, set ?*)
521 value "gcds [6, 9, 12] = 3"
522 value "gcds [6, -9, 12] = 3"
523 value "gcds [8, 12, 16] = 4"
524 value "gcds [-8, 12, -16] = 4"
526 (* prim_poly :: unipoly \<Rightarrow> unipoly
529 yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
530 fun primitive_up :: "unipoly \<Rightarrow> unipoly" where
531 "primitive_up [c] = (if c = 0 then [0] else [1])" |
535 if d = 1 then p else p %/ d)"
537 value "primitive_up [12, 16, 32, 44] = [3, 4, 8, 11]"
538 value "primitive_up [4, 5, 12] = [4, 5, 12]"
539 value "primitive_up [0] = [0]"
540 value "primitive_up [6] = [1]"
542 subsection {* gcd_up, code from [1] p.93 *}
543 (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
544 try_new_prime_up a b d M P g p = new_g
545 assumes d = gcd (lcoeff_up a, lcoeff_up b) \<and>
546 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
547 a is primitive \<and> b is primitive
548 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
550 argumentns "a b d M P g p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
551 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
553 "try_new_prime_up a b d M P g p =
554 (if P > M then g else
555 let p = p next_prime_not_dvd d;
556 g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
561 if deg_up g_p < deg_up g
563 if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p
565 if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p else
568 g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
569 in try_new_prime_up a b d M P g p)"
570 by pat_completeness auto \<comment> \<open>simp del: centr_up_def normalise.simps mod_up_gcd.simps\<close>
571 termination try_new_prime_up (*by lexicographic_order +PROOF primes? / by size_change LOOPS*)
574 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
575 HENSEL_lifting_up p1 p2 d M p = gcd
576 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
577 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and>
578 p1 is primitive \<and> p2 is primitive
579 yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
581 argumentns "a b d M p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
582 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where
583 "HENSEL_lifting_up a b d M p =
585 p = p next_prime_not_dvd d;
586 g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p (*see above*)
588 if deg_up g_p = 0 then [1] else
590 g = primitive_up (try_new_prime_up a b d M p g_p p)
592 if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"
593 by pat_completeness auto \<comment> \<open>simp del: centr_up_def normalise.simps mod_up_gcd.simps\<close>
594 termination HENSEL_lifting_up (*by lexicographic_order LOOPS +PROOF primes? / by size_change LOOPS*)
597 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
599 assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
600 a is primitive \<and> b is primitive
601 yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
602 function gcd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" where
604 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>
607 HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)"
608 by pat_completeness auto \<comment> \<open>simp del: lcoeff_up.simps ?+ others?\<close>
609 termination by lexicographic_order (*works*)
611 ML {*"----------- fun gcd_up ---------------------------------";*}
612 value "gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
613 definition "ASSERT_gcd_up1 \<longleftrightarrow>
614 gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
615 ML {* @{assert} @{code ASSERT_gcd_up1} *}
617 (* gcd (-1 + x^2) (x + x^2) = (1 + x) ...*)
618 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
619 definition "ASSERT_gcd_up2 \<longleftrightarrow> gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
620 ML {* @{assert} @{code ASSERT_gcd_up2} *}
624 declare [[simp_trace_depth_limit = 99]]
625 declare [[simp_trace = true]]
627 using [[simp_trace_depth_limit = 99]]
628 using [[simp_trace = true]]