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 |
|
neuper@52073
|
9 |
theory GCD_Poly_ML 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;
|
neuper@48870
|
100 |
|
neuper@48870
|
101 |
(* swaps the arguments of a function *)
|
neuper@48870
|
102 |
fun swapf f a b = f b a
|
meindl_diana@48797
|
103 |
*}
|
meindl_diana@48797
|
104 |
|
meindl_diana@48797
|
105 |
subsection {* modulo calculations for integers *}
|
meindl_diana@48797
|
106 |
ML {*
|
neuper@48836
|
107 |
(* mod_inv :: int \<Rightarrow> nat
|
neuper@48817
|
108 |
mod_inv r m = s
|
neuper@48836
|
109 |
assumes True
|
neuper@48836
|
110 |
yields r * s mod m = 1 *)
|
meindl_diana@48797
|
111 |
fun mod_inv r m =
|
meindl_diana@48797
|
112 |
let
|
meindl_diana@48797
|
113 |
fun modi (r, rold, m, rinv) =
|
neuper@48836
|
114 |
if m <= rinv then 0 else
|
meindl_diana@48797
|
115 |
if r mod m = 1
|
meindl_diana@48797
|
116 |
then rinv
|
meindl_diana@48797
|
117 |
else modi (rold * (rinv + 1), rold, m, rinv + 1)
|
neuper@48861
|
118 |
in modi (r, r, m, 1) end
|
neuper@48862
|
119 |
|
neuper@48865
|
120 |
(* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
|
neuper@48817
|
121 |
mod_div a b m = c
|
neuper@48836
|
122 |
assumes True
|
neuper@48865
|
123 |
yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
|
neuper@48862
|
124 |
fun mod_div a b m = a * (mod_inv b m) mod m
|
meindl_diana@48797
|
125 |
|
neuper@48836
|
126 |
(* required just for one approximation:
|
neuper@48836
|
127 |
approx_root :: nat \<Rightarrow> nat
|
neuper@48817
|
128 |
approx_root a = r
|
neuper@48836
|
129 |
assumes True
|
neuper@48836
|
130 |
yields r * r \<ge> a *)
|
neuper@48813
|
131 |
fun approx_root a =
|
neuper@48836
|
132 |
let
|
neuper@48836
|
133 |
fun root1 a n = if n * n < a then root1 a (n + 1) else n
|
neuper@48836
|
134 |
in root1 a 1 end
|
meindl_diana@48797
|
135 |
|
neuper@48862
|
136 |
(* chinese_remainder :: int * nat \<Rightarrow> int * nat \<Rightarrow> nat
|
neuper@48862
|
137 |
chinese_remainder (r1, m1) (r2, m2) = r
|
neuper@48836
|
138 |
assumes True
|
neuper@48836
|
139 |
yields r = r1 mod m1 \<and> r = r2 mod m2 *)
|
neuper@48862
|
140 |
fun chinese_remainder (r1, m1) (r2, m2) =
|
neuper@48862
|
141 |
(r1 mod m1) +
|
neuper@48875
|
142 |
(((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) *
|
neuper@48875
|
143 |
m1;
|
meindl_diana@48797
|
144 |
*}
|
meindl_diana@48797
|
145 |
|
neuper@48851
|
146 |
subsection {* creation of lists of primes for efficiency *}
|
meindl_diana@48797
|
147 |
ML{*
|
neuper@48830
|
148 |
(* is_prime :: int list \<Rightarrow> int \<Rightarrow> bool
|
neuper@48830
|
149 |
is_prime ps n = b
|
neuper@48830
|
150 |
assumes max ps < n \<and> n \<le> (max ps)^2 \<and>
|
neuper@48830
|
151 |
(* FIXME: really ^^^^^^^^^^^^^^^? *)
|
neuper@48831
|
152 |
(\<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
|
153 |
yields Primes.prime n *)
|
neuper@48834
|
154 |
fun is_prime [] _ = true
|
neuper@48834
|
155 |
| is_prime prs n =
|
neuper@48834
|
156 |
if n mod (List.hd prs) <> 0 then is_prime (List.tl prs) n else false
|
meindl_diana@48797
|
157 |
|
neuper@48830
|
158 |
(* make_primes is just local to primes_upto only:
|
neuper@48830
|
159 |
make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list
|
neuper@48830
|
160 |
make_primes ps last_p n = pps
|
neuper@48831
|
161 |
assumes last_p = maxs ps \<and> (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
162 |
(\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)
|
neuper@48836
|
163 |
yields n \<le> maxs pps \<and> (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
164 |
(\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
|
neuper@48834
|
165 |
fun make_primes ps last_p n =
|
neuper@48866
|
166 |
if n <= last ps then ps else
|
neuper@48813
|
167 |
if is_prime ps (last_p + 2)
|
neuper@48813
|
168 |
then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
|
neuper@48834
|
169 |
else make_primes ps (last_p + 2) n
|
neuper@48834
|
170 |
|
neuper@48831
|
171 |
(* primes_upto :: nat \<Rightarrow> nat list
|
neuper@48830
|
172 |
primes_upto n = ps
|
neuper@48831
|
173 |
assumes True
|
neuper@48836
|
174 |
yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
|
neuper@48831
|
175 |
n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
|
neuper@48835
|
176 |
(\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
|
neuper@48836
|
177 |
fun primes_upto n = if n < 3 then [2] else make_primes [2, 3] 3 n
|
neuper@48813
|
178 |
|
neuper@48831
|
179 |
(* find the next prime greater p not dividing the number n:
|
neuper@48827
|
180 |
next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
|
neuper@48827
|
181 |
n1 next_prime_not_dvd n2 = p
|
neuper@48831
|
182 |
assumes True
|
neuper@48837
|
183 |
yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
|
neuper@48837
|
184 |
(\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
|
neuper@48835
|
185 |
infix next_prime_not_dvd;
|
neuper@48835
|
186 |
fun n1 next_prime_not_dvd n2 =
|
neuper@48813
|
187 |
let
|
neuper@48835
|
188 |
val ps = primes_upto (n1 + 1)
|
neuper@48813
|
189 |
val nxt = nth ps (List.length ps - 1);
|
meindl_diana@48797
|
190 |
in
|
neuper@48836
|
191 |
if n2 mod nxt <> 0 then nxt else nxt next_prime_not_dvd n2
|
meindl_diana@48797
|
192 |
end
|
meindl_diana@48797
|
193 |
*}
|
meindl_diana@48797
|
194 |
|
neuper@48851
|
195 |
subsection {* basic notions for univariate polynomials *}
|
neuper@48813
|
196 |
ML {*
|
neuper@48857
|
197 |
(* leading coefficient, includes drop_lc0_up ?FIXME130507 *)
|
neuper@48866
|
198 |
fun lcoeff_up (p: unipoly) = (last o drop_tl_zeros) p
|
meindl_diana@48797
|
199 |
|
neuper@48857
|
200 |
(* drop leading coefficients equal 0 *)
|
neuper@52063
|
201 |
fun drop_lc0_up (p: unipoly) =drop_tl_zeros p : unipoly;
|
meindl_diana@48797
|
202 |
|
neuper@48857
|
203 |
(* degree, includes drop_lc0_up ?FIXME130507 *)
|
neuper@48857
|
204 |
fun deg_up (p: unipoly) = ((Integer.add ~1) o length o drop_lc0_up) p;
|
meindl_diana@48797
|
205 |
|
neuper@48846
|
206 |
(* normalise a unipoly such that lcoeff_up mod m = 1.
|
neuper@48817
|
207 |
normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
|
neuper@48817
|
208 |
normalise [p_0, .., p_k] m = [q_0, .., q_k]
|
neuper@48827
|
209 |
assumes True
|
neuper@48817
|
210 |
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
|
211 |
fun normalise (p: unipoly) (m: int) =
|
meindl_diana@48797
|
212 |
let
|
neuper@48846
|
213 |
val p = drop_lc0_up p
|
neuper@48846
|
214 |
val lcp = lcoeff_up p;
|
neuper@48815
|
215 |
fun norm nrm i =
|
neuper@48814
|
216 |
if i = 0
|
neuper@48815
|
217 |
then [mod_div (nth p i) lcp m] @ nrm
|
neuper@48815
|
218 |
else norm ([mod_div (nth p i) lcp m] @ nrm) (i - 1) ;
|
meindl_diana@48797
|
219 |
in
|
neuper@48815
|
220 |
if List.length p = 0 then p else norm [] (List.length p - 1)
|
neuper@48815
|
221 |
end;
|
neuper@48851
|
222 |
*}
|
meindl_diana@48797
|
223 |
|
neuper@48851
|
224 |
subsection {* addition, multiplication, division *}
|
neuper@48851
|
225 |
ML {*
|
neuper@48815
|
226 |
(* scalar multiplication *)
|
neuper@48815
|
227 |
infix %*
|
neuper@48839
|
228 |
fun (p: unipoly) %* n = map (Integer.mult n) p : unipoly
|
neuper@48815
|
229 |
|
neuper@48815
|
230 |
(* scalar divison *)
|
neuper@48815
|
231 |
infix %/
|
neuper@48870
|
232 |
fun (p: unipoly) %/ n = map (swapf (curry op div2) n) p : unipoly
|
meindl_diana@48797
|
233 |
|
neuper@48815
|
234 |
(* minus of polys *)
|
meindl_diana@48797
|
235 |
infix %-%
|
neuper@48839
|
236 |
fun (p1: unipoly) %-% (p2: unipoly) =
|
neuper@48839
|
237 |
map2 (curry op -) p1 (p2 @ replicate (List.length p1 - List.length p2) 0) : unipoly
|
meindl_diana@48797
|
238 |
|
neuper@48823
|
239 |
(* dvd for univariate polynomials *)
|
meindl_diana@48797
|
240 |
infix %|%
|
neuper@48847
|
241 |
fun [d] %|% [p] = (p mod d = 0)
|
neuper@48839
|
242 |
| (ds: unipoly) %|% (ps: unipoly) =
|
neuper@48837
|
243 |
let
|
neuper@48846
|
244 |
val ds = drop_lc0_up ds;
|
neuper@48846
|
245 |
val ps = drop_lc0_up ps;
|
neuper@48839
|
246 |
val d000 = (replicate (List.length ps - List.length ds) 0) @ ds ;
|
neuper@48846
|
247 |
val quot = (lcoeff_up ps) div2 (lcoeff_up d000);
|
neuper@52077
|
248 |
val remd = drop_lc0_up (ps %-% (d000 %* quot));
|
neuper@48845
|
249 |
in
|
neuper@52077
|
250 |
if remd = [] then true else
|
neuper@52077
|
251 |
if remd <> [] andalso quot <> 0 andalso List.length ds <= List.length remd
|
neuper@52077
|
252 |
then ds %|% remd
|
neuper@48846
|
253 |
else false
|
neuper@48837
|
254 |
end
|
neuper@48851
|
255 |
*}
|
neuper@48851
|
256 |
|
neuper@48851
|
257 |
subsection {* normalisation and Landau-Mignotte bound *}
|
neuper@48851
|
258 |
ML {*
|
neuper@48855
|
259 |
(* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
|
neuper@48826
|
260 |
normalise [p_0, .., p_k] m = [q_0, .., q_k]
|
neuper@48826
|
261 |
assumes
|
neuper@48836
|
262 |
yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
|
neuper@48836
|
263 |
(where |^ x ^| means round x up to the next greater integer) *)
|
neuper@48855
|
264 |
fun centr_up (p: unipoly) (m: int) =
|
meindl_diana@48797
|
265 |
let
|
neuper@48823
|
266 |
val mi = m div2 2;
|
neuper@48823
|
267 |
val mid = if m mod mi = 0 then mi else mi + 1;
|
neuper@48823
|
268 |
fun centr m mid p_i = if mid < p_i then p_i - m else p_i;
|
neuper@48839
|
269 |
in map (centr m mid) p : unipoly end;
|
meindl_diana@48797
|
270 |
|
neuper@48836
|
271 |
(*** landau mignotte bound (Winkler p91)
|
neuper@48831
|
272 |
every coefficient of the gcd of a and b in Z[x] is bounded in absolute value by
|
neuper@48827
|
273 |
2^{min(m,n)} * gcd(a_m,b_n) * min(1/|a_m|
|
neuper@48836
|
274 |
* root(sum_{i=0}^{m}a_i^2) ,1/|b_n| * root( sum_{i=0}^{n}b_i^2 ) ***)
|
neuper@48826
|
275 |
|
neuper@48855
|
276 |
(* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
|
neuper@48826
|
277 |
sum_lmb [p_0, .., p_k] e = s
|
neuper@48826
|
278 |
assumes exp >= 0
|
neuper@48836
|
279 |
yields. p_0^e + p_1^e + ... + p_k^e *)
|
neuper@48823
|
280 |
fun sum_lmb (p: unipoly) exp = fold ((curry op +) o (Integer.pow exp)) p 0;
|
meindl_diana@48810
|
281 |
|
neuper@48855
|
282 |
(* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly \<Rightarrow> unipoly \<Rightarrow> int
|
neuper@48830
|
283 |
LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
|
neuper@48827
|
284 |
assumes True
|
neuper@48836
|
285 |
yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
|
neuper@48827
|
286 |
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
|
287 |
fun LANDAU_MIGNOTTE_bound (p1: unipoly) (p2: unipoly) =
|
neuper@48846
|
288 |
(Integer.pow (Integer.min (deg_up p1) (deg_up p2)) 2) * (abs (Integer.gcd (lcoeff_up p1) (lcoeff_up p2))) *
|
neuper@48846
|
289 |
(Int.min (abs((approx_root (sum_lmb p1 2)) div2 ~(lcoeff_up p1)),
|
neuper@48846
|
290 |
abs((approx_root (sum_lmb p2 2)) div2 ~(lcoeff_up p2))));
|
meindl_diana@48797
|
291 |
*}
|
meindl_diana@48797
|
292 |
|
neuper@48851
|
293 |
subsection {* modulo calculations on polynomials *}
|
meindl_diana@48797
|
294 |
ML{*
|
neuper@48855
|
295 |
(* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
|
neuper@48855
|
296 |
chinese_remainder_up (m1, m2) (p1, p2) = p
|
meindl_diana@48838
|
297 |
assume m1, m2 relatively prime
|
neuper@48855
|
298 |
yields p1 = p mod m1 \<and> p2 = p mod m2 *)
|
neuper@48855
|
299 |
fun chinese_remainder_up (m1, m2) (p1: unipoly, p2: unipoly) =
|
meindl_diana@48797
|
300 |
let
|
neuper@48862
|
301 |
fun chinese_rem (m1, m2) (p_1i, p_2i) = chinese_remainder (p_1i, m1) (p_2i, m2)
|
neuper@48823
|
302 |
in map (chinese_rem (m1, m2)) (p1 ~~ p2): unipoly end
|
meindl_diana@48797
|
303 |
|
neuper@48858
|
304 |
(* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
neuper@48858
|
305 |
mod_up [p1, p2, ..., pk] m = up
|
meindl_diana@48838
|
306 |
assume m > 0
|
meindl_diana@48838
|
307 |
yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
|
neuper@48858
|
308 |
infix mod_up
|
neuper@48858
|
309 |
fun (p : unipoly) mod_up m =
|
neuper@48855
|
310 |
let
|
neuper@48863
|
311 |
fun mod' m p = (curry op mod) p m
|
neuper@48863
|
312 |
in map (mod' m) p : unipoly end
|
meindl_diana@48797
|
313 |
|
meindl_diana@48819
|
314 |
(* euclidean algoritm in Z_p[x].
|
neuper@48858
|
315 |
mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
neuper@48858
|
316 |
mod_up_gcd p1 p2 m = g
|
neuper@48836
|
317 |
assumes True
|
neuper@48836
|
318 |
yields gcd ((p1 mod m), (p2 mod m)) = g *)
|
neuper@48858
|
319 |
fun mod_up_gcd (p1: unipoly) (p2: unipoly) (m: int) =
|
meindl_diana@48797
|
320 |
let
|
neuper@48858
|
321 |
val p1m = p1 mod_up m;
|
neuper@48858
|
322 |
val p2m = drop_lc0_up (p2 mod_up m);
|
neuper@48823
|
323 |
val p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
|
neuper@48846
|
324 |
val quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
|
neuper@52077
|
325 |
val remd = drop_lc0_up ((p1m %-% (p2n %* quot)) mod_up m)
|
neuper@48836
|
326 |
in
|
neuper@52077
|
327 |
if remd = [] then p2 else
|
neuper@52077
|
328 |
if length remd < List.length p2
|
neuper@52077
|
329 |
then mod_up_gcd p2 remd m
|
neuper@52077
|
330 |
else mod_up_gcd remd p2 m
|
neuper@48823
|
331 |
end;
|
neuper@48823
|
332 |
|
neuper@48855
|
333 |
(* primitive_up :: unipoly \<Rightarrow> unipoly
|
neuper@48855
|
334 |
primitive_up p = pp
|
neuper@48831
|
335 |
assumes True
|
neuper@48836
|
336 |
yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
|
neuper@48863
|
337 |
fun primitive_up ([c] : unipoly) =
|
neuper@48863
|
338 |
if c = 0 then ([0] : unipoly) else ([1] : unipoly)
|
neuper@48863
|
339 |
| primitive_up p =
|
neuper@48831
|
340 |
let
|
neuper@48831
|
341 |
val d = abs (Integer.gcds p)
|
neuper@48831
|
342 |
in
|
neuper@48831
|
343 |
if d = 1 then p else p %/ d
|
neuper@48831
|
344 |
end
|
meindl_diana@48797
|
345 |
*}
|
meindl_diana@48797
|
346 |
|
neuper@48855
|
347 |
subsection {* gcd_up, code from [1] p.93 *}
|
meindl_diana@48797
|
348 |
ML {*
|
meindl_diana@48838
|
349 |
(* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
|
meindl_diana@48838
|
350 |
try_new_prime_up p1 p2 d M P g p = new_g
|
neuper@48855
|
351 |
assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
|
neuper@48855
|
352 |
M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
|
neuper@48855
|
353 |
p1 is primitive \<and> p2 is primitive
|
neuper@48855
|
354 |
yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
|
neuper@48855
|
355 |
|
neuper@48855
|
356 |
argumentns "a b d M P g p" named according to [1] p.93 *)
|
neuper@48835
|
357 |
fun try_new_prime_up a b d M P g p =
|
neuper@48874
|
358 |
if P > M then g else
|
neuper@48833
|
359 |
let
|
neuper@48865
|
360 |
val p = p next_prime_not_dvd d
|
neuper@48865
|
361 |
val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
|
neuper@48865
|
362 |
%* (d mod p))
|
neuper@48865
|
363 |
mod_up p)
|
neuper@48865
|
364 |
p
|
neuper@48827
|
365 |
in
|
neuper@48846
|
366 |
if deg_up g_p < deg_up g
|
neuper@48827
|
367 |
then
|
neuper@48874
|
368 |
if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p
|
neuper@48827
|
369 |
else
|
neuper@48865
|
370 |
if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else
|
neuper@48827
|
371 |
let
|
neuper@48865
|
372 |
val P = P * p
|
neuper@48865
|
373 |
val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
|
neuper@48865
|
374 |
in try_new_prime_up a b d M P g p end
|
neuper@48827
|
375 |
end
|
neuper@48874
|
376 |
|
neuper@52101
|
377 |
(* iterate_CHINESE_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
|
neuper@52101
|
378 |
iterate_CHINESE_up p1 p2 d M p = gcd
|
neuper@48855
|
379 |
assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
|
neuper@48855
|
380 |
M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and>
|
neuper@48855
|
381 |
p1 is primitive \<and> p2 is primitive
|
neuper@48855
|
382 |
yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
|
neuper@48855
|
383 |
|
neuper@48855
|
384 |
argumentns "a b d M p" named according to [1] p.93 *)
|
neuper@52101
|
385 |
fun iterate_CHINESE_up a b d M p =
|
neuper@48868
|
386 |
let
|
neuper@48835
|
387 |
val p = p next_prime_not_dvd d
|
neuper@48871
|
388 |
val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p
|
neuper@48871
|
389 |
(*.. nesting of functions see above*)
|
neuper@48823
|
390 |
in
|
neuper@48874
|
391 |
if deg_up g_p = 0 then [1] else
|
neuper@48835
|
392 |
let
|
neuper@48855
|
393 |
val g = primitive_up (try_new_prime_up a b d M p g_p p)
|
neuper@48827
|
394 |
in
|
neuper@52101
|
395 |
if (g %|% a) andalso (g %|% b) then g:unipoly else iterate_CHINESE_up a b d M p
|
neuper@48827
|
396 |
end
|
neuper@48827
|
397 |
end
|
neuper@48874
|
398 |
|
neuper@48855
|
399 |
(* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
|
neuper@48855
|
400 |
gcd_up a b = c
|
meindl_diana@48838
|
401 |
assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
|
meindl_diana@48838
|
402 |
a is primitive \<and> b is primitive
|
neuper@48836
|
403 |
yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
|
neuper@48855
|
404 |
fun gcd_up a b =
|
neuper@48855
|
405 |
let val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))
|
neuper@48855
|
406 |
in
|
neuper@48855
|
407 |
if a = b then a else
|
neuper@52101
|
408 |
iterate_CHINESE_up a b d (2 * d * LANDAU_MIGNOTTE_bound a b) 1
|
neuper@48855
|
409 |
end;
|
meindl_diana@48797
|
410 |
*}
|
meindl_diana@48797
|
411 |
|
neuper@48851
|
412 |
section {* gcd for multivariate polynomials *}
|
meindl_diana@48797
|
413 |
ML {*
|
meindl_diana@48797
|
414 |
type monom = (int * int list);
|
neuper@48846
|
415 |
type poly = monom list;
|
meindl_diana@48797
|
416 |
*}
|
meindl_diana@48797
|
417 |
|
neuper@48851
|
418 |
subsection {* auxiliary functions *}
|
meindl_diana@48797
|
419 |
ML {*
|
neuper@48846
|
420 |
fun land a b = a andalso b
|
neuper@48846
|
421 |
fun all_geq a b = fold land (map2 (curry op >=) a b) true
|
neuper@48850
|
422 |
fun all_geq' i es = fold land (map ((curry op <=) i) es) true
|
neuper@48850
|
423 |
fun maxs is = fold Integer.max is (hd is);
|
neuper@48850
|
424 |
fun div2_swapped d i = i div2 d
|
neuper@48850
|
425 |
fun nth_swapped i ls = nth ls i;
|
neuper@48849
|
426 |
fun drop_last ls = (rev o tl o rev) ls
|
neuper@48849
|
427 |
fun nth_ins n i xs = (*analogous to nth_drop: why not in library.ML?*)
|
neuper@48849
|
428 |
List.take (xs, n) @ [i] @ List.drop (xs, n);
|
neuper@52086
|
429 |
fun coeffs (p : poly) = map fst p
|
neuper@48851
|
430 |
*}
|
neuper@48850
|
431 |
|
neuper@48851
|
432 |
subsection {* basic notions for multivariate polynomials *}
|
neuper@48851
|
433 |
ML {*
|
neuper@48849
|
434 |
fun lcoeff (p: poly) = (fst o hd o rev) p
|
neuper@48850
|
435 |
fun lexp (p: poly) = (snd o hd o rev) p
|
neuper@48850
|
436 |
fun lmonom (p: poly) = (hd o rev) p
|
neuper@48850
|
437 |
|
neuper@48870
|
438 |
fun zero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;
|
neuper@48850
|
439 |
|
neuper@48851
|
440 |
(* drop leading coefficients equal 0 TODO: signature?*)
|
neuper@48870
|
441 |
fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then zero_poly (length es) else p
|
neuper@48851
|
442 |
| drop_lc0 p = if lcoeff p = 0 then drop_lc0 (drop_last p) else p
|
neuper@48851
|
443 |
|
neuper@48851
|
444 |
(* gcd of coefficients WN find right location in file *)
|
neuper@48851
|
445 |
fun gcd_coeff (up: unipoly) = abs (Integer.gcds up)
|
neuper@48851
|
446 |
|
neuper@48851
|
447 |
(* given p = 3x+4z = [(3,[1,0]),(4,[0,1])] with x<z
|
neuper@48851
|
448 |
add_variable p 2 => [(3,[1,0,0]),(4,[0,0,1])] with x<y<z *)
|
neuper@48870
|
449 |
fun add_variable pos (p: poly) = map (apsnd (nth_ins pos 0)) p : poly
|
neuper@48851
|
450 |
*}
|
neuper@48851
|
451 |
|
neuper@48851
|
452 |
subsection {* ordering and other auxiliary funcions *}
|
neuper@48851
|
453 |
ML {*
|
neuper@48857
|
454 |
(* monomial order: lexicographic order on exponents *)
|
neuper@48850
|
455 |
fun lex_ord ((_, a): monom) ((_, b): monom) =
|
neuper@48857
|
456 |
let val d = drop_tl_zeros (a %-% b)
|
neuper@48850
|
457 |
in
|
neuper@48866
|
458 |
if d = [] then false else last d > 0
|
neuper@48850
|
459 |
end
|
neuper@48849
|
460 |
fun lex_ord' ((_, a): monom, (_, b): monom) =
|
neuper@48857
|
461 |
let val d = drop_tl_zeros (a %-% b)
|
neuper@48849
|
462 |
in
|
neuper@48857
|
463 |
if d = [] then EQUAL
|
neuper@48866
|
464 |
else if last d > 0 then GREATER else LESS
|
neuper@48849
|
465 |
end
|
neuper@48850
|
466 |
|
neuper@48848
|
467 |
(* add monomials in ordered polynomal *)
|
neuper@48850
|
468 |
fun add_monoms (p: poly) =
|
neuper@48846
|
469 |
let
|
neuper@48870
|
470 |
fun add [(0, es)] [] = zero_poly (length es)
|
neuper@48847
|
471 |
| add [(0, _)] (new: poly) = new
|
neuper@48847
|
472 |
| add [m] (new: poly) = new @ [m]
|
neuper@48847
|
473 |
| add ((c1, es1) :: (c2, es2) :: ms) (new: poly) =
|
neuper@48846
|
474 |
if es1 = es2
|
neuper@48847
|
475 |
then add ([(c1 + c2, es1)] @ ms) new
|
neuper@48846
|
476 |
else
|
neuper@48846
|
477 |
if c1 = 0
|
neuper@48847
|
478 |
then add ((c2, es2) :: ms) new
|
neuper@48847
|
479 |
else add ((c2, es2) :: ms) (new @ [(c1, es1)])
|
neuper@48850
|
480 |
in add p [] : poly end
|
neuper@48850
|
481 |
|
neuper@48848
|
482 |
(* brings the polynom in correct order and adds monomials *)
|
neuper@48851
|
483 |
fun order (p: poly) = (add_monoms o (sort lex_ord')) p;
|
neuper@48850
|
484 |
|
neuper@48850
|
485 |
(* largest exponent of variable at location var *)
|
neuper@48851
|
486 |
fun max_deg_var (p: poly) var = (maxs o (map ((nth_swapped var) o snd))) p;
|
neuper@48851
|
487 |
*}
|
neuper@48850
|
488 |
|
neuper@48851
|
489 |
subsection {* evaluation, translation uni--multivariate, etc *}
|
neuper@48851
|
490 |
ML {*
|
neuper@48868
|
491 |
(* evaluate a polynomial in a certain variable and remove this var from the exp-list *)
|
neuper@48870
|
492 |
fun eval_monom var value ((c, es): monom) =
|
neuper@48870
|
493 |
(c * (Integer.pow (nth es var) value), nth_drop var es) : monom
|
neuper@48870
|
494 |
fun eval (p: poly) var value = order (map (eval_monom var value) p)
|
neuper@48851
|
495 |
|
neuper@48870
|
496 |
(* transform a multivariate to a univariate polynomial TODO map?*)
|
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@48870
|
508 |
(* transform a univariate to a multivariate polynomial TODO map *)
|
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@48870
|
516 |
e.g new z: (3x + 2y)*(z - 4) TODO mp %%*%% (uni_to_multi up) !TODO equal es ! *)
|
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 |
fun mult_monom ((c1, es1): monom) ((c2, es2): monom) =
|
neuper@48853
|
585 |
(c1 * c2, map2 Integer.add es1 es2): monom;
|
neuper@48853
|
586 |
fun mult_monoms (ms: poly) m = map (mult_monom m) ms : poly
|
neuper@48853
|
587 |
fun p1 %%*%% (p2 : poly) = order (fold (curry op @) (map (mult_monoms p1) p2) [])
|
neuper@48853
|
588 |
|
meindl_diana@48797
|
589 |
infix %%|%%
|
neuper@48851
|
590 |
(* the same algorithms as done by hand; TODO fold / map ? *)
|
neuper@48850
|
591 |
fun ([(c1, es1)]: poly) %%|%% ([(c2, es2)]: poly) = all_geq es2 es1 andalso c2 mod c1 = 0
|
neuper@48850
|
592 |
| _ %%|%% [(0,_)] = true
|
neuper@48850
|
593 |
| p1 %%|%% p2 =
|
neuper@48850
|
594 |
if [lmonom p1] %%|%% [lmonom p2]
|
neuper@48850
|
595 |
then
|
neuper@48850
|
596 |
p1 %%|%%
|
neuper@48850
|
597 |
(p2 %%-%%
|
neuper@48850
|
598 |
([((lcoeff p2) div2 (lcoeff p1), (lexp p2) %-% (lexp p1))]
|
neuper@48850
|
599 |
%%*%% p1))
|
neuper@48850
|
600 |
else false
|
neuper@48853
|
601 |
|
neuper@48851
|
602 |
(* %%/%% :: poly \<Rightarrow> poly \<Rightarrow> (poly \<times> poly)
|
neuper@52063
|
603 |
p1 %%/%% p2 = (q, r)
|
neuper@52063
|
604 |
assumes p2 \<noteq> [] \<and> deg p1 \<ge> deg p2
|
neuper@52063
|
605 |
yields THE q r. p1 = q *\<^isub>p p2 +\<^isub>p r *)
|
neuper@48851
|
606 |
infix %%/%%
|
neuper@48851
|
607 |
fun (p1: poly) %%/%% (p2: poly) =
|
neuper@48851
|
608 |
let
|
neuper@48851
|
609 |
fun div_up p1 p2 quot_old =
|
neuper@48851
|
610 |
let
|
neuper@48851
|
611 |
val p1 as (c, es)::_ = drop_lc0 (order p1);
|
neuper@48851
|
612 |
val p2 = drop_lc0 (order p2);
|
neuper@48870
|
613 |
val [c] = zero_poly (length es);
|
neuper@48851
|
614 |
val d = (replicate (List.length p1 - List.length p2) c) @ p2 ;
|
neuper@48851
|
615 |
val quot = [(lcoeff p1 div2 lcoeff d, lexp p1 %-% lexp p2)];
|
neuper@52077
|
616 |
val remd = drop_lc0 (p1 %%-%% (d %%*%% quot));
|
neuper@48851
|
617 |
in
|
neuper@52077
|
618 |
if remd <> [c] andalso all_geq' 0 (lexp remd %-% lexp p2)
|
neuper@52077
|
619 |
then div_up remd p2 (quot @ quot_old)
|
neuper@52077
|
620 |
else (quot @ quot_old : poly, remd)
|
neuper@48851
|
621 |
end
|
neuper@48851
|
622 |
in div_up p1 p2 [] end
|
meindl_diana@48797
|
623 |
*}
|
neuper@48848
|
624 |
|
neuper@48846
|
625 |
subsection {* Newton interpolation *}
|
meindl_diana@48797
|
626 |
ML{* (* first step *)
|
neuper@48870
|
627 |
fun newton_first (x: int list) (f: poly list) (pos: int) =
|
neuper@48851
|
628 |
let
|
neuper@48851
|
629 |
val new_vals = [(hd x) * ~1, 1];
|
neuper@48870
|
630 |
val p = (add_variable pos (hd f)) %%+%%
|
neuper@48851
|
631 |
((mult_with_new_var
|
neuper@48870
|
632 |
(((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x))) new_vals pos);
|
neuper@48851
|
633 |
val steps = [((nth f 1) %%-%% (hd f)) %%/ ((nth x 1) - (hd x))]
|
neuper@48851
|
634 |
in (p, new_vals, steps) end
|
meindl_diana@48797
|
635 |
|
neuper@48851
|
636 |
fun next_step ([]: poly list) (new_steps: poly list) (_: int list) = new_steps
|
neuper@48851
|
637 |
| next_step steps new_steps x =
|
neuper@48851
|
638 |
next_step (tl steps)
|
neuper@48851
|
639 |
(new_steps @
|
neuper@48866
|
640 |
[((last new_steps) %%-%% (hd steps)) %%/ ((last x) - (nth x (length x - 3)))])
|
neuper@48851
|
641 |
(nth_drop (length x -2) x)
|
neuper@48851
|
642 |
|
neuper@48870
|
643 |
fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) pos =
|
meindl_diana@48797
|
644 |
if length x = 2
|
neuper@48851
|
645 |
then
|
neuper@48870
|
646 |
let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] pos
|
neuper@48851
|
647 |
in (p', new_vals, steps) end
|
neuper@48851
|
648 |
else
|
neuper@48850
|
649 |
let
|
neuper@48851
|
650 |
val new_vals =
|
neuper@48851
|
651 |
multi_to_uni ((uni_to_multi up) %%*%% (uni_to_multi [(nth x (length x -2) ) * ~1, 1]));
|
neuper@48851
|
652 |
val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/
|
neuper@48866
|
653 |
((last x) - (nth x (length x - 2)))];
|
neuper@48851
|
654 |
val steps = next_step steps new_steps x;
|
neuper@48870
|
655 |
val p' = p %%+%% (mult_with_new_var (last steps) new_vals pos);
|
neuper@48851
|
656 |
in (order p', new_vals, steps) end
|
meindl_diana@48841
|
657 |
*}
|
meindl_diana@48797
|
658 |
|
neuper@48851
|
659 |
subsection {* gcd_poly algorithm, code for [1] p.95 *}
|
meindl_diana@48797
|
660 |
ML {*
|
neuper@52099
|
661 |
fun gcd_monom (c1, es1) (c2, es2) = (Integer.gcd c1 c2, map2 Integer.min es1 es2)
|
neuper@52099
|
662 |
|
neuper@48868
|
663 |
(* naming of n, M, m, r, ... follows Winkler p. 95
|
neuper@52099
|
664 |
assumes: n = length es1 = length es2
|
neuper@48868
|
665 |
r: *)
|
neuper@52099
|
666 |
fun gcd_poly' a [monom as (c, es)] _ _ = [fold gcd_monom a monom]
|
neuper@52099
|
667 |
| gcd_poly' [monom as (c, es)] b _ _ = [fold gcd_monom b monom]
|
neuper@52098
|
668 |
| gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r =
|
neuper@52098
|
669 |
if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else
|
neuper@52098
|
670 |
if n > 1
|
neuper@52098
|
671 |
then
|
neuper@52098
|
672 |
let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
|
neuper@52098
|
673 |
in try_new_r a b n M r [] [] end
|
neuper@52098
|
674 |
else
|
neuper@52098
|
675 |
let val (a', b') = (multi_to_uni a, multi_to_uni b)
|
neuper@52098
|
676 |
in
|
neuper@52098
|
677 |
(* a b are made primitive and this is compensated by multiplying with gcd of divisors *)
|
neuper@52098
|
678 |
uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %*
|
neuper@52098
|
679 |
(abs (Integer.gcd (gcd_coeff a') (gcd_coeff b'))))
|
neuper@52098
|
680 |
end
|
neuper@48867
|
681 |
|
neuper@52101
|
682 |
(* try_new_r :: poly -> poly -> int -> int -> int -> int list -> poly list -> poly
|
neuper@52101
|
683 |
try_new_r a b n M r r_list steps =
|
neuper@52099
|
684 |
assumes length a > 1 \<and> length b > 1
|
neuper@52101
|
685 |
yields TODO
|
neuper@52099
|
686 |
*)
|
neuper@48851
|
687 |
and try_new_r a b n M r r_list steps =
|
meindl_diana@48816
|
688 |
let
|
meindl_diana@48810
|
689 |
val r = find_new_r a b r;
|
meindl_diana@48810
|
690 |
val r_list = r_list @ [r];
|
neuper@48855
|
691 |
val g_r = gcd_poly' (order (eval a (n - 2) r))
|
neuper@48855
|
692 |
(order (eval b (n - 2) r)) (n - 1) 0
|
meindl_diana@48810
|
693 |
val steps = steps @ [g_r];
|
neuper@52101
|
694 |
in iterate_CHINESE a b n M 1 r r_list steps g_r (zero_poly n) [1] end
|
neuper@52101
|
695 |
|
neuper@52101
|
696 |
(* iterate_CHINESE :: poly -> poly -> int -> int -> int -> int -> int list ->
|
neuper@52101
|
697 |
| poly list -> poly -> poly -> int list -> poly
|
neuper@52101
|
698 |
iterate_CHINESE a | b n M m r r_list
|
neuper@52101
|
699 |
steps g g_n mult =
|
neuper@52099
|
700 |
assumes length a > 1 \<and> length b > 1
|
neuper@52101
|
701 |
yields TODO
|
neuper@52099
|
702 |
*)
|
neuper@52101
|
703 |
and iterate_CHINESE a b n M m r r_list steps g g_n mult =
|
neuper@48869
|
704 |
if m > M then if g_n %%|%% a andalso g_n %%|%% b then g_n else
|
neuper@52101
|
705 |
try_new_r a b n M r r_list steps
|
meindl_diana@48810
|
706 |
else
|
meindl_diana@48810
|
707 |
let
|
meindl_diana@48810
|
708 |
val r = find_new_r a b r;
|
meindl_diana@48810
|
709 |
val r_list = r_list @ [r];
|
neuper@48855
|
710 |
val g_r = gcd_poly' (order (eval a (n - 2) r))
|
neuper@48855
|
711 |
(order (eval b (n - 2) r)) (n - 1) 0
|
meindl_diana@48810
|
712 |
in
|
neuper@48849
|
713 |
if lex_ord (lmonom g) (lmonom g_r)
|
neuper@52101
|
714 |
then iterate_CHINESE a b n M 1 r r_list steps g g_n mult
|
neuper@52101
|
715 |
else if lexp g = lexp g_r
|
neuper@48867
|
716 |
then
|
neuper@48867
|
717 |
let val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
|
neuper@48867
|
718 |
in
|
neuper@52101
|
719 |
if nth steps (length steps - 1) = zero_poly (n - 1)
|
neuper@52101
|
720 |
then iterate_CHINESE a b n M (M + 1) r r_list steps g_r g_n new
|
neuper@52101
|
721 |
else iterate_CHINESE a b n M (m + 1) r r_list steps g_r g_n new
|
neuper@48867
|
722 |
end
|
neuper@48867
|
723 |
else (* \<not> lex_ord (lmonom g) (lmonom g_r) *)
|
neuper@52101
|
724 |
iterate_CHINESE a b n M (m + 1) r r_list steps g g_n mult
|
neuper@48868
|
725 |
end
|
neuper@48853
|
726 |
|
neuper@52101
|
727 |
fun majority_of_coefficients_is_negative_in a b c =
|
neuper@52086
|
728 |
let val cs = (coeffs a) @ (coeffs b) @ (coeffs c)
|
neuper@52086
|
729 |
in length (filter (curry op < 0) cs) < length cs div 2 end
|
neuper@52086
|
730 |
|
neuper@48851
|
731 |
(* gcd_poly :: poly \<Rightarrow> poly \<Rightarrow> poly
|
neuper@48851
|
732 |
gcd_poly a b = ((a', b'), c)
|
neuper@48851
|
733 |
assumes a \<noteq> [] \<and> b \<noteq> [] \<and>
|
neuper@52086
|
734 |
yields a = a' *\<^isub>p c \<and> b = b' *\<^isub>p c \<and> (\<forall>c'. (c' dvd\<^isub>p a \<and> c' dvd\<^isub>p b) \<longrightarrow> c' \<le>\<^isub>p c) \<and>
|
neuper@52086
|
735 |
\<not> gcds a' < 0 \<and> gcds a' < 0 *)
|
neuper@48851
|
736 |
fun gcd_poly (a as (_, es)::_ : poly) b =
|
neuper@48851
|
737 |
let val c = gcd_poly' a b (length es) 0
|
neuper@48851
|
738 |
val (a': poly, _) = a %%/%% c
|
neuper@48851
|
739 |
val (b': poly, _) = b %%/%% c
|
neuper@52086
|
740 |
in
|
neuper@52101
|
741 |
if majority_of_coefficients_is_negative_in a' b' c
|
neuper@52086
|
742 |
then ((a' %%* ~1, b' %%* ~1), c %%* ~1) (*makes results look nicer*)
|
neuper@52086
|
743 |
else ((a', b'), c)
|
neuper@52086
|
744 |
end
|
meindl_diana@48797
|
745 |
*}
|
meindl_diana@48797
|
746 |
|
neuper@52063
|
747 |
section {* example from the last mail to Tobias Nipkow *}
|
neuper@48851
|
748 |
ML {* (* test, see text below *)
|
neuper@52063
|
749 |
val a = [(1,[2, 0]), (~1,[0, 2])]; (* x^2 - y^2 *)
|
neuper@52063
|
750 |
val b = [(1,[2, 0]), (~1,[1, 1])]; (* x^2 - x*y *)
|
neuper@48851
|
751 |
val ((a', b'), c) = gcd_poly a b;
|
neuper@48851
|
752 |
|
neuper@48851
|
753 |
a' = [(~1, [1, 0]), (~1, [0, 1])];
|
neuper@48851
|
754 |
b' = [(~1, [1, 0])];
|
neuper@48867
|
755 |
c = [(~1, [1, 0]), (1, [0, 1])];
|
neuper@48867
|
756 |
(* checking the postcondition: *)
|
neuper@48867
|
757 |
a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)
|
neuper@48867
|
758 |
(* \<forall>c'. (c' dvd\<^isub>p a \<and> c' dvd\<^isub>p b) \<longrightarrow> c' \<le>\<^isub>p c) can NOT be checked this way, of course *)
|
neuper@48851
|
759 |
*}
|
neuper@48851
|
760 |
|
neuper@48851
|
761 |
text {* last mail to Tobias Nipkow:
|
neuper@48851
|
762 |
Eine erste Idee, wie die Integration der Diplomarbeit für einen Benutzer
|
neuper@48851
|
763 |
von Isabelle aussehen könnte, wäre zum Beispiel im
|
neuper@48851
|
764 |
|
neuper@48851
|
765 |
lemma cancel:
|
neuper@48851
|
766 |
assumes asm3: "x^2 - x*y \<noteq> 0" and asm4: "x \<noteq> 0"
|
neuper@48851
|
767 |
shows "(x^2 - y^2) / (x^2 - x*y) = (x + y) / (x::real)"
|
neuper@48851
|
768 |
apply (insert asm3 asm4)
|
neuper@48851
|
769 |
apply simp
|
neuper@48851
|
770 |
sorry
|
neuper@48851
|
771 |
|
neuper@48851
|
772 |
die Assumptions
|
neuper@48851
|
773 |
|
neuper@48851
|
774 |
asm1: "(x^2 - y^2) = (x + y) * (x - y)" and asm2: "x^2 - x*y = x * (x - y)"
|
neuper@48851
|
775 |
|
neuper@48851
|
776 |
im Hintergrund automatisch zu erzeugen (mit der Garantie, dass "(x - y)" der GCD ist)
|
neuper@48851
|
777 |
und sie dem Simplifier (für die Rule nonzero_mult_divide_mult_cancel_right) zur
|
neuper@48851
|
778 |
Verfügung zu stellen, sodass anstelle von "sorry" einfach "done" stehen kann.
|
neuper@48851
|
779 |
Und weiters wäre eventuell asm3 zu "x - y \<noteq> 0" zu vereinfachen, eine
|
neuper@48851
|
780 |
Rewriteorder zum Herstellen einer Normalform festzulegen, etc.
|
neuper@48851
|
781 |
*}
|
neuper@52063
|
782 |
|
neuper@52063
|
783 |
section {* Preparations for Lucas-Interpretation *}
|
neuper@52063
|
784 |
text {*
|
neuper@52063
|
785 |
Lucas-Interpretation shall make "gcd_poly" transparent to the user,
|
neuper@52063
|
786 |
such that step-wise interpretation gives some intuition of the algorithm
|
neuper@52063
|
787 |
and opportunities for interactive learning by trial and error.
|
neuper@52063
|
788 |
The vision is to provide a learning environment similar to a good
|
neuper@52063
|
789 |
chess program, where the matter (chess) is selfcontained and
|
neuper@52063
|
790 |
learning happens by interaction with the system.
|
neuper@52063
|
791 |
|
neuper@52063
|
792 |
GDC_Poly gives raise to a major case study on the above aims.
|
neuper@52063
|
793 |
Below pedagogical questions come first and respective technical
|
neuper@52063
|
794 |
requirements subsequently.
|
neuper@52063
|
795 |
*}
|
neuper@52063
|
796 |
|
neuper@52077
|
797 |
subsection {* Questions of learners *}
|
neuper@52063
|
798 |
text {*
|
neuper@52063
|
799 |
The earlier the phase in appropaching a new topic for a student, the more
|
neuper@52063
|
800 |
divergent the questions are. So here are just some which suggest themselves.
|
neuper@52063
|
801 |
*}
|
neuper@52063
|
802 |
subsubsection {* Why does naive transfer from Z to polynomials not work? *}
|
neuper@52063
|
803 |
text {*
|
neuper@52063
|
804 |
Given the Eucliden Algorithm in Z ...
|
neuper@52063
|
805 |
*}
|
neuper@52063
|
806 |
ML {*
|
neuper@52063
|
807 |
fun euclid_int a 0 = a
|
neuper@52063
|
808 |
| euclid_int a b =
|
neuper@52063
|
809 |
if a < b then euclid_int b a else
|
neuper@52063
|
810 |
let
|
neuper@52063
|
811 |
val (q, r) = (a div b, a mod b)
|
neuper@52063
|
812 |
val _ = writeln (PolyML.makestring a ^ " = " ^ PolyML.makestring q ^
|
neuper@52063
|
813 |
" * " ^ PolyML.makestring b ^ " + " ^ PolyML.makestring r)
|
neuper@52063
|
814 |
in euclid_int b r : int end;
|
neuper@52063
|
815 |
euclid_int 192 162; (* ..stepwise execution is essential for comprehension *)
|
neuper@52063
|
816 |
*}
|
neuper@52063
|
817 |
text {* ... let it naively transfer to the polynomial ring Z[x]: *}
|
neuper@52063
|
818 |
ML {*
|
neuper@52063
|
819 |
fun deg p = max_deg_var (p: poly) 0 (*TODO?*)
|
neuper@52063
|
820 |
fun EUCLID_naive a [(0, [0])] = a (*TODO: for Z[x,y] we would need [(0, [0, 0])] here*)
|
neuper@52063
|
821 |
| EUCLID_naive a b =
|
neuper@52063
|
822 |
if deg a < deg b then EUCLID_naive b a else
|
neuper@52063
|
823 |
let
|
neuper@52063
|
824 |
val (q, r) = a %%/%% b
|
neuper@52063
|
825 |
val _ = writeln (PolyML.makestring a ^ " == " ^ PolyML.makestring q ^
|
neuper@52063
|
826 |
" ** " ^ PolyML.makestring b ^ " ++ " ^ PolyML.makestring r)
|
neuper@52063
|
827 |
in EUCLID_naive b r : poly end;
|
neuper@52063
|
828 |
*}
|
neuper@52063
|
829 |
ML {*
|
neuper@52063
|
830 |
(*The example used in the inital mail to Tobias Nipkow, simplified to univariate*)
|
neuper@52063
|
831 |
"--------------------------------";
|
neuper@52063
|
832 |
val a = uni_to_multi [0,~1,1]; (* -x + x^2 = x *(-1 + x) *)
|
neuper@52063
|
833 |
val b = uni_to_multi [~1,0,1]; (* -1 + x^2 = (1+x)*(-1 + x) *)
|
neuper@52063
|
834 |
"--------------------------------";
|
neuper@52063
|
835 |
EUCLID_naive a b = [(1, [0]), (~1, [1])]; (* 1 - x .. is acceptable*)
|
neuper@52063
|
836 |
*}
|
neuper@52063
|
837 |
ML {*
|
neuper@52063
|
838 |
"--------------------------------";
|
neuper@52063
|
839 |
val a = uni_to_multi [~1,0,0,1]; (* -1 + x^3 = (1+x+x^2)*(-1 + x) *)
|
neuper@52063
|
840 |
val b = uni_to_multi [0,~1,1]; (* -x + x^2 = x *(-1 + x) *)
|
neuper@52063
|
841 |
"--------------------------------";
|
neuper@52063
|
842 |
EUCLID_naive a b = [(~1, [0]), (1, [1])]; (*-1 + x *)
|
neuper@52063
|
843 |
*}
|
neuper@52063
|
844 |
text {*
|
neuper@52063
|
845 |
EUCLID_naive, however, does not work in general, because division is
|
neuper@52063
|
846 |
not as general as in Z:
|
neuper@52063
|
847 |
*}
|
neuper@52063
|
848 |
ML {*
|
neuper@52063
|
849 |
"--------------------------------";
|
neuper@52063
|
850 |
val a = uni_to_multi [0,~1,0,0,1]; (* -x + x^4 = x*(1+x+x^2)*(-1 + x) *)
|
neuper@52063
|
851 |
val b = uni_to_multi [~2,2]; (* -x + x^2 = 2 *(-1 + x) *)
|
neuper@52063
|
852 |
"--------------------------------";
|
neuper@52063
|
853 |
(*EUCLID_naive a b = [(~1, [0]), (1, [1])] happens to loop*)
|
neuper@52063
|
854 |
*}
|
neuper@52063
|
855 |
text {*
|
neuper@52063
|
856 |
(* (x^4 - x) : (2*x - 2) = 1/2 * x^4 ...
|
neuper@52077
|
857 |
This division, equal to above division, is impossible in Z[x], rather, it requires Q[x].
|
neuper@52077
|
858 |
This problem of division with integer coefficients proceeds down from the least
|
neuper@52077
|
859 |
coefficient to the other coefficients.
|
neuper@52063
|
860 |
Nevertheless, below we are going to generalise division such that the
|
neuper@52063
|
861 |
Euclidean Algorithm can work on Z[x].
|
neuper@52063
|
862 |
*}
|
neuper@52063
|
863 |
|
neuper@52063
|
864 |
subsubsection {* Can division of polynomials be adapted to Euclid? *}
|
neuper@52063
|
865 |
text {*
|
neuper@52063
|
866 |
In principle it is simplest to go from Z[x] to Q[x], polynamials over rational numbers
|
neuper@52063
|
867 |
and thus, for instance, allow
|
neuper@52063
|
868 |
(x^4 - x) : (2*x - 2) = 1/2 * x^4 ...
|
neuper@52077
|
869 |
However, Q is computationally less efficient and the terms are more complicated.
|
neuper@52077
|
870 |
Fortunately, [1].p.83/84 tells, that division in Z[x] can be generalised for the
|
neuper@52063
|
871 |
Euclidean Algorithm.
|
neuper@52063
|
872 |
In order to comprehend stepwise execution of the division algorithm below
|
neuper@52063
|
873 |
(which is exactly as done by hand) we adapt to our representation of (univariate)
|
neuper@52063
|
874 |
polynomials as lists and take an instructuve example:
|
neuper@52063
|
875 |
|
neuper@52063
|
876 |
(2*x^4 + 2*x^3 + 2*x^2 + 2*x + 2) : (3*x + 4) = 2/3 * x^2 ...
|
neuper@52063
|
877 |
|
neuper@52063
|
878 |
is translated to a reversed arrangement of monomials
|
neuper@52063
|
879 |
|
neuper@52063
|
880 |
(2 + 2*x + 2*x^2 + 2*x^3 + 2*x^4) : (1 + 2*x + 3*x^2) = ...
|
neuper@52063
|
881 |
|
neuper@52063
|
882 |
and further to lists:
|
neuper@52063
|
883 |
|
neuper@52063
|
884 |
[ 2, 2, 2, 2, 2] : [1, 2, 3] =
|
neuper@52063
|
885 |
3*[ 2, 2, 2, 2, 2] *3
|
neuper@52063
|
886 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:
|
neuper@52063
|
887 |
[ 6, 6, 6, 6, 6] :
|
neuper@52063
|
888 |
--[ 0, 0, 2, 4, 6] <--........ * [ 2]
|
neuper@52063
|
889 |
------------------------ :
|
neuper@52063
|
890 |
[ 6, 6, 4, 2, 0] = [ 2]*3
|
neuper@52063
|
891 |
3*[ 6, 6, 4, 2] *3
|
neuper@52063
|
892 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:
|
neuper@52063
|
893 |
[ 18, 18, 12, 6] :
|
neuper@52063
|
894 |
--[ 0, 2, 4, 6] <-------........--- * [2]
|
neuper@52063
|
895 |
------------------------ :
|
neuper@52063
|
896 |
[ 18, 16, 8, 0] = [ 6, 2]
|
neuper@52063
|
897 |
3*[ 18, 16, 8] *3
|
neuper@52063
|
898 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:
|
neuper@52063
|
899 |
[ 54, 48, 24] :
|
neuper@52063
|
900 |
--[ 8, 16, 24] <------------........-------- * [8]
|
neuper@52063
|
901 |
------------------------ :
|
neuper@52063
|
902 |
[46, 32, 0] = [18, 6, 8]
|
neuper@52063
|
903 |
=======================================
|
neuper@52077
|
904 |
27*[ 2, 2, 2, 2, 2] : [1, 2, 3] = [18, 6, 8], remd [46, 32]
|
neuper@52063
|
905 |
|
neuper@52063
|
906 |
We implement this algorithm below and check it with these and other values.
|
neuper@52063
|
907 |
|
neuper@52077
|
908 |
First we need some auxiliary functions.
|
neuper@52063
|
909 |
*}
|
neuper@52063
|
910 |
ML {* (* auxiliary functions *)
|
neuper@52077
|
911 |
fun for_quot_regarding p1 p2 p1' quot remd =
|
neuper@52077
|
912 |
let
|
neuper@52077
|
913 |
val zeros = deg_up p1' - deg_up remd - 1(*length [q]*)
|
neuper@52077
|
914 |
val max_qot_deg = deg_up p1 - deg_up p2
|
neuper@52077
|
915 |
in
|
neuper@52077
|
916 |
if zeros + 1(*length [q]*) + deg_up quot > max_qot_deg (*possible in last step of div*)
|
neuper@52077
|
917 |
then max_qot_deg - deg_up quot - 1(*length [q]*) else zeros
|
neuper@52077
|
918 |
end
|
neuper@52063
|
919 |
|
neuper@52077
|
920 |
(* multiply polynomial p such that it has degree d with d = deg p*)
|
neuper@52077
|
921 |
fun mult_to_deg d p =
|
neuper@52077
|
922 |
let val d' = deg_up p
|
neuper@52077
|
923 |
in if d >= d' then (replicate (d - d') 0) @ p : unipoly else p end;
|
neuper@52063
|
924 |
|
neuper@52077
|
925 |
(* find a factor n and a quotient q such that n > 0 \<and> n*c1 = q*c2 *)
|
neuper@52063
|
926 |
fun fact_quot c1 c2 =
|
neuper@52077
|
927 |
let
|
neuper@52077
|
928 |
val d = abs (Integer.gcd c1 c2)
|
neuper@52077
|
929 |
val n = c2 div2 d
|
neuper@52077
|
930 |
val q = c1 div2 d
|
neuper@52077
|
931 |
in if n > 0 then (n, q) else (~1 * n, ~1 * q) end;
|
neuper@52063
|
932 |
|
neuper@52063
|
933 |
infix %+%
|
neuper@52063
|
934 |
fun p1 %+% p2 =
|
neuper@52063
|
935 |
if length p1 > length p2
|
neuper@52063
|
936 |
then map2 (curry op +) p1 (p2 @ (replicate (List.length p1 - List.length p2) 0))
|
neuper@52063
|
937 |
else map2 (curry op +) (p1 @ (replicate (List.length p2 - List.length p1) 0)) p2;
|
neuper@52063
|
938 |
|
neuper@52063
|
939 |
infix %*%
|
neuper@52063
|
940 |
fun p1 %*% p2 =multi_to_uni (uni_to_multi p1 %%*%% uni_to_multi p2);
|
neuper@52063
|
941 |
([1,2,1] %*% [1,3,3,1]) = [1, 5, 10, 10, 5, 1];
|
neuper@52063
|
942 |
|
neuper@52077
|
943 |
replicate 3 0 = [0, 0, 0];
|
neuper@52077
|
944 |
replicate 0 0 = [];
|
neuper@52077
|
945 |
|
neuper@52063
|
946 |
val trace_div = Unsynchronized.ref true;
|
neuper@52077
|
947 |
val trace_div_invariant = Unsynchronized.ref false; (*.. true expects !trace_div = false *)
|
neuper@52077
|
948 |
fun div_invariant2 p1 p2 n ns zeros q quot remd =
|
neuper@52077
|
949 |
(PolyML.makestring p1 ^ " * " ^ PolyML.makestring (n * ns) ^ " = " ^
|
neuper@52077
|
950 |
PolyML.makestring (mult_to_deg (deg_up p1 - deg_up p2)
|
neuper@52077
|
951 |
((zeros @ [q] @ (quot %* n)))) ^
|
neuper@52077
|
952 |
" ** " ^ PolyML.makestring p2 ^ " ++ " ^
|
neuper@52077
|
953 |
PolyML.makestring ((rev o (mult_to_deg (deg_up p1)) o rev) remd))
|
neuper@52077
|
954 |
fun writeln_trace p1 p1' p2 quot n q (zeros:int) remd (ns: int) =
|
neuper@52077
|
955 |
(if (! trace_div) then
|
neuper@52077
|
956 |
( writeln (" div " ^ PolyML.makestring p1' ^ " * " ^ PolyML.makestring n);
|
neuper@52077
|
957 |
writeln " div ~~~~~~~~~~~~~~~~~~~~~~~~~~~";
|
neuper@52077
|
958 |
writeln (" div " ^ PolyML.makestring (p1' %* n));
|
neuper@52077
|
959 |
writeln (" div-- " ^ PolyML.makestring (mult_to_deg (deg_up p1') p2 %* q) ^
|
neuper@52077
|
960 |
" <--- " ^ PolyML.makestring p2 ^ " * " ^ PolyML.makestring q);
|
neuper@52077
|
961 |
writeln " div ---------------------------";
|
neuper@52077
|
962 |
writeln (" div " ^ PolyML.makestring ((p1' %* n)%-%(mult_to_deg (deg_up p1') p2 %* q)) ^
|
neuper@52077
|
963 |
" quot = " ^
|
neuper@52077
|
964 |
PolyML.makestring (replicate zeros 0) ^ " @ [" ^
|
neuper@52077
|
965 |
PolyML.makestring q ^ "] @ (" ^ PolyML.makestring n ^
|
neuper@52077
|
966 |
" * " ^ PolyML.makestring quot ^ ")" );
|
neuper@52077
|
967 |
if (p1' %* n) = ((mult_to_deg (deg_up p1') (p2 %* q)) %+% remd) then ()
|
neuper@52077
|
968 |
else error ("invariant 1 does not hold: " ^ PolyML.makestring (p1' %* n) ^ " = " ^
|
neuper@52077
|
969 |
PolyML.makestring (mult_to_deg (deg_up p1') p2) ^ " * " ^ PolyML.makestring q ^ " ++ " ^
|
neuper@52077
|
970 |
PolyML.makestring (remd @ []))
|
neuper@52077
|
971 |
) else ();
|
neuper@52077
|
972 |
(*invariant 2: initial p1 *\<^isub>s ns = quot *\<^isub>p p2 +\<^isub>p remd *)
|
neuper@52077
|
973 |
if (p1 %* (n * ns)) =
|
neuper@52077
|
974 |
((mult_to_deg (deg_up p1 - deg_up p2)
|
neuper@52077
|
975 |
((replicate zeros 0) @ [q] @ (quot %* n)) %*% p2 %+% remd))
|
neuper@52077
|
976 |
then
|
neuper@52077
|
977 |
if (! trace_div_invariant)
|
neuper@52077
|
978 |
then writeln (" div " ^ div_invariant2 p1 p2 n ns (replicate zeros 0) q quot remd)
|
neuper@52077
|
979 |
else ()
|
neuper@52077
|
980 |
else error ("invariant 2 does not hold: " ^ div_invariant2 p1 p2 n ns (replicate zeros 0) q quot remd)
|
neuper@52077
|
981 |
)
|
neuper@52077
|
982 |
|
neuper@52063
|
983 |
(* generalized division: divide in Z[x] with multiplication of the dividend
|
neuper@52063
|
984 |
such that the quotient is again in Z[x], i.e. has integer coefficients.
|
neuper@52063
|
985 |
fact_div :: unipoly \<Rightarrow> unipoly \<Rightarrow> (unipoly, unipoly, nat)
|
neuper@52063
|
986 |
fact_div p1 p2 = (n, q, r)
|
neuper@52063
|
987 |
assumes p2 \<noteq> [0] \<and> deg p1 \<ge> deg p2
|
neuper@52077
|
988 |
yields p1 *\<^isub>s n = q *\<^isub>p p2 +\<^isub>p r
|
neuper@52077
|
989 |
|
neuper@52077
|
990 |
note 1: p1' changes with iterations, p1 remains equal -- only required for invariant 2
|
neuper@52077
|
991 |
note 2: test -- fun %*/% Subscript raised (line 509 of lib-- shows a bug. *)
|
neuper@52063
|
992 |
infix %*/%
|
neuper@52077
|
993 |
fun div_int_up p1 p1' p2 quot ns =
|
neuper@52063
|
994 |
let
|
neuper@52077
|
995 |
val (n, q) = fact_quot (lcoeff_up p1') (lcoeff_up p2);
|
neuper@52077
|
996 |
val remd = drop_tl_zeros ((p1' %* n) %-% (mult_to_deg (deg_up p1') p2 %* q));
|
neuper@52077
|
997 |
val zeros = for_quot_regarding p1 p2 p1' quot remd
|
neuper@52077
|
998 |
val _ = writeln_trace p1 p1' p2 quot n q zeros remd ns
|
neuper@52077
|
999 |
val quot = (replicate zeros 0) @ [q] @ (quot %* n)
|
neuper@52063
|
1000 |
(*div by hand uses q; * n multiplies the previous quot*)
|
neuper@52077
|
1001 |
val ns = n * ns
|
neuper@52063
|
1002 |
in
|
neuper@52077
|
1003 |
if deg_up remd >= deg_up p2
|
neuper@52077
|
1004 |
then div_int_up p1 remd p2 quot ns
|
neuper@52077
|
1005 |
else (ns, quot, remd)
|
neuper@52063
|
1006 |
end
|
neuper@52077
|
1007 |
|
neuper@52063
|
1008 |
fun p1 %*/% p2 =
|
neuper@52063
|
1009 |
let
|
neuper@52077
|
1010 |
val (n, q, r) = div_int_up p1 p1 p2 [] 1
|
neuper@52077
|
1011 |
val _ = if (! trace_div) then
|
neuper@52077
|
1012 |
(writeln " div .....................................................................";
|
neuper@52077
|
1013 |
writeln (" div " ^ PolyML.makestring n ^ " * " ^ PolyML.makestring p1 ^ " = " ^
|
neuper@52077
|
1014 |
PolyML.makestring q ^ " ** " ^ PolyML.makestring p2 ^ " ++ " ^ PolyML.makestring r);
|
neuper@52077
|
1015 |
writeln " div ....................................................................."
|
neuper@52077
|
1016 |
) else ()
|
neuper@52063
|
1017 |
in (n, q, r) end;
|
neuper@52063
|
1018 |
*}
|
neuper@52063
|
1019 |
text {* Here is the check of the algorithm with the values from above: *}
|
neuper@52063
|
1020 |
ML {*
|
neuper@52077
|
1021 |
trace_div := true;
|
neuper@52077
|
1022 |
trace_div_invariant := false;
|
neuper@52063
|
1023 |
val p1 = [2, 2, 2, 2, 2];
|
neuper@52063
|
1024 |
val p2 = [1, 2, 3];
|
neuper@52063
|
1025 |
val (n (* = 27*),
|
neuper@52063
|
1026 |
q (* = [8, 6, 18]*),
|
neuper@52063
|
1027 |
r (* = [46, 32]*)) = p1 %*/% p2;
|
neuper@52063
|
1028 |
if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed"
|
neuper@52063
|
1029 |
*}
|
neuper@52077
|
1030 |
text {* However, intermediate values increase exponentially: *}
|
neuper@52063
|
1031 |
ML {* (* example from [1].p.84 *)
|
neuper@52063
|
1032 |
val p1 = [~5,2,8,~3,~3,0,1,0,1];
|
neuper@52063
|
1033 |
val p2 = [21,~9,~4,5,0,3];
|
neuper@52063
|
1034 |
val (n (* = 27*),
|
neuper@52072
|
1035 |
q (* = [22, ~6, ~6, 9]*),
|
neuper@52072
|
1036 |
r (* = [~597, 378, 376, ~458, 6]*)) = p1 %*/% p2;
|
neuper@52077
|
1037 |
if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed 1"
|
neuper@52063
|
1038 |
*}
|
neuper@52077
|
1039 |
ML {* (* example used to develop algorithm by hand *)
|
neuper@52063
|
1040 |
val p1 = [2,2,2, 2,2,2];
|
neuper@52063
|
1041 |
val p2 = [7,8,6];
|
neuper@52063
|
1042 |
val (n (* = 162*),
|
neuper@52063
|
1043 |
q (* = [55, 15, ~18, 54]*),
|
neuper@52063
|
1044 |
r (* = [~61, ~221]*)) = p1 %*/% p2;
|
neuper@52077
|
1045 |
if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed 2"
|
neuper@52063
|
1046 |
*}
|
neuper@52063
|
1047 |
|
neuper@52063
|
1048 |
subsubsection {* The Euclidean Algorithm really works with polynomials? *}
|
neuper@52063
|
1049 |
text {*
|
neuper@52063
|
1050 |
The generalized division implemented above really can replace the
|
neuper@52063
|
1051 |
division of integers in the Eucliden Algorithm:
|
neuper@52063
|
1052 |
*}
|
neuper@52063
|
1053 |
ML {*
|
neuper@52077
|
1054 |
val trace_Euclid = Unsynchronized.ref true;
|
neuper@52077
|
1055 |
fun writeln_trace_Euclid (a: int list) (b: int list) (n: int) (q: int list) (r: int list) =
|
neuper@52077
|
1056 |
if (! trace_Euclid) then
|
neuper@52077
|
1057 |
writeln ("Euclid " ^ PolyML.makestring a ^ " * " ^ PolyML.makestring n ^
|
neuper@52077
|
1058 |
" == " ^ PolyML.makestring q ^
|
neuper@52077
|
1059 |
" ** " ^ PolyML.makestring b ^ " ++ " ^ PolyML.makestring r)
|
neuper@52077
|
1060 |
else ();
|
neuper@52077
|
1061 |
|
neuper@52077
|
1062 |
fun EUCLID_naive_up a [] = primitive_up a
|
neuper@52063
|
1063 |
| EUCLID_naive_up a b =
|
neuper@52063
|
1064 |
if deg_up a < deg_up b then EUCLID_naive_up b a else
|
neuper@52063
|
1065 |
let
|
neuper@52063
|
1066 |
val (n, q, r) = a %*/% b
|
neuper@52077
|
1067 |
val _ = writeln_trace_Euclid a b n q r
|
neuper@52063
|
1068 |
in EUCLID_naive_up b r : unipoly end;
|
neuper@52063
|
1069 |
*}
|
neuper@52063
|
1070 |
ML {*
|
neuper@52077
|
1071 |
(*The example used in the inital mail to Tobias Nipkow made univariate: y -> 1*)
|
neuper@52063
|
1072 |
trace_div := false;
|
neuper@52077
|
1073 |
trace_div_invariant := false;
|
neuper@52077
|
1074 |
trace_Euclid := true;
|
neuper@52063
|
1075 |
"--------------------------------";
|
neuper@52077
|
1076 |
val a = [0,~1,1]: unipoly; (* -x + x^2 = x *(-1 + x) *)
|
neuper@52077
|
1077 |
val b = [~1,0,1]: unipoly; (* -1 + x^2 = (1+x)*(-1 + x) *)
|
neuper@52063
|
1078 |
"--------------------------------";
|
neuper@52077
|
1079 |
EUCLID_naive_up a b; (* = [1, ~1]*) (*( 1 - x) *)
|
neuper@52063
|
1080 |
*}
|
neuper@52063
|
1081 |
ML {*
|
neuper@52063
|
1082 |
"--------------------------------";
|
neuper@52077
|
1083 |
val a = [~1,0,0,1]: unipoly; (* -1 + x^3 = (1+x+x^2)*(-1 + x) *)
|
neuper@52077
|
1084 |
val b = [0,~1,1]: unipoly; (* -x + x^2 = x *(-1 + x) *)
|
neuper@52063
|
1085 |
"--------------------------------";
|
neuper@52063
|
1086 |
EUCLID_naive_up a b; (* = [~1, 1]*) (*(-1 + x) *)
|
neuper@52063
|
1087 |
*}
|
neuper@52063
|
1088 |
ML {* (*NOW this works with generalized division*)
|
neuper@52063
|
1089 |
"--------------------------------";
|
neuper@52077
|
1090 |
val a = [0,~1,0,0,1]: unipoly; (* -x + x^4 = x*(1+x+x^2)*(-1 + x) *)
|
neuper@52077
|
1091 |
val b = [~2,2]: unipoly; (* -x + x^2 = 2 *(-1 + x) *)
|
neuper@52063
|
1092 |
"--------------------------------";
|
neuper@52063
|
1093 |
EUCLID_naive_up a b; (* = [~1, 1]*)
|
neuper@52063
|
1094 |
*}
|
neuper@52063
|
1095 |
ML {*
|
neuper@52077
|
1096 |
(* from [1].p.84; the generated coefficients are different for an unknown reason *)
|
neuper@52063
|
1097 |
"--------------------------------";
|
neuper@52077
|
1098 |
val a = [~5,2,8,~3,~3,0,1,0,1]: unipoly;
|
neuper@52077
|
1099 |
val b = [21,~9,~4,5,0,3]: unipoly;
|
neuper@52063
|
1100 |
"--------------------------------";
|
neuper@52063
|
1101 |
EUCLID_naive_up a b; (* = [1]*)
|
neuper@52063
|
1102 |
*}
|
neuper@52063
|
1103 |
ML {*
|
neuper@52072
|
1104 |
(* from [1].p.94 *)
|
neuper@52063
|
1105 |
"--------------------------------";
|
neuper@52077
|
1106 |
val a = [~18, ~15, ~20, 12, 20, ~13, 2]: unipoly;
|
neuper@52077
|
1107 |
val b = [8, 28, 22, ~11, ~14, 1, 2]: unipoly;
|
neuper@52063
|
1108 |
"--------------------------------";
|
neuper@52063
|
1109 |
EUCLID_naive_up a b; (* = [~2, ~1, 1]*)
|
neuper@52063
|
1110 |
*}
|
neuper@52063
|
1111 |
text {*
|
neuper@52077
|
1112 |
Note: The above algorithm still needs to take into account the case, where
|
neuper@52063
|
1113 |
constants can be factored out from polynomials, see [1].p.82-84.
|
neuper@52063
|
1114 |
*}
|
neuper@52063
|
1115 |
|
neuper@52078
|
1116 |
subsubsection {* Interactivity on Euclids Algorithm *}
|
neuper@52078
|
1117 |
text {*
|
neuper@52078
|
1118 |
Euclids Algorithm in modern algebraic notation explains itself if presented
|
neuper@52078
|
1119 |
this way:
|
neuper@52078
|
1120 |
*}
|
neuper@52078
|
1121 |
ML {*
|
neuper@52078
|
1122 |
trace_div := false;
|
neuper@52078
|
1123 |
trace_div_invariant := false;
|
neuper@52078
|
1124 |
trace_Euclid := true;
|
neuper@52078
|
1125 |
"--------------------------------";
|
neuper@52078
|
1126 |
val a = [~1,0,0,1]: unipoly; (* -1 + x^3 = (1+x+x^2)*(-1 + x) *)
|
neuper@52078
|
1127 |
val b = [0,~1,1]: unipoly; (* -x + x^2 = x *(-1 + x) *)
|
neuper@52078
|
1128 |
"--------------------------------";
|
neuper@52078
|
1129 |
EUCLID_naive_up a b; (* = [~1, 1]*) (*(-1 + x) *)
|
neuper@52078
|
1130 |
(*
|
neuper@52078
|
1131 |
Euclid [~1, 0, 0, 1] * 1 == [1, 1] ** [0, ~1, 1] ++ [~1, 1]
|
neuper@52078
|
1132 |
Euclid [0, ~1, 1] * 1 == [0, 1] ** [~1, 1] ++ []
|
neuper@52078
|
1133 |
*)
|
neuper@52078
|
1134 |
*}
|
neuper@52078
|
1135 |
ML {*
|
neuper@52078
|
1136 |
trace_Euclid := true;
|
neuper@52078
|
1137 |
(* from [1].p.94 *)
|
neuper@52078
|
1138 |
"--------------------------------";
|
neuper@52078
|
1139 |
val a = [~18, ~15, ~20, 12, 20, ~13, 2]: unipoly;
|
neuper@52078
|
1140 |
val b = [8, 28, 22, ~11, ~14, 1, 2]: unipoly;
|
neuper@52078
|
1141 |
"--------------------------------";
|
neuper@52078
|
1142 |
EUCLID_naive_up a b; (* = [~2, ~1, 1]*)
|
neuper@52078
|
1143 |
(*
|
neuper@52078
|
1144 |
Euclid [~18, ~15, ~20, 12, 20, ~13, 2] * 1 == [1] ** [8, 28, 22, ~11, ~14, 1, 2] ++ [~26, ~43, ~42, 23, 34, ~14]
|
neuper@52078
|
1145 |
Euclid [8, 28, 22, ~11, ~14, 1, 2] * 98 == [~41, ~14] ** [~26, ~43, ~42, 23, 34, ~14] ++ [~282, 617, ~168, ~723, 344]
|
neuper@52078
|
1146 |
Euclid [~26, ~43, ~42, 23, 34, ~14] * 59168 == [787, ~2408] ** [~282, 617, ~168, ~723, 344] ++ [~1316434, ~3708859, ~867104, 1525321]
|
neuper@52078
|
1147 |
Euclid [~282, 617, ~168, ~723, 344] * 6783102487 == [~2345549, 1529768] ** [~1316434, ~3708859, ~867104, 1525321] ++ [~5000595353600, ~2500297676800, 2500297676800]
|
neuper@52078
|
1148 |
Euclid [~1316434, ~3708859, ~867104, 1525321] * 7289497600 == [1919, 4447] ** [~5000595353600, ~2500297676800, 2500297676800] ++ []
|
neuper@52078
|
1149 |
*)
|
neuper@52078
|
1150 |
*}
|
neuper@52078
|
1151 |
text {*
|
neuper@52078
|
1152 |
A student only needs to observe, how the second factor on the right-hand side
|
neuper@52078
|
1153 |
shifts to the left-hand side; and he or she has to observe the remainder []
|
neuper@52078
|
1154 |
in the last line and take the preceeding remainder as the gcd.
|
neuper@52078
|
1155 |
|
neuper@52078
|
1156 |
In order to see that this recursion really computes the gcd, one needs to
|
neuper@52078
|
1157 |
know the fixpoint (or loop invariant):
|
neuper@52078
|
1158 |
*}
|
neuper@52079
|
1159 |
(*
|
neuper@52078
|
1160 |
lemma "a = (q *\<^isub>P b +\<^isub>P r) *\<^isub>s n \<and> degree b > degree r \<longrightarrow> gcd a b = gcd b r"
|
neuper@52078
|
1161 |
sorry
|
neuper@52079
|
1162 |
*)
|
neuper@52078
|
1163 |
text {*
|
neuper@52078
|
1164 |
The above lemma is analogous to
|
neuper@52078
|
1165 |
@{thm gcd_add_mult_int} = "gcd ?m (?k * ?m + ?n) = gcd ?m ?n" from HOL/GCD.thy
|
neuper@52078
|
1166 |
*}
|
neuper@52078
|
1167 |
|
neuper@52078
|
1168 |
subsubsection {* Interactivity and/or understanding of polynomial division *}
|
neuper@52078
|
1169 |
text {*
|
neuper@52078
|
1170 |
Division raises two requirements wrt. interactivity:
|
neuper@52078
|
1171 |
(1) a student wants to roughly understand what happens in the algorithm
|
neuper@52078
|
1172 |
(2) a student wants to reproduce the algorithm by hand.
|
neuper@52078
|
1173 |
|
neuper@52078
|
1174 |
Case (1) again calls for the fixpoint:
|
neuper@52078
|
1175 |
*}
|
neuper@52078
|
1176 |
ML {*
|
neuper@52078
|
1177 |
trace_div := false;
|
neuper@52078
|
1178 |
trace_div_invariant := true;
|
neuper@52078
|
1179 |
val p1 = [2, 2, 2, 2, 2];
|
neuper@52078
|
1180 |
val p2 = [1, 2, 3];
|
neuper@52078
|
1181 |
val (n (* = 27*),
|
neuper@52078
|
1182 |
q (* = [8, 6, 18]*),
|
neuper@52078
|
1183 |
r (* = [46, 32]*)) = p1 %*/% p2;
|
neuper@52078
|
1184 |
(*
|
neuper@52078
|
1185 |
div [2, 2, 2, 2, 2] * 3 = [0, 0, 2] ** [1, 2, 3] ++ [6, 6, 4, 2, 0]
|
neuper@52078
|
1186 |
div [2, 2, 2, 2, 2] * 9 = [0, 2, 6] ** [1, 2, 3] ++ [18, 16, 8, 0, 0]
|
neuper@52078
|
1187 |
div [2, 2, 2, 2, 2] * 27 = [8, 6, 18] ** [1, 2, 3] ++ [46, 32, 0, 0, 0]
|
neuper@52078
|
1188 |
*)
|
neuper@52078
|
1189 |
*}
|
neuper@52078
|
1190 |
text {*
|
neuper@52078
|
1191 |
The trace shows how the quotient, i.e. the first factor on the right-hand side,
|
neuper@52078
|
1192 |
is being increased step by step until the remainder has a degree lower than the
|
neuper@52078
|
1193 |
divisor.
|
neuper@52078
|
1194 |
|
neuper@52078
|
1195 |
This representation, however, is inappropriate for reproduction by hand.
|
neuper@52078
|
1196 |
Dividing polynomials by hand is taught in this format:
|
neuper@52078
|
1197 |
*}
|
neuper@52078
|
1198 |
ML {*
|
neuper@52078
|
1199 |
trace_div := true;
|
neuper@52078
|
1200 |
trace_div_invariant := false;
|
neuper@52078
|
1201 |
val p1 = [2, 2, 2, 2, 2];
|
neuper@52078
|
1202 |
val p2 = [1, 2, 3];
|
neuper@52078
|
1203 |
val (n (* = 27*),
|
neuper@52078
|
1204 |
q (* = [8, 6, 18]*),
|
neuper@52078
|
1205 |
r (* = [46, 32]*)) = p1 %*/% p2;
|
neuper@52078
|
1206 |
(*
|
neuper@52078
|
1207 |
div [2, 2, 2, 2, 2] * 3
|
neuper@52078
|
1208 |
div ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
neuper@52078
|
1209 |
div [6, 6, 6, 6, 6]
|
neuper@52078
|
1210 |
div-- [0, 0, 2, 4, 6] <--- [1, 2, 3] * 2
|
neuper@52078
|
1211 |
div ---------------------------
|
neuper@52078
|
1212 |
div [6, 6, 4, 2, 0] quot = [] @ [2] @ (3 * [])
|
neuper@52078
|
1213 |
div [6, 6, 4, 2] * 3
|
neuper@52078
|
1214 |
div ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
neuper@52078
|
1215 |
div [18, 18, 12, 6]
|
neuper@52078
|
1216 |
div-- [0, 2, 4, 6] <--- [1, 2, 3] * 2
|
neuper@52078
|
1217 |
div ---------------------------
|
neuper@52078
|
1218 |
div [18, 16, 8, 0] quot = [] @ [2] @ (3 * [2])
|
neuper@52078
|
1219 |
div [18, 16, 8] * 3
|
neuper@52078
|
1220 |
div ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
neuper@52078
|
1221 |
div [54, 48, 24]
|
neuper@52078
|
1222 |
div-- [8, 16, 24] <--- [1, 2, 3] * 8
|
neuper@52078
|
1223 |
div ---------------------------
|
neuper@52078
|
1224 |
div [46, 32, 0] quot = [] @ [8] @ (3 * [2, 6])
|
neuper@52078
|
1225 |
div .....................................................................
|
neuper@52078
|
1226 |
div 27 * [2, 2, 2, 2, 2] = [8, 6, 18] ** [1, 2, 3] ++ [46, 32]
|
neuper@52078
|
1227 |
div .....................................................................
|
neuper@52078
|
1228 |
*)
|
neuper@52078
|
1229 |
*}
|
neuper@52078
|
1230 |
text {*
|
neuper@52078
|
1231 |
The above algorithm is unusual in that the dividend first is multiplied
|
neuper@52078
|
1232 |
such that the quotient can have integer coefficients ("generalised division").
|
neuper@52078
|
1233 |
Thus the quotient needs to be multiplied subsequently as well.
|
neuper@52078
|
1234 |
The algorithm also has to handle cases, where zeros occur in the quotient:
|
neuper@52078
|
1235 |
*}
|
neuper@52078
|
1236 |
ML {*
|
neuper@52078
|
1237 |
trace_div := true;
|
neuper@52078
|
1238 |
val p1 = [3,2,2,2,1,1,1,1];
|
neuper@52078
|
1239 |
val p2 = [1,1,1,1];
|
neuper@52078
|
1240 |
val (n (* = 1*),
|
neuper@52078
|
1241 |
q (* = [2, 0, 0, 0, 1]*),
|
neuper@52078
|
1242 |
r (* = [1]*)) = p1 %*/% p2;
|
neuper@52078
|
1243 |
(*
|
neuper@52078
|
1244 |
div [3, 2, 2, 2, 1, 1, 1, 1] * 1
|
neuper@52078
|
1245 |
div ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
neuper@52078
|
1246 |
div [3, 2, 2, 2, 1, 1, 1, 1]
|
neuper@52078
|
1247 |
div-- [0, 0, 0, 0, 1, 1, 1, 1] <--- [1, 1, 1, 1] * 1
|
neuper@52078
|
1248 |
div ---------------------------
|
neuper@52078
|
1249 |
div [3, 2, 2, 2, 0, 0, 0, 0] quot = [0, 0, 0] @ [1] @ (1 * [])
|
neuper@52078
|
1250 |
div [3, 2, 2, 2] * 1
|
neuper@52078
|
1251 |
div ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
neuper@52078
|
1252 |
div [3, 2, 2, 2]
|
neuper@52078
|
1253 |
div-- [2, 2, 2, 2] <--- [1, 1, 1, 1] * 2
|
neuper@52078
|
1254 |
div ---------------------------
|
neuper@52078
|
1255 |
div [1, 0, 0, 0] quot = [] @ [2] @ (1 * [0, 0, 0, 1])
|
neuper@52078
|
1256 |
div .....................................................................
|
neuper@52078
|
1257 |
div 1 * [3, 2, 2, 2, 1, 1, 1, 1] = [2, 0, 0, 0, 1] ** [1, 1, 1, 1] ++ [1]
|
neuper@52078
|
1258 |
div .....................................................................
|
neuper@52078
|
1259 |
*)
|
neuper@52078
|
1260 |
*}
|
neuper@52078
|
1261 |
|
neuper@52063
|
1262 |
subsubsection {* *}
|
neuper@52063
|
1263 |
text {*
|
neuper@52063
|
1264 |
|
neuper@52063
|
1265 |
*}
|
neuper@52063
|
1266 |
ML {*
|
neuper@52063
|
1267 |
|
neuper@52063
|
1268 |
*}
|
neuper@52063
|
1269 |
|
neuper@52063
|
1270 |
subsection {* Technical requirements raised by pedagogical questions *}
|
neuper@52078
|
1271 |
subsubsection {* Calculations might be different from interpretation *}
|
neuper@52078
|
1272 |
text {*
|
neuper@52078
|
1273 |
Until now the demonstration cases for Lucas-Interpretation were
|
neuper@52078
|
1274 |
simplification, equation solving and ad-hoc algorithms in Signal Processing
|
neuper@52078
|
1275 |
and Structural Engineering. In these cases execution of the functional program
|
neuper@52078
|
1276 |
coincided with the calculation generated as a side-effect of interpretation:
|
neuper@52078
|
1277 |
The tactics "Take", "Rewrite", "Rewrite_Set" and "Substitute", when executed,
|
neuper@52078
|
1278 |
immediately created respective output and handed over control to the user.
|
neuper@52078
|
1279 |
Appropriate user-input finally created a "calculation" very close to
|
neuper@52078
|
1280 |
what would be written on paper.
|
neuper@52078
|
1281 |
|
neuper@52078
|
1282 |
Now, in Computer Algebra, the programs contain lots of sophisticated functions.
|
neuper@52078
|
1283 |
Simple traces, for instance arguments and values of function calls, comprise too
|
neuper@52078
|
1284 |
many data to be intuitively comprehensible.
|
neuper@52078
|
1285 |
|
neuper@52078
|
1286 |
The examples in the case study suggest to have mechanisms for creating
|
neuper@52078
|
1287 |
steps for calculations, which are somehow abstracted from pure execution.
|
neuper@52078
|
1288 |
Interestingly these abstractions are close to the fixpoints/invariants
|
neuper@52078
|
1289 |
of the respective recursive functions --- see
|
neuper@52078
|
1290 |
section "Interactivity on Euclids Algorithm" and section "Interactivity and/or
|
neuper@52078
|
1291 |
understanding of polynomial division" above.
|
neuper@52078
|
1292 |
|
neuper@52078
|
1293 |
Such mechanisms might be specific tactics, which take data from the
|
neuper@52078
|
1294 |
local environment and the local context for presentation to the student.
|
neuper@52078
|
1295 |
Thus one such tactic is related to several lines of code in a function.
|
neuper@52078
|
1296 |
These tactics would not contribute to automated execution as done by "value",
|
neuper@52078
|
1297 |
for instance.
|
neuper@52078
|
1298 |
However, the effect of these tactics for resuming execution upon completion
|
neuper@52078
|
1299 |
of input needs to be clarified.
|
neuper@52078
|
1300 |
*}
|
neuper@52078
|
1301 |
|
neuper@52078
|
1302 |
subsubsection {* One function might generate various representations *}
|
neuper@52078
|
1303 |
text {*
|
neuper@52078
|
1304 |
The section "Interactivity and/or understanding of polynomial division" suggests
|
neuper@52078
|
1305 |
to have at least two kinds of representation for "calculations":
|
neuper@52078
|
1306 |
(1) a student wants to roughly understand what happens in the algorithm
|
neuper@52078
|
1307 |
(2) a student wants to reproduce the algorithm by hand.
|
neuper@52078
|
1308 |
|
neuper@52078
|
1309 |
A solution to this requirement of variety is to create related functions
|
neuper@52078
|
1310 |
with different combinations of the specific tactics as mentioned above.
|
neuper@52078
|
1311 |
These functions can be selected by the dialogue module. The computational
|
neuper@52078
|
1312 |
effect of such related functions is the same, the interactive behaviour
|
neuper@52078
|
1313 |
is different (much more different as variations of "Rewrite", for instance.)
|
neuper@52078
|
1314 |
*}
|
neuper@52078
|
1315 |
|
neuper@52078
|
1316 |
subsubsection {* One "specific tactic" might create several lines of calculation *}
|
neuper@52078
|
1317 |
text {*
|
neuper@52078
|
1318 |
The case of "division of polynomial" suggests to have one "specific tactic"
|
neuper@52078
|
1319 |
as mentioned above, which has one location in the function (e.g. "writeln_trace"
|
neuper@52078
|
1320 |
in "div_int_up"), but creates several lines of a calculation; see the example in
|
neuper@52078
|
1321 |
section "Can division of polynomials be adapted":
|
neuper@52078
|
1322 |
|
neuper@52078
|
1323 |
[ 2, 2, 2, 2, 2] : [1, 2, 3] =
|
neuper@52078
|
1324 |
3*[ 2, 2, 2, 2, 2] *3
|
neuper@52078
|
1325 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:
|
neuper@52078
|
1326 |
[ 6, 6, 6, 6, 6] :
|
neuper@52078
|
1327 |
--[ 0, 0, 2, 4, 6] <--........ * [ 2]
|
neuper@52078
|
1328 |
------------------------ :
|
neuper@52078
|
1329 |
[ 6, 6, 4, 2, 0] = [ 2]*3
|
neuper@52078
|
1330 |
:
|
neuper@52078
|
1331 |
:
|
neuper@52078
|
1332 |
|
neuper@52078
|
1333 |
Interaction on these lines should resemble calculation with paper and pencil.
|
neuper@52078
|
1334 |
*}
|
neuper@52063
|
1335 |
|
neuper@52063
|
1336 |
subsubsection {* *}
|
neuper@52063
|
1337 |
text {*
|
neuper@52063
|
1338 |
|
neuper@52063
|
1339 |
*}
|
neuper@52063
|
1340 |
ML {*
|
neuper@52063
|
1341 |
|
neuper@52063
|
1342 |
*}
|
neuper@52063
|
1343 |
|
meindl_diana@48797
|
1344 |
end
|
meindl_diana@48797
|
1345 |
|
meindl_diana@48797
|
1346 |
|