meindl_diana@48797
|
1 |
(* rationals, i.e. fractions of multivariate polynomials over the real field
|
meindl_diana@48797
|
2 |
author: Diana Meindl
|
meindl_diana@48797
|
3 |
(c) Diana Meindl
|
meindl_diana@48797
|
4 |
Use is subject to license terms.
|
meindl_diana@48797
|
5 |
1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
|
meindl_diana@48797
|
6 |
10 20 30 40 50 60 70 80 90 100
|
meindl_diana@48797
|
7 |
*)
|
meindl_diana@48797
|
8 |
|
meindl_diana@48797
|
9 |
theory GCD_Poly imports Complex_Main begin
|
meindl_diana@48797
|
10 |
|
meindl_diana@48797
|
11 |
text {*
|
meindl_diana@48797
|
12 |
Below we reference
|
neuper@48852
|
13 |
F. Winkler, Polynomial algorithyms in computer algebra. Springer 1996.
|
neuper@48831
|
14 |
by page 93 and 95. This is the final version of the Isabelle theory,
|
neuper@48831
|
15 |
which conforms with Isabelle's coding standards and with requirements
|
neuper@48831
|
16 |
in terms of logics appropriate for a theorem prover.
|
neuper@48831
|
17 |
*}
|
neuper@48831
|
18 |
|
neuper@48831
|
19 |
section {* Predicates for the specifications *}
|
neuper@48831
|
20 |
text {*
|
neuper@48831
|
21 |
This thesis provides program code to be integrated into the Isabelle
|
neuper@48831
|
22 |
distribution. Thus the code will be verified against specifications, and
|
neuper@48831
|
23 |
consequently the code, as \emph{explicit} specification of how to calculate
|
neuper@48831
|
24 |
results, is accompanied by \emph{implicit} specifications describing the
|
neuper@48831
|
25 |
properties of the program code.
|
neuper@48831
|
26 |
|
neuper@48831
|
27 |
For describing the properties of code the format of Haskell will be used.
|
neuper@48831
|
28 |
Haskell \cite{TODO} is front-of-the-wave in research on functional
|
neuper@48831
|
29 |
programming and as such advances beyond SML. The advancement relevant for
|
neuper@48831
|
30 |
this thesis is the concise combination of explicit and implicit specification
|
neuper@48831
|
31 |
of a function. The final result of the code developed within this thesis is
|
neuper@48831
|
32 |
described as follows in Haskell format:
|
neuper@48831
|
33 |
|
neuper@48831
|
34 |
\begin{verbatim} % TODO: find respective Isabelle/LaTeX features
|
neuper@48831
|
35 |
gcd_poly :: int poly \<Rightarrow> int poly \<Rightarrow> int poly
|
neuper@48831
|
36 |
gcd_poly p1 p2 = d
|
neuper@48831
|
37 |
assumes and p1 \<noteq> [:0:] and p2 \<noteq> [:0:]
|
neuper@48851
|
38 |
yields d dvd p1 \<and> d dvd p2 \<and> \<forall>d'. (d' dvd p1 \<and> d' dvd p2) \<longrightarrow> d' \<le> d
|
neuper@48831
|
39 |
\end{verbatim
|
neuper@48831
|
40 |
|
neuper@48831
|
41 |
The above specification combines Haskell format with the term language of
|
neuper@48831
|
42 |
Isabelle\HOL. The latter provides a large and ever growing library of
|
neuper@48831
|
43 |
mechanised mathematics knowledge, which this thesis re-uses in order to
|
neuper@48831
|
44 |
prepare for integration of the code.
|
neuper@48831
|
45 |
|
neuper@48831
|
46 |
In the sequel there is a brief account of the notions required for the
|
neuper@48831
|
47 |
specification of the functions implemented within this thesis. *}
|
neuper@48831
|
48 |
|
neuper@48831
|
49 |
subsection {* The Isabelle types required *}
|
neuper@48831
|
50 |
text {*
|
neuper@48851
|
51 |
bool, nat, int, option, THE, list, 'a poly
|
neuper@48831
|
52 |
*}
|
neuper@48851
|
53 |
|
neuper@48831
|
54 |
subsection {* The connectives required *}
|
neuper@48831
|
55 |
text {*
|
neuper@48831
|
56 |
logic ~~/doc/tutorial.pdf p.203
|
neuper@48851
|
57 |
arithmetic < \<le> + - * ^ div \<bar>_\<bar>
|
neuper@48831
|
58 |
*}
|
neuper@48831
|
59 |
|
neuper@48831
|
60 |
subsection {* The Isabelle functions required *}
|
neuper@48831
|
61 |
text {*
|
neuper@48831
|
62 |
number theory
|
neuper@48831
|
63 |
dvd, mod, gcd
|
neuper@48831
|
64 |
Primes.prime
|
neuper@48831
|
65 |
Primes.next_prime_bound
|
neuper@48831
|
66 |
|
neuper@48831
|
67 |
lists
|
neuper@48831
|
68 |
List.hd, List.tl
|
neuper@48831
|
69 |
List.length, List.member
|
neuper@48831
|
70 |
List.fold, List.map
|
neuper@48831
|
71 |
List.nth, List.take, List.drop
|
neuper@48831
|
72 |
List.replicate
|
neuper@48831
|
73 |
|
neuper@48831
|
74 |
others
|
neuper@48831
|
75 |
Fact.fact
|
meindl_diana@48797
|
76 |
*}
|
meindl_diana@48797
|
77 |
|
neuper@48813
|
78 |
section {* gcd for univariate polynomials *}
|
meindl_diana@48797
|
79 |
ML {*
|
meindl_diana@48799
|
80 |
type unipoly = int list;
|
meindl_diana@48799
|
81 |
val a = [~18, ~15, ~20, 12, 20, ~13, 2];
|
meindl_diana@48799
|
82 |
val b = [8, 28, 22, ~11, ~14, 1, 2];
|
meindl_diana@48797
|
83 |
*}
|
meindl_diana@48797
|
84 |
|
neuper@48857
|
85 |
subsection {* auxiliary functions *}
|
neuper@48857
|
86 |
ML {*
|
neuper@48857
|
87 |
(* a variant for div:
|
neuper@48857
|
88 |
5 div 2 = 2; ~5 div 2 = ~3; BUT WEE NEED ~5 div2 2 = ~2; *)
|
neuper@48813
|
89 |
infix div2
|
neuper@48813
|
90 |
fun a div2 b =
|
neuper@48851
|
91 |
if a div b < 0 then (abs a) div (abs b) * ~1 else a div b
|
neuper@48857
|
92 |
|
neuper@48866
|
93 |
(* why has "last" been removed since 2002, although "last" is in List.thy ? *)
|
neuper@48866
|
94 |
fun last ls = (hd o rev) ls
|
neuper@48857
|
95 |
|
neuper@48857
|
96 |
(* drop tail elements equal 0 *)
|
neuper@48857
|
97 |
fun drop_hd_zeros (0 :: ps) = drop_hd_zeros ps
|
neuper@48857
|
98 |
| drop_hd_zeros p = p
|
neuper@48857
|
99 |
fun drop_tl_zeros p = (rev o drop_hd_zeros o rev) p;
|
meindl_diana@48797
|
100 |
*}
|
meindl_diana@48797
|
101 |
|
meindl_diana@48797
|
102 |
subsection {* modulo calculations for integers *}
|
meindl_diana@48797
|
103 |
ML {*
|
neuper@48836
|
104 |
(* mod_inv :: int \<Rightarrow> nat
|
neuper@48817
|
105 |
mod_inv r m = s
|
neuper@48836
|
106 |
assumes True
|
neuper@48836
|
107 |
yields r * s mod m = 1 *)
|
meindl_diana@48797
|
108 |
fun mod_inv r m =
|
meindl_diana@48797
|
109 |
let
|
meindl_diana@48797
|
110 |
fun modi (r, rold, m, rinv) =
|
neuper@48836
|
111 |
if m <= rinv then 0 else
|
meindl_diana@48797
|
112 |
if r mod m = 1
|
meindl_diana@48797
|
113 |
then rinv
|
meindl_diana@48797
|
114 |
else modi (rold * (rinv + 1), rold, m, rinv + 1)
|
neuper@48861
|
115 |
in modi (r, r, m, 1) end
|
neuper@48862
|
116 |
|
neuper@48865
|
117 |
(* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
|
neuper@48817
|
118 |
mod_div a b m = c
|
neuper@48836
|
119 |
assumes True
|
neuper@48865
|
120 |
yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
|
neuper@48862
|
121 |
fun mod_div a b m = a * (mod_inv b m) mod m
|
meindl_diana@48797
|
122 |
|
neuper@48836
|
123 |
(* required just for one approximation:
|
neuper@48836
|
124 |
approx_root :: nat \<Rightarrow> nat
|
neuper@48817
|
125 |
approx_root a = r
|
neuper@48836
|
126 |
assumes True
|
neuper@48836
|
127 |
yields r * r \<ge> a *)
|
neuper@48813
|
128 |
fun approx_root a =
|
neuper@48836
|
129 |
let
|
neuper@48836
|
130 |
fun root1 a n = if n * n < a then root1 a (n + 1) else n
|
neuper@48836
|
131 |
in root1 a 1 end
|
meindl_diana@48797
|
132 |
|
neuper@48862
|
133 |
(* chinese_remainder :: int * nat \<Rightarrow> int * nat \<Rightarrow> nat
|
neuper@48862
|
134 |
chinese_remainder (r1, m1) (r2, m2) = r
|
neuper@48836
|
135 |
assumes True
|
neuper@48836
|
136 |
yields r = r1 mod m1 \<and> r = r2 mod m2 *)
|
neuper@48862
|
137 |
fun chinese_remainder (r1, m1) (r2, m2) =
|
neuper@48862
|
138 |
(r1 mod m1) +
|
neuper@48862
|
139 |
(((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) * m1;
|
meindl_diana@48797
|
140 |
*}
|
meindl_diana@48797
|
141 |
|
neuper@48851
|
142 |
subsection {* creation of lists of primes for efficiency *}
|
meindl_diana@48797
|
143 |
ML{*
|
neuper@48830
|
144 |
(* is_prime :: int list \<Rightarrow> int \<Rightarrow> bool
|
neuper@48830
|
145 |
is_prime ps n = b
|
neuper@48830
|
146 |
assumes max ps < n \<and> n \<le> (max ps)^2 \<and>
|
neuper@48830
|
147 |
(* FIXME: really ^^^^^^^^^^^^^^^? *)
|
neuper@48831
|
148 |
(\<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
|
149 |
yields Primes.prime n *)
|
neuper@48834
|
150 |
fun is_prime [] _ = true
|
neuper@48834
|
151 |
| is_prime prs n =
|
neuper@48834
|
152 |
if n mod (List.hd prs) <> 0 then is_prime (List.tl prs) n else false
|
meindl_diana@48797
|
153 |
|
neuper@48830
|
154 |
(* make_primes is just local to primes_upto only:
|
neuper@48830
|
155 |
make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list
|
neuper@48830
|
156 |
make_primes ps last_p n = pps
|
neuper@48831
|
157 |
assumes last_p = maxs ps \<and> (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
158 |
(\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)
|
neuper@48836
|
159 |
yields n \<le> maxs pps \<and> (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
160 |
(\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
|
neuper@48834
|
161 |
fun make_primes ps last_p n =
|
neuper@48866
|
162 |
if n <= last ps then ps else
|
neuper@48813
|
163 |
if is_prime ps (last_p + 2)
|
neuper@48813
|
164 |
then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
|
neuper@48834
|
165 |
else make_primes ps (last_p + 2) n
|
neuper@48834
|
166 |
|
neuper@48831
|
167 |
(* primes_upto :: nat \<Rightarrow> nat list
|
neuper@48830
|
168 |
primes_upto n = ps
|
neuper@48831
|
169 |
assumes True
|
neuper@48836
|
170 |
yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
171 |
n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
|
neuper@48835
|
172 |
(\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
|
neuper@48836
|
173 |
fun primes_upto n = if n < 3 then [2] else make_primes [2, 3] 3 n
|
neuper@48813
|
174 |
|
neuper@48831
|
175 |
(* find the next prime greater p not dividing the number n:
|
neuper@48827
|
176 |
next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
|
neuper@48827
|
177 |
n1 next_prime_not_dvd n2 = p
|
neuper@48831
|
178 |
assumes True
|
neuper@48837
|
179 |
yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
|
neuper@48837
|
180 |
(\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
|
neuper@48835
|
181 |
infix next_prime_not_dvd;
|
neuper@48835
|
182 |
fun n1 next_prime_not_dvd n2 =
|
neuper@48813
|
183 |
let
|
neuper@48835
|
184 |
val ps = primes_upto (n1 + 1)
|
neuper@48813
|
185 |
val nxt = nth ps (List.length ps - 1);
|
meindl_diana@48797
|
186 |
in
|
neuper@48836
|
187 |
if n2 mod nxt <> 0 then nxt else nxt next_prime_not_dvd n2
|
meindl_diana@48797
|
188 |
end
|
meindl_diana@48797
|
189 |
*}
|
meindl_diana@48797
|
190 |
|
neuper@48851
|
191 |
subsection {* basic notions for univariate polynomials *}
|
neuper@48813
|
192 |
ML {*
|
neuper@48857
|
193 |
(* leading coefficient, includes drop_lc0_up ?FIXME130507 *)
|
neuper@48866
|
194 |
fun lcoeff_up (p: unipoly) = (last o drop_tl_zeros) p
|
meindl_diana@48797
|
195 |
|
neuper@48857
|
196 |
(* drop leading coefficients equal 0 *)
|
neuper@48857
|
197 |
fun drop_lc0_up (p: unipoly) = drop_tl_zeros p : unipoly;
|
meindl_diana@48797
|
198 |
|
neuper@48857
|
199 |
(* degree, includes drop_lc0_up ?FIXME130507 *)
|
neuper@48857
|
200 |
fun deg_up (p: unipoly) = ((Integer.add ~1) o length o drop_lc0_up) p;
|
meindl_diana@48797
|
201 |
|
neuper@48846
|
202 |
(* normalise a unipoly such that lcoeff_up mod m = 1.
|
neuper@48817
|
203 |
normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
|
neuper@48817
|
204 |
normalise [p_0, .., p_k] m = [q_0, .., q_k]
|
neuper@48827
|
205 |
assumes True
|
neuper@48817
|
206 |
yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i \<and> (p_k * t) mod m = 1 *)
|
neuper@48815
|
207 |
fun normalise (p: unipoly) (m: int) =
|
meindl_diana@48797
|
208 |
let
|
neuper@48846
|
209 |
val p = drop_lc0_up p
|
neuper@48846
|
210 |
val lcp = lcoeff_up p;
|
neuper@48815
|
211 |
fun norm nrm i =
|
neuper@48814
|
212 |
if i = 0
|
neuper@48815
|
213 |
then [mod_div (nth p i) lcp m] @ nrm
|
neuper@48815
|
214 |
else norm ([mod_div (nth p i) lcp m] @ nrm) (i - 1) ;
|
meindl_diana@48797
|
215 |
in
|
neuper@48815
|
216 |
if List.length p = 0 then p else norm [] (List.length p - 1)
|
neuper@48815
|
217 |
end;
|
neuper@48851
|
218 |
*}
|
meindl_diana@48797
|
219 |
|
neuper@48851
|
220 |
subsection {* addition, multiplication, division *}
|
neuper@48851
|
221 |
ML {*
|
neuper@48815
|
222 |
(* scalar multiplication *)
|
neuper@48815
|
223 |
infix %*
|
neuper@48839
|
224 |
fun (p: unipoly) %* n = map (Integer.mult n) p : unipoly
|
neuper@48815
|
225 |
|
neuper@48815
|
226 |
(* scalar divison *)
|
neuper@48815
|
227 |
infix %/
|
neuper@48817
|
228 |
fun (p: unipoly) %/ n =
|
neuper@48817
|
229 |
let fun swapf f a b = f b a (* swaps the arguments of a function *)
|
neuper@48839
|
230 |
in map (swapf (curry op div2) n) p : unipoly end;
|
meindl_diana@48797
|
231 |
|
neuper@48815
|
232 |
(* minus of polys *)
|
meindl_diana@48797
|
233 |
infix %-%
|
neuper@48839
|
234 |
fun (p1: unipoly) %-% (p2: unipoly) =
|
neuper@48839
|
235 |
map2 (curry op -) p1 (p2 @ replicate (List.length p1 - List.length p2) 0) : unipoly
|
meindl_diana@48797
|
236 |
|
neuper@48823
|
237 |
(* dvd for univariate polynomials *)
|
meindl_diana@48797
|
238 |
infix %|%
|
neuper@48845
|
239 |
(**)
|
neuper@48847
|
240 |
fun [d] %|% [p] = (p mod d = 0)
|
neuper@48839
|
241 |
| (ds: unipoly) %|% (ps: unipoly) =
|
neuper@48837
|
242 |
let
|
neuper@48846
|
243 |
val ds = drop_lc0_up ds;
|
neuper@48846
|
244 |
val ps = drop_lc0_up ps;
|
neuper@48839
|
245 |
val d000 = (replicate (List.length ps - List.length ds) 0) @ ds ;
|
neuper@48846
|
246 |
val quot = (lcoeff_up ps) div2 (lcoeff_up d000);
|
neuper@48846
|
247 |
val rest = drop_lc0_up (ps %-% (d000 %* quot));
|
neuper@48845
|
248 |
in
|
neuper@48846
|
249 |
if rest = [] then true else
|
neuper@48845
|
250 |
if rest <> [] andalso quot <> 0 andalso List.length ds <= List.length rest
|
neuper@48846
|
251 |
then ds %|% rest
|
neuper@48846
|
252 |
else false
|
neuper@48837
|
253 |
end
|
neuper@48851
|
254 |
*}
|
neuper@48851
|
255 |
|
neuper@48851
|
256 |
subsection {* normalisation and Landau-Mignotte bound *}
|
neuper@48851
|
257 |
ML {*
|
neuper@48855
|
258 |
(* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
|
neuper@48826
|
259 |
normalise [p_0, .., p_k] m = [q_0, .., q_k]
|
neuper@48826
|
260 |
assumes
|
neuper@48836
|
261 |
yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
|
neuper@48836
|
262 |
(where |^ x ^| means round x up to the next greater integer) *)
|
neuper@48855
|
263 |
fun centr_up (p: unipoly) (m: int) =
|
meindl_diana@48797
|
264 |
let
|
neuper@48823
|
265 |
val mi = m div2 2;
|
neuper@48823
|
266 |
val mid = if m mod mi = 0 then mi else mi + 1;
|
neuper@48823
|
267 |
fun centr m mid p_i = if mid < p_i then p_i - m else p_i;
|
neuper@48839
|
268 |
in map (centr m mid) p : unipoly end;
|
meindl_diana@48797
|
269 |
|
neuper@48836
|
270 |
(*** landau mignotte bound (Winkler p91)
|
neuper@48831
|
271 |
every coefficient of the gcd of a and b in Z[x] is bounded in absolute value by
|
neuper@48827
|
272 |
2^{min(m,n)} * gcd(a_m,b_n) * min(1/|a_m|
|
neuper@48836
|
273 |
* root(sum_{i=0}^{m}a_i^2) ,1/|b_n| * root( sum_{i=0}^{n}b_i^2 ) ***)
|
neuper@48826
|
274 |
|
neuper@48855
|
275 |
(* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
|
neuper@48826
|
276 |
sum_lmb [p_0, .., p_k] e = s
|
neuper@48826
|
277 |
assumes exp >= 0
|
neuper@48836
|
278 |
yields. p_0^e + p_1^e + ... + p_k^e *)
|
neuper@48823
|
279 |
fun sum_lmb (p: unipoly) exp = fold ((curry op +) o (Integer.pow exp)) p 0;
|
meindl_diana@48810
|
280 |
|
neuper@48855
|
281 |
(* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly \<Rightarrow> unipoly \<Rightarrow> int
|
neuper@48830
|
282 |
LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
|
neuper@48827
|
283 |
assumes True
|
neuper@48836
|
284 |
yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
|
neuper@48827
|
285 |
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
|
286 |
fun LANDAU_MIGNOTTE_bound (p1: unipoly) (p2: unipoly) =
|
neuper@48846
|
287 |
(Integer.pow (Integer.min (deg_up p1) (deg_up p2)) 2) * (abs (Integer.gcd (lcoeff_up p1) (lcoeff_up p2))) *
|
neuper@48846
|
288 |
(Int.min (abs((approx_root (sum_lmb p1 2)) div2 ~(lcoeff_up p1)),
|
neuper@48846
|
289 |
abs((approx_root (sum_lmb p2 2)) div2 ~(lcoeff_up p2))));
|
meindl_diana@48797
|
290 |
*}
|
meindl_diana@48797
|
291 |
|
neuper@48851
|
292 |
subsection {* modulo calculations on polynomials *}
|
meindl_diana@48797
|
293 |
ML{*
|
neuper@48855
|
294 |
(* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
|
neuper@48855
|
295 |
chinese_remainder_up (m1, m2) (p1, p2) = p
|
meindl_diana@48838
|
296 |
assume m1, m2 relatively prime
|
neuper@48855
|
297 |
yields p1 = p mod m1 \<and> p2 = p mod m2 *)
|
neuper@48855
|
298 |
fun chinese_remainder_up (m1, m2) (p1: unipoly, p2: unipoly) =
|
meindl_diana@48797
|
299 |
let
|
neuper@48862
|
300 |
fun chinese_rem (m1, m2) (p_1i, p_2i) = chinese_remainder (p_1i, m1) (p_2i, m2)
|
neuper@48823
|
301 |
in map (chinese_rem (m1, m2)) (p1 ~~ p2): unipoly end
|
meindl_diana@48797
|
302 |
|
neuper@48858
|
303 |
(* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
neuper@48858
|
304 |
mod_up [p1, p2, ..., pk] m = up
|
meindl_diana@48838
|
305 |
assume m > 0
|
meindl_diana@48838
|
306 |
yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
|
neuper@48858
|
307 |
infix mod_up
|
neuper@48858
|
308 |
fun (p : unipoly) mod_up m =
|
neuper@48855
|
309 |
let
|
neuper@48863
|
310 |
fun mod' m p = (curry op mod) p m
|
neuper@48863
|
311 |
in map (mod' m) p : unipoly end
|
meindl_diana@48797
|
312 |
|
meindl_diana@48819
|
313 |
(* euclidean algoritm in Z_p[x].
|
neuper@48858
|
314 |
mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
neuper@48858
|
315 |
mod_up_gcd p1 p2 m = g
|
neuper@48836
|
316 |
assumes True
|
neuper@48836
|
317 |
yields gcd ((p1 mod m), (p2 mod m)) = g *)
|
neuper@48858
|
318 |
fun mod_up_gcd (p1: unipoly) (p2: unipoly) (m: int) =
|
meindl_diana@48797
|
319 |
let
|
neuper@48858
|
320 |
val p1m = p1 mod_up m;
|
neuper@48858
|
321 |
val p2m = drop_lc0_up (p2 mod_up m);
|
neuper@48823
|
322 |
val p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
|
neuper@48846
|
323 |
val quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
|
neuper@48858
|
324 |
val rest = drop_lc0_up ((p1m %-% (p2n %* quot)) mod_up m)
|
neuper@48836
|
325 |
in
|
neuper@48836
|
326 |
if rest = [] then p2 else
|
neuper@48836
|
327 |
if length rest < List.length p2
|
neuper@48858
|
328 |
then mod_up_gcd p2 rest m
|
neuper@48858
|
329 |
else mod_up_gcd rest p2 m
|
neuper@48823
|
330 |
end;
|
neuper@48823
|
331 |
|
neuper@48855
|
332 |
(* primitive_up :: unipoly \<Rightarrow> unipoly
|
neuper@48855
|
333 |
primitive_up p = pp
|
neuper@48831
|
334 |
assumes True
|
neuper@48836
|
335 |
yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
|
neuper@48863
|
336 |
fun primitive_up ([c] : unipoly) =
|
neuper@48863
|
337 |
if c = 0 then ([0] : unipoly) else ([1] : unipoly)
|
neuper@48863
|
338 |
| primitive_up p =
|
neuper@48831
|
339 |
let
|
neuper@48831
|
340 |
val d = abs (Integer.gcds p)
|
neuper@48831
|
341 |
in
|
neuper@48831
|
342 |
if d = 1 then p else p %/ d
|
neuper@48831
|
343 |
end
|
meindl_diana@48797
|
344 |
*}
|
meindl_diana@48797
|
345 |
|
neuper@48855
|
346 |
subsection {* gcd_up, code from [1] p.93 *}
|
meindl_diana@48797
|
347 |
ML {*
|
meindl_diana@48838
|
348 |
(* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
meindl_diana@48838
|
349 |
try_new_prime_up p1 p2 d M P g p = new_g
|
neuper@48855
|
350 |
assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
|
neuper@48855
|
351 |
M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
|
neuper@48855
|
352 |
p1 is primitive \<and> p2 is primitive
|
neuper@48855
|
353 |
yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
|
neuper@48855
|
354 |
|
neuper@48855
|
355 |
argumentns "a b d M P g p" named according to [1] p.93 *)
|
neuper@48835
|
356 |
fun try_new_prime_up a b d M P g p =
|
neuper@48836
|
357 |
if P > M then g else
|
neuper@48833
|
358 |
let
|
neuper@48865
|
359 |
val p = p next_prime_not_dvd d
|
neuper@48865
|
360 |
val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
|
neuper@48865
|
361 |
%* (d mod p))
|
neuper@48865
|
362 |
mod_up p)
|
neuper@48865
|
363 |
p
|
neuper@48827
|
364 |
in
|
neuper@48846
|
365 |
if deg_up g_p < deg_up g
|
neuper@48827
|
366 |
then
|
neuper@48865
|
367 |
if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p
|
neuper@48827
|
368 |
else
|
neuper@48865
|
369 |
if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else
|
neuper@48827
|
370 |
let
|
neuper@48865
|
371 |
val P = P * p
|
neuper@48865
|
372 |
val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
|
neuper@48865
|
373 |
in try_new_prime_up a b d M P g p end
|
neuper@48827
|
374 |
end
|
neuper@48834
|
375 |
|
meindl_diana@48838
|
376 |
(* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
|
meindl_diana@48838
|
377 |
HENSEL_lifting_up p1 p2 d M p = gcd
|
neuper@48855
|
378 |
assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
|
neuper@48855
|
379 |
M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and>
|
neuper@48855
|
380 |
p1 is primitive \<and> p2 is primitive
|
neuper@48855
|
381 |
yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
|
neuper@48855
|
382 |
|
neuper@48855
|
383 |
argumentns "a b d M p" named according to [1] p.93 *)
|
neuper@48835
|
384 |
fun HENSEL_lifting_up a b d M p =
|
neuper@48823
|
385 |
let
|
neuper@48835
|
386 |
val p = p next_prime_not_dvd d
|
neuper@48865
|
387 |
val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
|
neuper@48823
|
388 |
in
|
neuper@48846
|
389 |
if deg_up g_p = 0 then [1] else
|
neuper@48835
|
390 |
let
|
neuper@48855
|
391 |
val g = primitive_up (try_new_prime_up a b d M p g_p p)
|
neuper@48827
|
392 |
in
|
neuper@48846
|
393 |
if (g %|% a) andalso (g %|% b) then g:unipoly else HENSEL_lifting_up a b d M p
|
neuper@48827
|
394 |
end
|
neuper@48827
|
395 |
end
|
neuper@48857
|
396 |
|
neuper@48855
|
397 |
(* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
|
neuper@48855
|
398 |
gcd_up a b = c
|
meindl_diana@48838
|
399 |
assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
|
meindl_diana@48838
|
400 |
a is primitive \<and> b is primitive
|
neuper@48836
|
401 |
yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
|
neuper@48855
|
402 |
fun gcd_up a b =
|
neuper@48855
|
403 |
let val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))
|
neuper@48855
|
404 |
in
|
neuper@48855
|
405 |
if a = b then a else
|
neuper@48855
|
406 |
HENSEL_lifting_up a b d (2 * d * LANDAU_MIGNOTTE_bound a b) 1
|
neuper@48855
|
407 |
end;
|
meindl_diana@48797
|
408 |
*}
|
meindl_diana@48797
|
409 |
|
neuper@48851
|
410 |
section {* gcd for multivariate polynomials *}
|
meindl_diana@48797
|
411 |
ML {*
|
meindl_diana@48797
|
412 |
type monom = (int * int list);
|
neuper@48846
|
413 |
type poly = monom list;
|
meindl_diana@48797
|
414 |
*}
|
meindl_diana@48797
|
415 |
|
neuper@48851
|
416 |
subsection {* auxiliary functions *}
|
meindl_diana@48797
|
417 |
ML {*
|
neuper@48846
|
418 |
fun land a b = a andalso b
|
neuper@48846
|
419 |
fun all_geq a b = fold land (map2 (curry op >=) a b) true
|
neuper@48850
|
420 |
fun all_geq' i es = fold land (map ((curry op <=) i) es) true
|
neuper@48850
|
421 |
fun maxs is = fold Integer.max is (hd is);
|
neuper@48850
|
422 |
fun div2_swapped d i = i div2 d
|
neuper@48850
|
423 |
fun nth_swapped i ls = nth ls i;
|
neuper@48849
|
424 |
fun drop_last ls = (rev o tl o rev) ls
|
neuper@48849
|
425 |
fun nth_ins n i xs = (*analogous to nth_drop: why not in library.ML?*)
|
neuper@48849
|
426 |
List.take (xs, n) @ [i] @ List.drop (xs, n);
|
neuper@48851
|
427 |
*}
|
neuper@48850
|
428 |
|
neuper@48851
|
429 |
subsection {* basic notions for multivariate polynomials *}
|
neuper@48851
|
430 |
ML {*
|
neuper@48849
|
431 |
fun lcoeff (p: poly) = (fst o hd o rev) p
|
neuper@48850
|
432 |
fun lexp (p: poly) = (snd o hd o rev) p
|
neuper@48850
|
433 |
fun lmonom (p: poly) = (hd o rev) p
|
neuper@48850
|
434 |
|
neuper@48849
|
435 |
fun cero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;
|
neuper@48850
|
436 |
|
neuper@48851
|
437 |
(* drop leading coefficients equal 0 TODO: signature?*)
|
neuper@48851
|
438 |
fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then cero_poly (length es) else p
|
neuper@48851
|
439 |
| drop_lc0 p = if lcoeff p = 0 then drop_lc0 (drop_last p) else p
|
neuper@48851
|
440 |
|
neuper@48851
|
441 |
(* gcd of coefficients WN find right location in file *)
|
neuper@48851
|
442 |
fun gcd_coeff (up: unipoly) = abs (Integer.gcds up)
|
neuper@48851
|
443 |
|
neuper@48851
|
444 |
(* given p = 3x+4z = [(3,[1,0]),(4,[0,1])] with x<z
|
neuper@48851
|
445 |
add_variable p 2 => [(3,[1,0,0]),(4,[0,0,1])] with x<y<z *)
|
neuper@48851
|
446 |
fun add_variable (p: poly) var = map (apsnd (nth_ins var 0)) p : poly
|
neuper@48851
|
447 |
*}
|
neuper@48851
|
448 |
|
neuper@48851
|
449 |
subsection {* ordering and other auxiliary funcions *}
|
neuper@48851
|
450 |
ML {*
|
neuper@48857
|
451 |
(* monomial order: lexicographic order on exponents *)
|
neuper@48850
|
452 |
fun lex_ord ((_, a): monom) ((_, b): monom) =
|
neuper@48857
|
453 |
let val d = drop_tl_zeros (a %-% b)
|
neuper@48850
|
454 |
in
|
neuper@48866
|
455 |
if d = [] then false else last d > 0
|
neuper@48850
|
456 |
end
|
neuper@48849
|
457 |
fun lex_ord' ((_, a): monom, (_, b): monom) =
|
neuper@48857
|
458 |
let val d = drop_tl_zeros (a %-% b)
|
neuper@48849
|
459 |
in
|
neuper@48857
|
460 |
if d = [] then EQUAL
|
neuper@48866
|
461 |
else if last d > 0 then GREATER else LESS
|
neuper@48849
|
462 |
end
|
neuper@48850
|
463 |
|
neuper@48848
|
464 |
(* add monomials in ordered polynomal *)
|
neuper@48850
|
465 |
fun add_monoms (p: poly) =
|
neuper@48846
|
466 |
let
|
neuper@48847
|
467 |
fun add [(0, es)] [] = cero_poly (length es)
|
neuper@48847
|
468 |
| add [(0, _)] (new: poly) = new
|
neuper@48847
|
469 |
| add [m] (new: poly) = new @ [m]
|
neuper@48847
|
470 |
| add ((c1, es1) :: (c2, es2) :: ms) (new: poly) =
|
neuper@48846
|
471 |
if es1 = es2
|
neuper@48847
|
472 |
then add ([(c1 + c2, es1)] @ ms) new
|
neuper@48846
|
473 |
else
|
neuper@48846
|
474 |
if c1 = 0
|
neuper@48847
|
475 |
then add ((c2, es2) :: ms) new
|
neuper@48847
|
476 |
else add ((c2, es2) :: ms) (new @ [(c1, es1)])
|
neuper@48850
|
477 |
in add p [] : poly end
|
neuper@48850
|
478 |
|
neuper@48848
|
479 |
(* brings the polynom in correct order and adds monomials *)
|
neuper@48851
|
480 |
fun order (p: poly) = (add_monoms o (sort lex_ord')) p;
|
neuper@48850
|
481 |
|
neuper@48850
|
482 |
(* largest exponent of variable at location var *)
|
neuper@48851
|
483 |
fun max_deg_var (p: poly) var = (maxs o (map ((nth_swapped var) o snd))) p;
|
neuper@48851
|
484 |
*}
|
neuper@48850
|
485 |
|
neuper@48851
|
486 |
subsection {* evaluation, translation uni--multivariate, etc *}
|
neuper@48851
|
487 |
ML {*
|
neuper@48851
|
488 |
(* evaluate a polynomial in a certain variable and remove this from the exp-list *)
|
neuper@48855
|
489 |
fun eval (p: poly) var value =
|
neuper@48851
|
490 |
let
|
neuper@48851
|
491 |
fun eval [] _ _ new = order new
|
neuper@48851
|
492 |
| eval ((c, es)::ps) var value new =
|
neuper@48851
|
493 |
eval ps var value (new @ [((Integer.pow (nth es var) value) * c, nth_drop var es)]);
|
neuper@48851
|
494 |
in eval p var value [] end
|
neuper@48851
|
495 |
|
neuper@48851
|
496 |
(* transform a multivariate to a univariate polynomial *)
|
neuper@48851
|
497 |
fun multi_to_uni (p as (_, es)::_ : poly) =
|
neuper@48851
|
498 |
if length es = 1
|
neuper@48851
|
499 |
then
|
neuper@48851
|
500 |
let fun transfer ([]: poly) (up: unipoly) = up
|
neuper@48851
|
501 |
| transfer (mp as (c, es)::ps) up =
|
neuper@48851
|
502 |
if length up = hd es
|
neuper@48851
|
503 |
then transfer ps (up @ [c])
|
neuper@48851
|
504 |
else transfer mp (up @ [0])
|
neuper@48851
|
505 |
in transfer (order p) [] end
|
neuper@48851
|
506 |
else error "Polynom has more than one variable!";
|
neuper@48851
|
507 |
|
neuper@48851
|
508 |
(* transform a univariate to a multivariate polynomial *)
|
neuper@48851
|
509 |
fun uni_to_multi (up: unipoly) =
|
neuper@48851
|
510 |
let fun transfer ([]: unipoly) (mp: poly) (_: int) = add_monoms mp
|
neuper@48851
|
511 |
| transfer up mp var =
|
neuper@48851
|
512 |
transfer (tl up) (mp @ [(hd up, [var])]) (var + 1)
|
neuper@48851
|
513 |
in transfer up [] 0 end
|
neuper@48853
|
514 |
|
neuper@48851
|
515 |
(* multiply a poly with a variable represented by a unipoly (at position var)
|
neuper@48851
|
516 |
e.g new z: (3x + 2y)*(z - 4) TODO fold / map ? *)
|
neuper@48851
|
517 |
fun mult_with_new_var ([]: poly) (_: unipoly) (_: int) = []
|
neuper@48851
|
518 |
| mult_with_new_var mp up var =
|
neuper@48851
|
519 |
let
|
neuper@48851
|
520 |
fun mult ([]: poly) (_: unipoly) (_: int) (new: poly) (_: int) = add_monoms new
|
neuper@48851
|
521 |
| mult mp up var new _ =
|
neuper@48851
|
522 |
let
|
neuper@48851
|
523 |
(* the inner loop *)
|
neuper@48851
|
524 |
fun mult' (_: poly) ([]: unipoly) (_) (new': poly) (_: int) = order new'
|
neuper@48851
|
525 |
| mult' (mp' as (c, es)::ps) up' var new' c2' =
|
neuper@48851
|
526 |
let val (first, last) = chop var es
|
neuper@48851
|
527 |
in mult' mp' (tl up') var
|
neuper@48851
|
528 |
(new' @ [(c * hd up', first @ [c2'] @ last)]) (c2' + 1)
|
neuper@48851
|
529 |
end
|
neuper@48851
|
530 |
in mult (tl mp) up var (new @ (mult' mp up var [] 0)) (0) end
|
neuper@48851
|
531 |
in mult mp up var [] 0 : poly end
|
neuper@48853
|
532 |
|
neuper@48851
|
533 |
fun greater_deg (es1: int list) (es2: int list) =
|
neuper@48851
|
534 |
if es1 = [] andalso es2 =[] then 0
|
neuper@48851
|
535 |
else if (nth es1 (length es1 -1)) = (nth es2 (length es1 -1) )
|
neuper@48851
|
536 |
then greater_deg (nth_drop (length es1 -1) es1) (nth_drop (length es1 -1) es2)
|
neuper@48851
|
537 |
else if (nth es1 (length es1 -1)) > (nth es2 (length es1 -1))
|
neuper@48851
|
538 |
then 1 else 2
|
neuper@48853
|
539 |
|
neuper@48851
|
540 |
(* Winkler p.95 *)
|
neuper@48851
|
541 |
fun find_new_r (p1 as (_, es1)::_ : poly) (p2 as (_, es2)::_ : poly) old =
|
neuper@48855
|
542 |
let val p1' as (_, es1')::_ = eval p1 (length es1 - 2) (old + 1)
|
neuper@48855
|
543 |
val p2' as (_, es2')::_ = eval p2 (length es2 - 2) (old + 1);
|
neuper@48851
|
544 |
in
|
neuper@48851
|
545 |
if max_deg_var p1' (length es1' - 1) = max_deg_var p1 (length es1 - 1) orelse
|
neuper@48851
|
546 |
max_deg_var p2' (length es2' - 1) = max_deg_var p2 (length es1 - 1)
|
neuper@48851
|
547 |
then old + 1
|
neuper@48851
|
548 |
else find_new_r p1 p2 (old + 1)
|
neuper@48851
|
549 |
end
|
neuper@48851
|
550 |
*}
|
neuper@48851
|
551 |
|
neuper@48851
|
552 |
subsection {* addition, multiplication, division *}
|
neuper@48851
|
553 |
ML {*
|
neuper@48850
|
554 |
(* resulting 0 coefficients are removed by add_monoms *)
|
meindl_diana@48797
|
555 |
infix %%/
|
neuper@48850
|
556 |
fun (p: poly) %%/ d = (add_monoms o (map (apfst (div2_swapped d)))) p;
|
meindl_diana@48797
|
557 |
|
meindl_diana@48841
|
558 |
infix %%*
|
neuper@48850
|
559 |
fun (p: poly) %%* i = (add_monoms o (map (apfst (Integer.mult i)))) p;
|
meindl_diana@48841
|
560 |
|
meindl_diana@48797
|
561 |
infix %%+%%
|
neuper@48853
|
562 |
(* the same algorithms as done by hand
|
neuper@48846
|
563 |
fun ([]: poly) %%+%% (mvp2: poly) = mvp2
|
neuper@48850
|
564 |
| p1 %%+%% [] = p1
|
neuper@48850
|
565 |
| p1 %%+%% p2 =
|
neuper@48850
|
566 |
let
|
neuper@48851
|
567 |
fun plus [] p2 new = order (new @ p2)
|
neuper@48851
|
568 |
| plus p1 [] new = order (new @ p1)
|
neuper@48850
|
569 |
| plus (p1 as (c1, es1)::ps1) (p2 as (c2, es2)::ps2) new =
|
neuper@48850
|
570 |
if greater_deg es1 es2 = 0
|
neuper@48850
|
571 |
then plus ps1 ps2 (new @ [((c1 + c2), es1)])
|
neuper@48850
|
572 |
else
|
neuper@48850
|
573 |
if greater_deg es1 es2 = 1
|
neuper@48850
|
574 |
then plus p1 ps2 (new @ [(c2, es2)])
|
neuper@48850
|
575 |
else plus ps1 p2 (new @ [(c1, es1)])
|
neuper@48850
|
576 |
in plus p1 p2 [] end
|
neuper@48853
|
577 |
*)
|
neuper@48853
|
578 |
fun (p1 : poly) %%+%% (p2 : poly) = order (p1 @ p2)
|
neuper@48853
|
579 |
|
meindl_diana@48797
|
580 |
infix %%-%%
|
neuper@48850
|
581 |
fun p1 %%-%% (p2: poly) = p1 %%+%% (map (apfst (Integer.mult ~1)) p2)
|
neuper@48851
|
582 |
|
meindl_diana@48797
|
583 |
infix %%*%%
|
neuper@48853
|
584 |
(* the same algorithms as done by hand
|
neuper@48850
|
585 |
fun (p1: poly) %%*%% (p2: poly) =
|
neuper@48850
|
586 |
let
|
neuper@48851
|
587 |
fun mult ([]: poly) (_: poly) (_: poly) (new: poly) = order new
|
neuper@48850
|
588 |
| mult p1 [] regular_p2 new = mult (tl p1) regular_p2 regular_p2 new
|
neuper@48850
|
589 |
| mult (p1 as (c1, es1)::_) ((c2, es2)::ps2) regular_p2 (new: poly) =
|
neuper@48850
|
590 |
mult p1 ps2 regular_p2 (new @ [(c1 * c2, map2 Integer.add es1 es2)])
|
neuper@48850
|
591 |
in
|
neuper@48850
|
592 |
if length p1 > length p2
|
neuper@48850
|
593 |
then mult p1 p2 p2 []
|
neuper@48850
|
594 |
else mult p2 p1 p1 []
|
neuper@48850
|
595 |
end
|
neuper@48853
|
596 |
*)
|
neuper@48853
|
597 |
fun mult_monom ((c1, es1): monom) ((c2, es2): monom) =
|
neuper@48853
|
598 |
(c1 * c2, map2 Integer.add es1 es2): monom;
|
neuper@48853
|
599 |
fun mult_monoms (ms: poly) m = map (mult_monom m) ms : poly
|
neuper@48853
|
600 |
fun p1 %%*%% (p2 : poly) = order (fold (curry op @) (map (mult_monoms p1) p2) [])
|
neuper@48853
|
601 |
|
meindl_diana@48797
|
602 |
infix %%|%%
|
neuper@48851
|
603 |
(* the same algorithms as done by hand; TODO fold / map ? *)
|
neuper@48850
|
604 |
fun ([(c1, es1)]: poly) %%|%% ([(c2, es2)]: poly) = all_geq es2 es1 andalso c2 mod c1 = 0
|
neuper@48850
|
605 |
| _ %%|%% [(0,_)] = true
|
neuper@48850
|
606 |
| p1 %%|%% p2 =
|
neuper@48850
|
607 |
if [lmonom p1] %%|%% [lmonom p2]
|
neuper@48850
|
608 |
then
|
neuper@48850
|
609 |
p1 %%|%%
|
neuper@48850
|
610 |
(p2 %%-%%
|
neuper@48850
|
611 |
([((lcoeff p2) div2 (lcoeff p1), (lexp p2) %-% (lexp p1))]
|
neuper@48850
|
612 |
%%*%% p1))
|
neuper@48850
|
613 |
else false
|
neuper@48853
|
614 |
|
neuper@48851
|
615 |
(* %%/%% :: poly \<Rightarrow> poly \<Rightarrow> (poly \<times> poly)
|
neuper@48851
|
616 |
a %%/%% b = (c, r)
|
neuper@48851
|
617 |
assumes b \<noteq> [] \<and>
|
neuper@48851
|
618 |
yields THE c r. c *\<^isub>p b +\<^isub>p r = a *)
|
neuper@48851
|
619 |
infix %%/%%
|
neuper@48851
|
620 |
fun (p1: poly) %%/%% (p2: poly) =
|
neuper@48851
|
621 |
let
|
neuper@48851
|
622 |
fun div_up p1 p2 quot_old =
|
neuper@48851
|
623 |
let
|
neuper@48851
|
624 |
val p1 as (c, es)::_ = drop_lc0 (order p1);
|
neuper@48851
|
625 |
val p2 = drop_lc0 (order p2);
|
neuper@48851
|
626 |
val [c] = cero_poly (length es);
|
neuper@48851
|
627 |
val d = (replicate (List.length p1 - List.length p2) c) @ p2 ;
|
neuper@48851
|
628 |
val quot = [(lcoeff p1 div2 lcoeff d, lexp p1 %-% lexp p2)];
|
neuper@48851
|
629 |
val rest = drop_lc0 (p1 %%-%% (d %%*%% quot));
|
neuper@48851
|
630 |
in
|
neuper@48851
|
631 |
if rest <> [c] andalso all_geq' 0 (lexp rest %-% lexp p2)
|
neuper@48851
|
632 |
then div_up rest p2 (quot @ quot_old)
|
neuper@48851
|
633 |
else (quot @ quot_old : poly, rest)
|
neuper@48851
|
634 |
end
|
neuper@48851
|
635 |
in div_up p1 p2 [] end
|
meindl_diana@48797
|
636 |
*}
|
neuper@48848
|
637 |
|
neuper@48846
|
638 |
subsection {* Newton interpolation *}
|
meindl_diana@48797
|
639 |
ML{* (* first step *)
|
neuper@48851
|
640 |
fun newton_first (x: int list) (f: poly list) (ord: int) =
|
neuper@48851
|
641 |
let
|
neuper@48851
|
642 |
val new_vals = [(hd x) * ~1, 1];
|
neuper@48851
|
643 |
val p = (add_variable (hd f) ord) %%+%%
|
neuper@48851
|
644 |
((mult_with_new_var
|
neuper@48851
|
645 |
(((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x))) new_vals ord);
|
neuper@48851
|
646 |
val steps = [((nth f 1) %%-%% (hd f)) %%/ ((nth x 1) - (hd x))]
|
neuper@48851
|
647 |
in (p, new_vals, steps) end
|
meindl_diana@48797
|
648 |
|
neuper@48851
|
649 |
fun next_step ([]: poly list) (new_steps: poly list) (_: int list) = new_steps
|
neuper@48851
|
650 |
| next_step steps new_steps x =
|
neuper@48851
|
651 |
next_step (tl steps)
|
neuper@48851
|
652 |
(new_steps @
|
neuper@48866
|
653 |
[((last new_steps) %%-%% (hd steps)) %%/ ((last x) - (nth x (length x - 3)))])
|
neuper@48851
|
654 |
(nth_drop (length x -2) x)
|
neuper@48851
|
655 |
|
neuper@48851
|
656 |
fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) ord =
|
meindl_diana@48797
|
657 |
if length x = 2
|
neuper@48851
|
658 |
then
|
neuper@48851
|
659 |
let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] ord
|
neuper@48851
|
660 |
in (p', new_vals, steps) end
|
neuper@48851
|
661 |
else
|
neuper@48850
|
662 |
let
|
neuper@48851
|
663 |
val new_vals =
|
neuper@48851
|
664 |
multi_to_uni ((uni_to_multi up) %%*%% (uni_to_multi [(nth x (length x -2) ) * ~1, 1]));
|
neuper@48851
|
665 |
val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/
|
neuper@48866
|
666 |
((last x) - (nth x (length x - 2)))];
|
neuper@48851
|
667 |
val steps = next_step steps new_steps x;
|
neuper@48866
|
668 |
val p' = p %%+%% (mult_with_new_var (last steps) new_vals ord);
|
neuper@48851
|
669 |
in (order p', new_vals, steps) end
|
meindl_diana@48841
|
670 |
*}
|
meindl_diana@48797
|
671 |
|
neuper@48851
|
672 |
subsection {* gcd_poly algorithm, code for [1] p.95 *}
|
meindl_diana@48797
|
673 |
ML {*
|
neuper@48851
|
674 |
(* naming of n, M, m, r, ... follows Winkler p. 95 *)
|
neuper@48851
|
675 |
fun gcd_poly' (a as (c1, es1)::ps1 : poly) (b as (c2, es2)::ps2 : poly) n r =
|
neuper@48849
|
676 |
if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else
|
neuper@48851
|
677 |
if n - 1 = 0
|
neuper@48851
|
678 |
then
|
neuper@48851
|
679 |
let val (a', b') = (multi_to_uni a, multi_to_uni b)
|
neuper@48851
|
680 |
in
|
neuper@48851
|
681 |
(* a b are made primitive and this is compensated by multiplying with gcd of divisors *)
|
neuper@48855
|
682 |
uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %*
|
neuper@48851
|
683 |
(abs (Integer.gcd (gcd_coeff a') (gcd_coeff b'))))
|
neuper@48851
|
684 |
end
|
neuper@48848
|
685 |
else
|
neuper@48851
|
686 |
let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
|
neuper@48848
|
687 |
in try_new_r a b n M r [] [] end
|
meindl_diana@48816
|
688 |
|
neuper@48851
|
689 |
and try_new_r a b n M r r_list steps =
|
meindl_diana@48816
|
690 |
let
|
meindl_diana@48810
|
691 |
val r = find_new_r a b r;
|
meindl_diana@48810
|
692 |
val r_list = r_list @ [r];
|
neuper@48855
|
693 |
val g_r = gcd_poly' (order (eval a (n - 2) r))
|
neuper@48855
|
694 |
(order (eval b (n - 2) r)) (n - 1) 0
|
meindl_diana@48810
|
695 |
val steps = steps @ [g_r];
|
neuper@48851
|
696 |
in HENSEL_lifting a b n M 1 r r_list steps g_r ( cero_poly n ) [1] end
|
meindl_diana@48816
|
697 |
|
neuper@48851
|
698 |
and HENSEL_lifting a b n M m r r_list steps g g_n mult =
|
meindl_diana@48816
|
699 |
if m > M then if g_n %%|%% a andalso g_n %%|%% b then g_n else
|
meindl_diana@48816
|
700 |
try_new_r a b n M r r_list steps
|
meindl_diana@48810
|
701 |
else
|
meindl_diana@48810
|
702 |
let
|
meindl_diana@48810
|
703 |
val r = find_new_r a b r;
|
meindl_diana@48810
|
704 |
val r_list = r_list @ [r];
|
neuper@48855
|
705 |
val g_r = gcd_poly' (order (eval a (n - 2) r))
|
neuper@48855
|
706 |
(order (eval b (n - 2) r)) (n - 1) 0
|
meindl_diana@48810
|
707 |
in
|
neuper@48849
|
708 |
if lex_ord (lmonom g) (lmonom g_r)
|
neuper@48855
|
709 |
then HENSEL_lifting a b n M 1 r r_list steps g g_n mult
|
neuper@48850
|
710 |
else if (lexp g) = (lexp g_r)
|
meindl_diana@48810
|
711 |
then let
|
neuper@48851
|
712 |
val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
|
meindl_diana@48810
|
713 |
in
|
neuper@48851
|
714 |
if (nth steps (length steps - 1)) = (cero_poly (n - 1))
|
neuper@48851
|
715 |
then HENSEL_lifting a b n M (M + 1) r r_list steps g_r g_n new
|
neuper@48851
|
716 |
else HENSEL_lifting a b n M (m + 1) r r_list steps g_r g_n new
|
meindl_diana@48810
|
717 |
end
|
neuper@48851
|
718 |
else HENSEL_lifting a b n M (m + 1) r r_list steps g g_n mult
|
meindl_diana@48810
|
719 |
end
|
neuper@48853
|
720 |
|
neuper@48851
|
721 |
(* gcd_poly :: poly \<Rightarrow> poly \<Rightarrow> poly
|
neuper@48851
|
722 |
gcd_poly a b = ((a', b'), c)
|
neuper@48851
|
723 |
assumes a \<noteq> [] \<and> b \<noteq> [] \<and>
|
neuper@48851
|
724 |
yields a = a' *\<^isub>p c \<and> b = b' *\<^isub>p c \<and>
|
neuper@48851
|
725 |
(\<forall>c'. (c' dvd\<^isub>p a \<and> c' dvd\<^isub>p b) \<longrightarrow> c' \<le>\<^isub>p c) *)
|
neuper@48851
|
726 |
fun gcd_poly (a as (_, es)::_ : poly) b =
|
neuper@48851
|
727 |
let val c = gcd_poly' a b (length es) 0
|
neuper@48851
|
728 |
val (a': poly, _) = a %%/%% c
|
neuper@48851
|
729 |
val (b': poly, _) = b %%/%% c
|
neuper@48851
|
730 |
in ((a', b'), c) end
|
meindl_diana@48797
|
731 |
*}
|
meindl_diana@48797
|
732 |
|
neuper@48851
|
733 |
section {* example from the last mail *}
|
neuper@48851
|
734 |
ML {* (* test, see text below *)
|
neuper@48851
|
735 |
val a = [(1,[2, 0]), (~1,[0, 2])];
|
neuper@48851
|
736 |
val b = [(1,[2, 0]), (~1,[1, 1])] ;
|
neuper@48851
|
737 |
val ((a', b'), c) = gcd_poly a b;
|
neuper@48851
|
738 |
|
neuper@48851
|
739 |
a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)
|
neuper@48851
|
740 |
c = [(~1, [1, 0]), (1, [0, 1])];
|
neuper@48851
|
741 |
a' = [(~1, [1, 0]), (~1, [0, 1])];
|
neuper@48851
|
742 |
b' = [(~1, [1, 0])];
|
neuper@48851
|
743 |
*}
|
neuper@48851
|
744 |
|
neuper@48851
|
745 |
text {* last mail to Tobias Nipkow:
|
neuper@48851
|
746 |
Eine erste Idee, wie die Integration der Diplomarbeit für einen Benutzer
|
neuper@48851
|
747 |
von Isabelle aussehen könnte, wäre zum Beispiel im
|
neuper@48851
|
748 |
|
neuper@48851
|
749 |
lemma cancel:
|
neuper@48851
|
750 |
assumes asm3: "x^2 - x*y \<noteq> 0" and asm4: "x \<noteq> 0"
|
neuper@48851
|
751 |
shows "(x^2 - y^2) / (x^2 - x*y) = (x + y) / (x::real)"
|
neuper@48851
|
752 |
apply (insert asm3 asm4)
|
neuper@48851
|
753 |
apply simp
|
neuper@48851
|
754 |
sorry
|
neuper@48851
|
755 |
|
neuper@48851
|
756 |
die Assumptions
|
neuper@48851
|
757 |
|
neuper@48851
|
758 |
asm1: "(x^2 - y^2) = (x + y) * (x - y)" and asm2: "x^2 - x*y = x * (x - y)"
|
neuper@48851
|
759 |
|
neuper@48851
|
760 |
im Hintergrund automatisch zu erzeugen (mit der Garantie, dass "(x - y)" der GCD ist)
|
neuper@48851
|
761 |
und sie dem Simplifier (für die Rule nonzero_mult_divide_mult_cancel_right) zur
|
neuper@48851
|
762 |
Verfügung zu stellen, sodass anstelle von "sorry" einfach "done" stehen kann.
|
neuper@48851
|
763 |
Und weiters wäre eventuell asm3 zu "x - y \<noteq> 0" zu vereinfachen, eine
|
neuper@48851
|
764 |
Rewriteorder zum Herstellen einer Normalform festzulegen, etc.
|
neuper@48851
|
765 |
*}
|
meindl_diana@48797
|
766 |
end
|
meindl_diana@48797
|
767 |
|
meindl_diana@48797
|
768 |
|