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_OLD.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_old.sml;
11 test/../gcd_poly_ml.sml documents the changes towards GCD_Poly_ML.thy;
12 the style of GCD_Poly_OLD.thy has been adapted to the function package.
15 section \<open>Isabelle's predefined polynomials\<close>
16 \<comment> \<open>TODO\<close>
18 section \<open>gcd for univariate polynomials\<close>
20 type_synonym unipoly = "int list" (*TODO: compare Polynomial.thy*)
21 value "[0, 1, 2, 3, 4, 5] :: unipoly"
23 subsection \<open>auxiliary functions\<close>
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 \<open>modulo calculations for integers\<close>
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 \<open>@{assert} @{code ASSERT_mod_div1}\<close>
86 definition "ASSERT_mod_div2 \<longleftrightarrow> mod_div 1 4 5 = 4" ML \<open>@{assert} @{code ASSERT_mod_div2}\<close>
87 definition "ASSERT_mod_div3 \<longleftrightarrow> mod_div 0 4 5 = 0" ML \<open>@{assert} @{code ASSERT_mod_div3}\<close>
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 \<open>creation of lists of primes for efficiency\<close>
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 "next_prime_not_dvd n1 n2 =
242 ps = primes_upto (n1 + 1);
245 if n2 mod nxt \<noteq> 0
247 else next_prime_not_dvd nxt 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 "next_prime_not_dvd 1 15 = 2"
252 value "next_prime_not_dvd 2 15 = 7"
253 value "next_prime_not_dvd 3 15 = 7"
254 value "next_prime_not_dvd 4 15 = 7"
255 value "next_prime_not_dvd 5 15 = 7"
256 value "next_prime_not_dvd 6 15 = 7"
257 value "next_prime_not_dvd 7 15 =11"
259 subsection \<open>basic notions for univariate polynomials\<close>
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 \<open>addition, multiplication, division\<close>
311 (* scalar multiplication *)
312 definition mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
313 "p %* m = List.map (\<lambda>i. m * i) 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 div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%'/" 70) where
328 "p %/ m = map (\<lambda>i. i div2 m) p"
330 value "[4, 3, 2, 5, 6] %/ 3 = [1, 1, 0, 1, 2]"
331 value "[4, 3, 2, 0] %/ 3 = [1, 1, 0, 0]"
333 (* not in List.thy, copy from library.ML *)
334 fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where
335 "map2 _ [] [] = []" |
336 "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys" |
337 "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
340 definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "%-%" 70) where
341 "p1 %-% p2 = map2 (\<lambda>i j. i - j) p1 p2"
343 value "[8, -7, 0, 1] %-% [-2, 2, 3, 0] = [10, -9, -3, 1]"
344 value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
346 function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where
347 "[d] %|% [p] \<longleftrightarrow> (\<bar>d\<bar> \<le> \<bar>p\<bar>) \<and> (p mod d = 0)" |
348 "ds %|% ps \<longleftrightarrow> \<comment> \<open>a hint by FH\<close>
350 ds = drop_tl_zeros ds; ps = drop_tl_zeros ps;
351 d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
352 quot = (lcoeff_up ps) div2 (lcoeff_up d000);
353 rest = drop_tl_zeros (ps %-% (d000 %* quot))
355 rest = [] \<or> (quot \<noteq> 0 \<and> List.length ds \<le> List.length rest \<and> ds %|% rest))"
356 apply pat_completeness
358 done (* > 1 sec IMPROVED BY FLORIAN'S drop_tl_zeros AND declare simp del:
359 centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)
360 termination (*dvd_up: by lexicographic_order LOOPS +PROOF primes? / size_change LOOPS*)
361 using [[linarith_split_limit = 999]]
362 apply (relation "measure (\<lambda>(_, ps). length ps)") (*a hint by FH*)
366 value "[4] %|% [6] = False"
367 value "[8] %|% [16, 0] = True"
368 value "[3, 2] %|% [0, 0, 9, 12, 4] = True"
369 value "[8, 0] %|% [16] = True"
371 subsection \<open>normalisation and Landau-Mignotte bound\<close>
373 (* centr is just local to centr_up *)
374 definition centr :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where
375 "centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"
377 (* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
378 normalise [p_0, .., p_k] m = [q_0, .., q_k]
380 yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
381 (where |^ x ^| means round x up to the next greater integer) *)
382 definition centr_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
386 mid = if (int m) mod mi = 0 then mi else mi + 1
387 in map (centr m mid) p)"
389 value "centr_up [7, 3, 5, 8, 1, 3] 10 = [-3, 3, 5, -2, 1, 3]"
390 value "centr_up [1, 2, 3, 4, 5] 2 = [1, 0, 1, 2, 3]"
391 value "centr_up [1, 2, 3, 4, 5] 3 = [1, -1, 0, 1, 2]"
392 value "centr_up [1, 2, 3, 4, 5] 4 = [1, 2, -1, 0, 1]"
393 value "centr_up [1, 2, 3, 4, 5] 5 = [1, 2, 3, -1, 0]"
394 value "centr_up [1, 2, 3, 4, 5] 6 = [1, 2, 3, -2, -1]"
395 value "centr_up [1, 2, 3, 4, 5] 7 = [1, 2, 3, 4, -2]"
396 value "centr_up [1, 2, 3, 4, 5] 8 = [1, 2, 3, 4, -3]"
397 value "centr_up [1, 2, 3, 4, 5] 9 = [1, 2, 3, 4, 5]"
398 value "centr_up [1, 2, 3, 4, 5] 10 = [1, 2, 3, 4, 5]"
399 value "centr_up [1, 2, 3, 4, 5] 11 = [1, 2, 3, 4, 5]"
401 (* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
402 sum_lmb [p_0, .., p_k] e = s
404 yields. p_0^e + p_1^e + ... + p_k^e *)
405 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where
406 "sum_lmb p e = List.fold ((\<lambda>i. (+) i) o (\<lambda>i. i ^ e)) p 0"
408 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"
409 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"
410 value "sum_lmb [-1, 2, -3, 4, -5] 3 = -81"
411 value "sum_lmb [-1, 2, -3, 4, -5] 4 = 979"
412 value "sum_lmb [-1, 2, -3, 4, -5] 5 = -2313"
413 value "sum_lmb [-1, 2, -3, 4, -5] 6 = 20515"
415 (* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly => unipoly \<Rightarrow> int
416 LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
418 yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
419 min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)
420 definition LANDAU_MIGNOTTE_bound :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat" where
421 "LANDAU_MIGNOTTE_bound p1 p2 =
422 ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) *
423 (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))
424 (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"
426 value "LANDAU_MIGNOTTE_bound [1] [4, 5] = 1"
427 value "LANDAU_MIGNOTTE_bound [1, 2] [4, 5] = 2"
428 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5] = 2"
429 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4] = 1"
430 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5] = 2"
431 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5, 6] = 12"
433 value "LANDAU_MIGNOTTE_bound [-1] [4, 5] = 1"
434 value "LANDAU_MIGNOTTE_bound [-1, 2] [4, 5] = 2"
435 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5] = 2"
436 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4] = 1"
437 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5] = 2"
438 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"
440 subsection \<open>modulo calculations for polynomials\<close>
442 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)
443 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where
444 "([] pair []) = []" |
445 "([] pair ys) = []" | (*raise ListPair.UnequalLengths*)
446 "(xs pair []) = []" | (*raise ListPair.UnequalLengths*)
447 "((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
448 fun chinese_rem :: "nat \<times> nat \<Rightarrow> int \<times> int \<Rightarrow> int" where
449 "chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"
451 (* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
452 chinese_remainder_up (m1, m2) (p1, p2) = p
453 assume m1, m2 relatively prime
454 yields p1 = p mod m1 \<and> p2 = p mod m2 *)
455 fun chinese_remainder_up :: "nat \<times> nat \<Rightarrow> unipoly \<times> unipoly \<Rightarrow> unipoly" where
456 "chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"
458 value "chinese_remainder_up (5, 7) ([2, 2, 4, 3], [3, 2, 3, 5]) = [17, 2, 24, 33]"
460 (* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
461 mod_up [p1, p2, ..., pk] m = up
463 yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
464 definition mod' :: "nat \<Rightarrow> int \<Rightarrow> int" where "mod' m i = i mod (int m)"
465 definition mod_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" (*infixl "mod_up" 70*) where
466 "mod_up p m = map (mod' m) p"
468 value "mod_up [5, 4, 7, 8, 1] 5 = [0, 4, 2, 3, 1]"
469 value "mod_up [5, 4,-7, 8,-1] 5 = [0, 4, 3, 3, 4]"
471 (* euclidean algoritm in Z_p[x/m].
472 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly
473 mod_up_gcd p1 p2 m = g
475 yields gcd (p1 mod m) (p2 mod m) = g *)
476 function mod_up_gcd :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
477 "mod_up_gcd p1 p2 m =
480 p2m = drop_tl_zeros (mod_up p2 m);
481 p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
482 quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
483 rest = drop_tl_zeros (mod_up (p1m %-% (p2n %* (int quot))) m)
485 if rest = [] then p2 else
486 if List.length rest < List.length p2
487 then mod_up_gcd p2 rest m
488 else mod_up_gcd rest p2 m)"
490 termination mod_up_gcd (*by lexicographic_order +PROOF primes? / size_change LOOPS*)
492 declare mod_up_gcd.simps [simp del] \<comment> \<open>HENSEL_lifting_up\<close>
494 value "mod_up_gcd [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] 7 = [2, 6, 0, 2, 6]"
495 value "mod_up_gcd [8, 28, 22, -11, -14, 1, 2] [2, 6, 0, 2, 6] 7 = [2, 6, 0, 2, 6]"
496 value "mod_up_gcd [20, 15, 8, 6] [8, -2, -2, 3] 2 = [0, 1]"
497 value "[20, 15, 8, 6] %|% [0, 1] = False"
498 value "[8, -2, -2, 3] %|% [0, 1] = False"
500 (* analogous to Integer.gcds
501 gcds :: int list \<Rightarrow> int
504 yields THE d. ((\<forall>n. List.member ns n \<longrightarrow> d dvd n) \<and>
505 (\<forall>d'. (\<forall>n. List.member ns n \<and> d' dvd n) \<longrightarrow> d'modp \<le> d)) *)
506 fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)" (*FH Gcd, set ?*)
508 value "gcds [6, 9, 12] = 3"
509 value "gcds [6, -9, 12] = 3"
510 value "gcds [8, 12, 16] = 4"
511 value "gcds [-8, 12, -16] = 4"
513 (* prim_poly :: unipoly \<Rightarrow> unipoly
516 yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
517 fun primitive_up :: "unipoly \<Rightarrow> unipoly" where
518 "primitive_up [c] = (if c = 0 then [0] else [1])" |
522 if d = 1 then p else p %/ d)"
524 value "primitive_up [12, 16, 32, 44] = [3, 4, 8, 11]"
525 value "primitive_up [4, 5, 12] = [4, 5, 12]"
526 value "primitive_up [0] = [0]"
527 value "primitive_up [6] = [1]"
529 subsection \<open>gcd_up, code from [1] p.93\<close>
530 (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
531 try_new_prime_up a b d M P g p = new_g
532 assumes d = gcd (lcoeff_up a, lcoeff_up b) \<and>
533 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
534 a is primitive \<and> b is primitive
535 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
537 argumentns "a b d M P g p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
538 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
540 "try_new_prime_up a b d M P g p =
541 (if P > M then g else
542 let p = next_prime_not_dvd p d;
543 g_p = centr_up (mod_up ( (normalise (mod_up_gcd a b p) p)
544 %* (int (d mod p)))p)
547 if deg_up g_p < deg_up g
549 if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p
551 if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p else
554 g = centr_up (mod_up (chinese_remainder_up (P, p) (g, g_p)) P) P
555 in try_new_prime_up a b d M P g p)"
556 by pat_completeness auto \<comment> \<open>simp del: centr_up_def normalise.simps mod_up_gcd.simps\<close>
557 termination try_new_prime_up (*by lexicographic_order +PROOF primes? / by size_change LOOPS*)
560 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
561 HENSEL_lifting_up p1 p2 d M p = gcd
562 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
563 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and>
564 p1 is primitive \<and> p2 is primitive
565 yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
567 argumentns "a b d M p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
568 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where
569 "HENSEL_lifting_up a b d M p =
571 p = next_prime_not_dvd p d;
572 g_p = centr_up (mod_up ((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) p) p \<comment> \<open>see above\<close>
574 if deg_up g_p = 0 then [1] else
576 g = primitive_up (try_new_prime_up a b d M p g_p p)
578 if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"
579 by pat_completeness auto \<comment> \<open>simp del: centr_up_def normalise.simps mod_up_gcd.simps\<close>
580 termination HENSEL_lifting_up (*by lexicographic_order LOOPS +PROOF primes? / by size_change LOOPS*)
583 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
585 assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
586 a is primitive \<and> b is primitive
587 yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
588 function gcd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" where
590 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>
593 HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)"
594 by pat_completeness auto \<comment> \<open>simp del: lcoeff_up.simps ?+ others?\<close>
595 termination by lexicographic_order (*works*)
597 ML \<open>"----------- fun gcd_up ---------------------------------";\<close>
598 value "gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
599 definition "ASSERT_gcd_up1 \<longleftrightarrow>
600 gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
601 ML \<open>@{assert} @{code ASSERT_gcd_up1}\<close>
603 (* gcd (-1 + x^2) (x + x^2) = (1 + x) ...*)
604 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
605 definition "ASSERT_gcd_up2 \<longleftrightarrow> gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
606 ML \<open>@{assert} @{code ASSERT_gcd_up2}\<close>
610 declare [[simp_trace_depth_limit = 99]]
611 declare [[simp_trace = true]]
613 using [[simp_trace_depth_limit = 99]]
614 using [[simp_trace = true]]