11 test/../gcd_poly.sml documents the changes towards GCD_Poly.thy; |
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. |
12 the style of GCD_Poly.thy has been adapted to the function package. |
13 *} |
13 *} |
14 |
14 |
15 section {* Isabelle's predefined polynomials *} |
15 section {* Isabelle's predefined polynomials *} |
16 --"TODO" |
16 \<comment> \<open>TODO\<close> |
17 |
17 |
18 section {* gcd for univariate polynomials *} |
18 section {* gcd for univariate polynomials *} |
19 |
19 |
20 type_synonym unipoly = "int list" (*TODO: compare Polynomial.thy*) |
20 type_synonym unipoly = "int list" (*TODO: compare Polynomial.thy*) |
21 value "[0, 1, 2, 3, 4, 5] :: unipoly" |
21 value "[0, 1, 2, 3, 4, 5] :: unipoly" |
145 then |
145 then |
146 if (n mod (List.hd ps)) = 0 |
146 if (n mod (List.hd ps)) = 0 |
147 then False |
147 then False |
148 else is_prime (List.tl ps) n |
148 else is_prime (List.tl ps) n |
149 else True)" |
149 else True)" |
150 declare is_prime.simps [simp del] -- "make_primes, next_prime_not_dvd" |
150 declare is_prime.simps [simp del] \<comment> \<open>make_primes, next_prime_not_dvd\<close> |
151 |
151 |
152 value "is_prime [2, 3] 2 = False" --"... precondition!" |
152 value "is_prime [2, 3] 2 = False" \<comment> \<open>... precondition!\<close> |
153 value "is_prime [2, 3] 3 = False" --"... precondition!" |
153 value "is_prime [2, 3] 3 = False" \<comment> \<open>... precondition!\<close> |
154 value "is_prime [2, 3] 4 = False" |
154 value "is_prime [2, 3] 4 = False" |
155 value "is_prime [2, 3] 5 = True" |
155 value "is_prime [2, 3] 5 = True" |
156 value "is_prime [2, 3, 5] 5 = False" --"... precondition!" |
156 value "is_prime [2, 3, 5] 5 = False" \<comment> \<open>... precondition!\<close> |
157 value "is_prime [2, 3] 6 = False" |
157 value "is_prime [2, 3] 6 = False" |
158 value "is_prime [2, 3] 7 = True" |
158 value "is_prime [2, 3] 7 = True" |
159 value "is_prime [2, 3] 25 = True" -- "... because 5 not in list" |
159 value "is_prime [2, 3] 25 = True" \<comment> \<open>... because 5 not in list\<close> |
160 |
160 |
161 (* make_primes is just local to primes_upto only: |
161 (* make_primes is just local to primes_upto only: |
162 make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list |
162 make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list |
163 make_primes ps last_p n = pps |
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> |
164 assumes last_p = maxs ps \<and> (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and> |
168 function make_primes :: "nat \<Rightarrow> nat \<Rightarrow> nat list \<Rightarrow> nat list" where |
168 function make_primes :: "nat \<Rightarrow> nat \<Rightarrow> nat list \<Rightarrow> nat list" where |
169 "make_primes last_p n ps = |
169 "make_primes last_p n ps = |
170 (if n \<le> last ps then ps |
170 (if n \<le> last ps then ps |
171 else make_primes (last_p + 2) n |
171 else make_primes (last_p + 2) n |
172 (if is_prime ps (last_p + 2) then ps @ [last_p + 2] else ps))" |
172 (if is_prime ps (last_p + 2) then ps @ [last_p + 2] else ps))" |
173 by pat_completeness auto --"simp del: is_prime.simps <-- declare" |
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*) |
174 termination make_primes (*by lexicographic_order +PROOF primes? / size_change LOOPS*) |
175 sorry -- {* FIXME proof needs semantic properties of primes themselves *} |
175 sorry \<comment> \<open>FIXME proof needs semantic properties of primes themselves\<close> |
176 |
176 |
177 declare make_primes.simps [simp del] -- "next_prime_not_dvd" |
177 declare make_primes.simps [simp del] \<comment> \<open>next_prime_not_dvd\<close> |
178 |
178 |
179 value "make_primes 3 3 [2, 3] = [2, 3]" |
179 value "make_primes 3 3 [2, 3] = [2, 3]" |
180 value "make_primes 3 4 [2, 3] = [2, 3, 5]" |
180 value "make_primes 3 4 [2, 3] = [2, 3, 5]" |
181 value "make_primes 3 5 [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]" |
182 value "make_primes 3 6 [2, 3] = [2, 3, 5, 7]" |
243 nxt = last ps |
243 nxt = last ps |
244 in |
244 in |
245 if n2 mod nxt \<noteq> 0 |
245 if n2 mod nxt \<noteq> 0 |
246 then nxt |
246 then nxt |
247 else nxt next_prime_not_dvd n2)" |
247 else nxt next_prime_not_dvd n2)" |
248 by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare" |
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*) |
249 termination sorry (*next_prime_not_dvd: lexicographic_order +PROOF primes? / size_change: Failed*) |
250 |
250 |
251 value "1 next_prime_not_dvd 15 = 2" |
251 value "1 next_prime_not_dvd 15 = 2" |
252 value "2 next_prime_not_dvd 15 = 7" |
252 value "2 next_prime_not_dvd 15 = 7" |
253 value "3 next_prime_not_dvd 15 = 7" |
253 value "3 next_prime_not_dvd 15 = 7" |
299 (let |
299 (let |
300 p = drop_tl_zeros p; |
300 p = drop_tl_zeros p; |
301 lcp = lcoeff_up p |
301 lcp = lcoeff_up p |
302 in |
302 in |
303 if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))" |
303 if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))" |
304 declare normalise.simps [simp del] --"HENSEL_lifting_up" |
304 declare normalise.simps [simp del] \<comment> \<open>HENSEL_lifting_up\<close> |
305 |
305 |
306 value "normalise [-18, -15, -20, 12, 20, -13, 2] 5 = [1, 0, 0, 1, 0, 1, 1]" |
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]" |
307 value "normalise [9, 6, 3] 10 = [3, 2, 1]" |
308 |
308 |
309 subsection {* addition, multiplication, division *} |
309 subsection {* addition, multiplication, division *} |
500 then mod_up_gcd p2 rest m |
500 then mod_up_gcd p2 rest m |
501 else mod_up_gcd rest p2 m)" |
501 else mod_up_gcd rest p2 m)" |
502 by auto |
502 by auto |
503 termination mod_up_gcd (*by lexicographic_order +PROOF primes? / size_change LOOPS*) |
503 termination mod_up_gcd (*by lexicographic_order +PROOF primes? / size_change LOOPS*) |
504 sorry |
504 sorry |
505 declare mod_up_gcd.simps [simp del] --"HENSEL_lifting_up" |
505 declare mod_up_gcd.simps [simp del] \<comment> \<open>HENSEL_lifting_up\<close> |
506 |
506 |
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]" |
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]" |
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]" |
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" |
510 value "[20, 15, 8, 6] %|% [0, 1] = False" |
565 if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p else |
565 if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p else |
566 let |
566 let |
567 P = P * p; |
567 P = P * p; |
568 g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P |
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)" |
569 in try_new_prime_up a b d M P g p)" |
570 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps" |
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*) |
571 termination try_new_prime_up (*by lexicographic_order +PROOF primes? / by size_change LOOPS*) |
572 sorry |
572 sorry |
573 |
573 |
574 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly |
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 |
575 HENSEL_lifting_up p1 p2 d M p = gcd |
588 if deg_up g_p = 0 then [1] else |
588 if deg_up g_p = 0 then [1] else |
589 (let |
589 (let |
590 g = primitive_up (try_new_prime_up a b d M p g_p p) |
590 g = primitive_up (try_new_prime_up a b d M p g_p p) |
591 in |
591 in |
592 if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))" |
592 if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))" |
593 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps" |
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*) |
594 termination HENSEL_lifting_up (*by lexicographic_order LOOPS +PROOF primes? / by size_change LOOPS*) |
595 sorry |
595 sorry |
596 |
596 |
597 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly |
597 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly |
598 gcd_up a b = c |
598 gcd_up a b = c |
603 "gcd_up a b = |
603 "gcd_up a b = |
604 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar> |
604 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar> |
605 in |
605 in |
606 if a = b then a else |
606 if a = b then a else |
607 HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)" |
607 HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)" |
608 by pat_completeness auto --"simp del: lcoeff_up.simps ?+ others?" |
608 by pat_completeness auto \<comment> \<open>simp del: lcoeff_up.simps ?+ others?\<close> |
609 termination by lexicographic_order (*works*) |
609 termination by lexicographic_order (*works*) |
610 |
610 |
611 ML {*"----------- fun gcd_up ---------------------------------";*} |
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]" |
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> |
613 definition "ASSERT_gcd_up1 \<longleftrightarrow> |