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