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