1 (* rationals, i.e. fractions of multivariate polynomials over the real field
4 Use is subject to license terms.
5 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
6 10 20 30 40 50 60 70 80 90 100
9 theory GCD_Poly imports Complex_Main begin
13 F. Winkler, Polynomial algorithyms in computer algebra. Springer 1996.
14 by page 93 and 95. This is the final version of the Isabelle theory,
15 which conforms with Isabelle's coding standards and with requirements
16 in terms of logics appropriate for a theorem prover.
19 section {* Predicates for the specifications *}
21 This thesis provides program code to be integrated into the Isabelle
22 distribution. Thus the code will be verified against specifications, and
23 consequently the code, as \emph{explicit} specification of how to calculate
24 results, is accompanied by \emph{implicit} specifications describing the
25 properties of the program code.
27 For describing the properties of code the format of Haskell will be used.
28 Haskell \cite{TODO} is front-of-the-wave in research on functional
29 programming and as such advances beyond SML. The advancement relevant for
30 this thesis is the concise combination of explicit and implicit specification
31 of a function. The final result of the code developed within this thesis is
32 described as follows in Haskell format:
34 \begin{verbatim} % TODO: find respective Isabelle/LaTeX features
35 gcd_poly :: int poly \<Rightarrow> int poly \<Rightarrow> int poly
37 assumes and p1 \<noteq> [:0:] and p2 \<noteq> [:0:]
38 yields d dvd p1 \<and> d dvd p2 \<and> \<forall>d'. (d' dvd p1 \<and> d' dvd p2) \<longrightarrow> d' \<le> d
41 The above specification combines Haskell format with the term language of
42 Isabelle\HOL. The latter provides a large and ever growing library of
43 mechanised mathematics knowledge, which this thesis re-uses in order to
44 prepare for integration of the code.
46 In the sequel there is a brief account of the notions required for the
47 specification of the functions implemented within this thesis. *}
49 subsection {* The Isabelle types required *}
51 bool, nat, int, option, THE, list, 'a poly
54 subsection {* The connectives required *}
56 logic ~~/doc/tutorial.pdf p.203
57 arithmetic < \<le> + - * ^ div \<bar>_\<bar>
60 subsection {* The Isabelle functions required *}
65 Primes.next_prime_bound
69 List.length, List.member
71 List.nth, List.take, List.drop
78 section {* gcd for univariate polynomials *}
80 type unipoly = int list;
81 val a = [~18, ~15, ~20, 12, 20, ~13, 2];
82 val b = [8, 28, 22, ~11, ~14, 1, 2];
85 subsection {* auxiliary functions *}
88 5 div 2 = 2; ~5 div 2 = ~3; BUT WEE NEED ~5 div2 2 = ~2; *)
91 if a div b < 0 then (abs a) div (abs b) * ~1 else a div b
93 (* why has "last" been removed since 2002, although "last" is in List.thy ? *)
94 fun last ls = (hd o rev) ls
96 (* drop tail elements equal 0 *)
97 fun drop_hd_zeros (0 :: ps) = drop_hd_zeros ps
99 fun drop_tl_zeros p = (rev o drop_hd_zeros o rev) p;
102 subsection {* modulo calculations for integers *}
104 (* mod_inv :: int \<Rightarrow> nat
107 yields r * s mod m = 1 *)
110 fun modi (r, rold, m, rinv) =
111 if m <= rinv then 0 else
114 else modi (rold * (rinv + 1), rold, m, rinv + 1)
115 in modi (r, r, m, 1) end
117 (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO
120 yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*)
121 fun mod_div a b m = a * (mod_inv b m) mod m
123 (* required just for one approximation:
124 approx_root :: nat \<Rightarrow> nat
127 yields r * r \<ge> a *)
130 fun root1 a n = if n * n < a then root1 a (n + 1) else n
133 (* chinese_remainder :: int * nat \<Rightarrow> int * nat \<Rightarrow> nat
134 chinese_remainder (r1, m1) (r2, m2) = r
136 yields r = r1 mod m1 \<and> r = r2 mod m2 *)
137 fun chinese_remainder (r1, m1) (r2, m2) =
139 (((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) * m1;
142 subsection {* creation of lists of primes for efficiency *}
144 (* is_prime :: int list \<Rightarrow> int \<Rightarrow> bool
146 assumes max ps < n \<and> n \<le> (max ps)^2 \<and>
147 (* FIXME: really ^^^^^^^^^^^^^^^? *)
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)
149 yields Primes.prime n *)
150 fun is_prime [] _ = true
152 if n mod (List.hd prs) <> 0 then is_prime (List.tl prs) n else false
154 (* make_primes is just local to primes_upto only:
155 make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list
156 make_primes ps last_p n = pps
157 assumes last_p = maxs ps \<and> (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
158 (\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)
159 yields n \<le> maxs pps \<and> (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p) \<and>
160 (\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
161 fun make_primes ps last_p n =
162 if n <= last ps then ps else
163 if is_prime ps (last_p + 2)
164 then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
165 else make_primes ps (last_p + 2) n
167 (* primes_upto :: nat \<Rightarrow> nat list
170 yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
171 n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
172 (\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
173 fun primes_upto n = if n < 3 then [2] else make_primes [2, 3] 3 n
175 (* find the next prime greater p not dividing the number n:
176 next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
177 n1 next_prime_not_dvd n2 = p
179 yields p is_prime \<and> n1 < p \<and> \<not> p dvd n2 \<and> (* smallest with these properties... *)
180 (\<forall> p'. (p' is_prime \<and> n1 < p' \<and> \<not> p' dvd n2) \<longrightarrow> p \<le> p') *)
181 infix next_prime_not_dvd;
182 fun n1 next_prime_not_dvd n2 =
184 val ps = primes_upto (n1 + 1)
185 val nxt = nth ps (List.length ps - 1);
187 if n2 mod nxt <> 0 then nxt else nxt next_prime_not_dvd n2
191 subsection {* basic notions for univariate polynomials *}
193 (* leading coefficient, includes drop_lc0_up ?FIXME130507 *)
194 fun lcoeff_up (p: unipoly) = (last o drop_tl_zeros) p
196 (* drop leading coefficients equal 0 *)
197 fun drop_lc0_up (p: unipoly) = drop_tl_zeros p : unipoly;
199 (* degree, includes drop_lc0_up ?FIXME130507 *)
200 fun deg_up (p: unipoly) = ((Integer.add ~1) o length o drop_lc0_up) p;
202 (* normalise a unipoly such that lcoeff_up mod m = 1.
203 normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
204 normalise [p_0, .., p_k] m = [q_0, .., q_k]
206 yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i \<and> (p_k * t) mod m = 1 *)
207 fun normalise (p: unipoly) (m: int) =
209 val p = drop_lc0_up p
210 val lcp = lcoeff_up p;
213 then [mod_div (nth p i) lcp m] @ nrm
214 else norm ([mod_div (nth p i) lcp m] @ nrm) (i - 1) ;
216 if List.length p = 0 then p else norm [] (List.length p - 1)
220 subsection {* addition, multiplication, division *}
222 (* scalar multiplication *)
224 fun (p: unipoly) %* n = map (Integer.mult n) p : unipoly
228 fun (p: unipoly) %/ n =
229 let fun swapf f a b = f b a (* swaps the arguments of a function *)
230 in map (swapf (curry op div2) n) p : unipoly end;
234 fun (p1: unipoly) %-% (p2: unipoly) =
235 map2 (curry op -) p1 (p2 @ replicate (List.length p1 - List.length p2) 0) : unipoly
237 (* dvd for univariate polynomials *)
240 fun [d] %|% [p] = (p mod d = 0)
241 | (ds: unipoly) %|% (ps: unipoly) =
243 val ds = drop_lc0_up ds;
244 val ps = drop_lc0_up ps;
245 val d000 = (replicate (List.length ps - List.length ds) 0) @ ds ;
246 val quot = (lcoeff_up ps) div2 (lcoeff_up d000);
247 val rest = drop_lc0_up (ps %-% (d000 %* quot));
249 if rest = [] then true else
250 if rest <> [] andalso quot <> 0 andalso List.length ds <= List.length rest
256 subsection {* normalisation and Landau-Mignotte bound *}
258 (* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
259 normalise [p_0, .., p_k] m = [q_0, .., q_k]
261 yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
262 (where |^ x ^| means round x up to the next greater integer) *)
263 fun centr_up (p: unipoly) (m: int) =
266 val mid = if m mod mi = 0 then mi else mi + 1;
267 fun centr m mid p_i = if mid < p_i then p_i - m else p_i;
268 in map (centr m mid) p : unipoly end;
270 (*** landau mignotte bound (Winkler p91)
271 every coefficient of the gcd of a and b in Z[x] is bounded in absolute value by
272 2^{min(m,n)} * gcd(a_m,b_n) * min(1/|a_m|
273 * root(sum_{i=0}^{m}a_i^2) ,1/|b_n| * root( sum_{i=0}^{n}b_i^2 ) ***)
275 (* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
276 sum_lmb [p_0, .., p_k] e = s
278 yields. p_0^e + p_1^e + ... + p_k^e *)
279 fun sum_lmb (p: unipoly) exp = fold ((curry op +) o (Integer.pow exp)) p 0;
281 (* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly \<Rightarrow> unipoly \<Rightarrow> int
282 LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
284 yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
285 min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)
286 fun LANDAU_MIGNOTTE_bound (p1: unipoly) (p2: unipoly) =
287 (Integer.pow (Integer.min (deg_up p1) (deg_up p2)) 2) * (abs (Integer.gcd (lcoeff_up p1) (lcoeff_up p2))) *
288 (Int.min (abs((approx_root (sum_lmb p1 2)) div2 ~(lcoeff_up p1)),
289 abs((approx_root (sum_lmb p2 2)) div2 ~(lcoeff_up p2))));
292 subsection {* modulo calculations on polynomials *}
294 (* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
295 chinese_remainder_up (m1, m2) (p1, p2) = p
296 assume m1, m2 relatively prime
297 yields p1 = p mod m1 \<and> p2 = p mod m2 *)
298 fun chinese_remainder_up (m1, m2) (p1: unipoly, p2: unipoly) =
300 fun chinese_rem (m1, m2) (p_1i, p_2i) = chinese_remainder (p_1i, m1) (p_2i, m2)
301 in map (chinese_rem (m1, m2)) (p1 ~~ p2): unipoly end
303 (* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
304 mod_up [p1, p2, ..., pk] m = up
306 yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
308 fun (p : unipoly) mod_up m =
310 fun mod' m p = (curry op mod) p m
311 in map (mod' m) p : unipoly end
313 (* euclidean algoritm in Z_p[x].
314 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
315 mod_up_gcd p1 p2 m = g
317 yields gcd ((p1 mod m), (p2 mod m)) = g *)
318 fun mod_up_gcd (p1: unipoly) (p2: unipoly) (m: int) =
320 val p1m = p1 mod_up m;
321 val p2m = drop_lc0_up (p2 mod_up m);
322 val p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
323 val quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
324 val rest = drop_lc0_up ((p1m %-% (p2n %* quot)) mod_up m)
326 if rest = [] then p2 else
327 if length rest < List.length p2
328 then mod_up_gcd p2 rest m
329 else mod_up_gcd rest p2 m
332 (* primitive_up :: unipoly \<Rightarrow> unipoly
335 yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
336 fun primitive_up ([c] : unipoly) =
337 if c = 0 then ([0] : unipoly) else ([1] : unipoly)
340 val d = abs (Integer.gcds p)
342 if d = 1 then p else p %/ d
346 subsection {* gcd_up, code from [1] p.93 *}
348 (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
349 try_new_prime_up p1 p2 d M P g p = new_g
350 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
351 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
352 p1 is primitive \<and> p2 is primitive
353 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
355 argumentns "a b d M P g p" named according to [1] p.93 *)
356 fun try_new_prime_up a b d M P g p =
359 val p = p next_prime_not_dvd d
360 val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
365 if deg_up g_p < deg_up g
367 if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p
369 if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else
372 val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
373 in try_new_prime_up a b d M P g p end
376 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
377 HENSEL_lifting_up p1 p2 d M p = gcd
378 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
379 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and>
380 p1 is primitive \<and> p2 is primitive
381 yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
383 argumentns "a b d M p" named according to [1] p.93 *)
384 fun HENSEL_lifting_up a b d M p =
386 val _ = writeln ("HENSEL-univ d = " ^ PolyML.makestring d ^
387 ", M = " ^ PolyML.makestring (2 * d * LANDAU_MIGNOTTE_bound a b) ^
388 ", p = " ^ PolyML.makestring p)
389 val p = p next_prime_not_dvd d
390 val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
392 if deg_up g_p = 0 then [1] else
394 val g = primitive_up (try_new_prime_up a b d M p g_p p)
396 if (g %|% a) andalso (g %|% b) then
397 (writeln ("HENSEL-univ -------------------> " ^ PolyML.makestring g);
400 else HENSEL_lifting_up a b d M p
404 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
406 assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
407 a is primitive \<and> b is primitive
408 yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
410 let val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))
413 (writeln ("gcd_up calls HENSEL-univ");
414 HENSEL_lifting_up a b d (2 * d * LANDAU_MIGNOTTE_bound a b) 1
419 section {* gcd for multivariate polynomials *}
421 type monom = (int * int list);
422 type poly = monom list;
425 subsection {* auxiliary functions *}
427 fun land a b = a andalso b
428 fun all_geq a b = fold land (map2 (curry op >=) a b) true
429 fun all_geq' i es = fold land (map ((curry op <=) i) es) true
430 fun maxs is = fold Integer.max is (hd is);
431 fun div2_swapped d i = i div2 d
432 fun nth_swapped i ls = nth ls i;
433 fun drop_last ls = (rev o tl o rev) ls
434 fun nth_ins n i xs = (*analogous to nth_drop: why not in library.ML?*)
435 List.take (xs, n) @ [i] @ List.drop (xs, n);
438 subsection {* basic notions for multivariate polynomials *}
440 fun lcoeff (p: poly) = (fst o hd o rev) p
441 fun lexp (p: poly) = (snd o hd o rev) p
442 fun lmonom (p: poly) = (hd o rev) p
444 fun cero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;
446 (* drop leading coefficients equal 0 TODO: signature?*)
447 fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then cero_poly (length es) else p
448 | drop_lc0 p = if lcoeff p = 0 then drop_lc0 (drop_last p) else p
450 (* gcd of coefficients WN find right location in file *)
451 fun gcd_coeff (up: unipoly) = abs (Integer.gcds up)
453 (* given p = 3x+4z = [(3,[1,0]),(4,[0,1])] with x<z
454 add_variable p 2 => [(3,[1,0,0]),(4,[0,0,1])] with x<y<z *)
455 fun add_variable (p: poly) var = map (apsnd (nth_ins var 0)) p : poly
458 subsection {* ordering and other auxiliary funcions *}
460 (* monomial order: lexicographic order on exponents *)
461 fun lex_ord ((_, a): monom) ((_, b): monom) =
462 let val d = drop_tl_zeros (a %-% b)
464 if d = [] then false else last d > 0
466 fun lex_ord' ((_, a): monom, (_, b): monom) =
467 let val d = drop_tl_zeros (a %-% b)
470 else if last d > 0 then GREATER else LESS
473 (* add monomials in ordered polynomal *)
474 fun add_monoms (p: poly) =
476 fun add [(0, es)] [] = cero_poly (length es)
477 | add [(0, _)] (new: poly) = new
478 | add [m] (new: poly) = new @ [m]
479 | add ((c1, es1) :: (c2, es2) :: ms) (new: poly) =
481 then add ([(c1 + c2, es1)] @ ms) new
484 then add ((c2, es2) :: ms) new
485 else add ((c2, es2) :: ms) (new @ [(c1, es1)])
486 in add p [] : poly end
488 (* brings the polynom in correct order and adds monomials *)
489 fun order (p: poly) = (add_monoms o (sort lex_ord')) p;
491 (* largest exponent of variable at location var *)
492 fun max_deg_var (p: poly) var = (maxs o (map ((nth_swapped var) o snd))) p;
495 subsection {* evaluation, translation uni--multivariate, etc *}
497 (* evaluate a polynomial in a certain variable and remove this var from the exp-list *)
498 fun eval (p: poly) var value = (* TODO use "map" instead*)
500 fun eval [] _ _ new = order new
501 | eval ((c, es)::ps) var value new =
502 eval ps var value (new @ [((Integer.pow (nth es var) value) * c, nth_drop var es)]);
503 in eval p var value [] end
505 (* transform a multivariate to a univariate polynomial *)
506 fun multi_to_uni (p as (_, es)::_ : poly) =
509 let fun transfer ([]: poly) (up: unipoly) = up
510 | transfer (mp as (c, es)::ps) up =
512 then transfer ps (up @ [c])
513 else transfer mp (up @ [0])
514 in transfer (order p) [] end
515 else error "Polynom has more than one variable!";
517 (* transform a univariate to a multivariate polynomial *)
518 fun uni_to_multi (up: unipoly) =
519 let fun transfer ([]: unipoly) (mp: poly) (_: int) = add_monoms mp
520 | transfer up mp var =
521 transfer (tl up) (mp @ [(hd up, [var])]) (var + 1)
522 in transfer up [] 0 end
524 (* multiply a poly with a variable represented by a unipoly (at position var)
525 e.g new z: (3x + 2y)*(z - 4) TODO fold / map ? *)
526 fun mult_with_new_var ([]: poly) (_: unipoly) (_: int) = []
527 | mult_with_new_var mp up var =
529 fun mult ([]: poly) (_: unipoly) (_: int) (new: poly) (_: int) = add_monoms new
530 | mult mp up var new _ =
533 fun mult' (_: poly) ([]: unipoly) (_) (new': poly) (_: int) = order new'
534 | mult' (mp' as (c, es)::ps) up' var new' c2' =
535 let val (first, last) = chop var es
536 in mult' mp' (tl up') var
537 (new' @ [(c * hd up', first @ [c2'] @ last)]) (c2' + 1)
539 in mult (tl mp) up var (new @ (mult' mp up var [] 0)) (0) end
540 in mult mp up var [] 0 : poly end
542 fun greater_deg (es1: int list) (es2: int list) =
543 if es1 = [] andalso es2 =[] then 0
544 else if (nth es1 (length es1 -1)) = (nth es2 (length es1 -1) )
545 then greater_deg (nth_drop (length es1 -1) es1) (nth_drop (length es1 -1) es2)
546 else if (nth es1 (length es1 -1)) > (nth es2 (length es1 -1))
550 fun find_new_r (p1 as (_, es1)::_ : poly) (p2 as (_, es2)::_ : poly) old =
551 let val p1' as (_, es1')::_ = eval p1 (length es1 - 2) (old + 1)
552 val p2' as (_, es2')::_ = eval p2 (length es2 - 2) (old + 1);
554 if max_deg_var p1' (length es1' - 1) = max_deg_var p1 (length es1 - 1) orelse
555 max_deg_var p2' (length es2' - 1) = max_deg_var p2 (length es1 - 1)
557 else find_new_r p1 p2 (old + 1)
561 subsection {* addition, multiplication, division *}
563 (* resulting 0 coefficients are removed by add_monoms *)
565 fun (p: poly) %%/ d = (add_monoms o (map (apfst (div2_swapped d)))) p;
568 fun (p: poly) %%* i = (add_monoms o (map (apfst (Integer.mult i)))) p;
571 (* the same algorithms as done by hand
572 fun ([]: poly) %%+%% (mvp2: poly) = mvp2
576 fun plus [] p2 new = order (new @ p2)
577 | plus p1 [] new = order (new @ p1)
578 | plus (p1 as (c1, es1)::ps1) (p2 as (c2, es2)::ps2) new =
579 if greater_deg es1 es2 = 0
580 then plus ps1 ps2 (new @ [((c1 + c2), es1)])
582 if greater_deg es1 es2 = 1
583 then plus p1 ps2 (new @ [(c2, es2)])
584 else plus ps1 p2 (new @ [(c1, es1)])
587 fun (p1 : poly) %%+%% (p2 : poly) = order (p1 @ p2)
590 fun p1 %%-%% (p2: poly) = p1 %%+%% (map (apfst (Integer.mult ~1)) p2)
593 (* the same algorithms as done by hand
594 fun (p1: poly) %%*%% (p2: poly) =
596 fun mult ([]: poly) (_: poly) (_: poly) (new: poly) = order new
597 | mult p1 [] regular_p2 new = mult (tl p1) regular_p2 regular_p2 new
598 | mult (p1 as (c1, es1)::_) ((c2, es2)::ps2) regular_p2 (new: poly) =
599 mult p1 ps2 regular_p2 (new @ [(c1 * c2, map2 Integer.add es1 es2)])
601 if length p1 > length p2
602 then mult p1 p2 p2 []
603 else mult p2 p1 p1 []
606 fun mult_monom ((c1, es1): monom) ((c2, es2): monom) =
607 (c1 * c2, map2 Integer.add es1 es2): monom;
608 fun mult_monoms (ms: poly) m = map (mult_monom m) ms : poly
609 fun p1 %%*%% (p2 : poly) = order (fold (curry op @) (map (mult_monoms p1) p2) [])
612 (* the same algorithms as done by hand; TODO fold / map ? *)
613 fun ([(c1, es1)]: poly) %%|%% ([(c2, es2)]: poly) = all_geq es2 es1 andalso c2 mod c1 = 0
614 | _ %%|%% [(0,_)] = true
616 if [lmonom p1] %%|%% [lmonom p2]
620 ([((lcoeff p2) div2 (lcoeff p1), (lexp p2) %-% (lexp p1))]
624 (* %%/%% :: poly \<Rightarrow> poly \<Rightarrow> (poly \<times> poly)
626 assumes b \<noteq> [] \<and>
627 yields THE c r. c *\<^isub>p b +\<^isub>p r = a *)
629 fun (p1: poly) %%/%% (p2: poly) =
631 fun div_up p1 p2 quot_old =
633 val p1 as (c, es)::_ = drop_lc0 (order p1);
634 val p2 = drop_lc0 (order p2);
635 val [c] = cero_poly (length es);
636 val d = (replicate (List.length p1 - List.length p2) c) @ p2 ;
637 val quot = [(lcoeff p1 div2 lcoeff d, lexp p1 %-% lexp p2)];
638 val rest = drop_lc0 (p1 %%-%% (d %%*%% quot));
640 if rest <> [c] andalso all_geq' 0 (lexp rest %-% lexp p2)
641 then div_up rest p2 (quot @ quot_old)
642 else (quot @ quot_old : poly, rest)
644 in div_up p1 p2 [] end
647 subsection {* Newton interpolation *}
648 ML{* (* first step *)
649 fun newton_first (x: int list) (f: poly list) (ord: int) =
651 val new_vals = [(hd x) * ~1, 1];
652 val p = (add_variable (hd f) ord) %%+%%
654 (((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x))) new_vals ord);
655 val steps = [((nth f 1) %%-%% (hd f)) %%/ ((nth x 1) - (hd x))]
656 in (p, new_vals, steps) end
658 fun next_step ([]: poly list) (new_steps: poly list) (_: int list) = new_steps
659 | next_step steps new_steps x =
662 [((last new_steps) %%-%% (hd steps)) %%/ ((last x) - (nth x (length x - 3)))])
663 (nth_drop (length x -2) x)
665 fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) ord =
668 let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] ord
669 in (p', new_vals, steps) end
673 multi_to_uni ((uni_to_multi up) %%*%% (uni_to_multi [(nth x (length x -2) ) * ~1, 1]));
674 val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/
675 ((last x) - (nth x (length x - 2)))];
676 val steps = next_step steps new_steps x;
677 val p' = p %%+%% (mult_with_new_var (last steps) new_vals ord);
678 in (order p', new_vals, steps) end
681 subsection {* gcd_poly algorithm, code for [1] p.95 *}
683 (* naming of n, M, m, r, ... follows Winkler p. 95
684 n: amount of variables
686 fun gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r =
687 if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else
690 let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
691 in try_new_r a b n M r [] [] end
693 let val (a', b') = (multi_to_uni a, multi_to_uni b)
695 (* a b are made primitive and this is compensated by multiplying with gcd of divisors *)
696 uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %*
697 (abs (Integer.gcd (gcd_coeff a') (gcd_coeff b'))))
700 (* fn: poly -> poly -> int ->
701 int -> int -> int list -> poly list -> poly*)
702 and try_new_r a b n M r r_list steps =
704 val r = find_new_r a b r;
705 val r_list = r_list @ [r];
706 val g_r = gcd_poly' (order (eval a (n - 2) r))
707 (order (eval b (n - 2) r)) (n - 1) 0
708 val steps = steps @ [g_r];
710 (writeln ("try_new_r calls HENSEL-poly");
711 HENSEL_lifting a b n M 1 r r_list steps g_r ( cero_poly n ) [1]
714 (* fn: poly -> poly -> int ->
715 int -> int -> int -> int list -> poly list ->
716 | poly -> poly -> int list -> poly*)
717 and HENSEL_lifting a b n M m r r_list steps g g_n mult =
718 ((*begin function body*)
719 writeln ("HENSEL-poly n = " ^ PolyML.makestring n ^
720 ", M = " ^ PolyML.makestring M ^ ", m = " ^ PolyML.makestring m ^
721 ", r = " ^ PolyML.makestring r ^ ", r_list = " ^ PolyML.makestring r_list^
722 ", steps = " ^ PolyML.makestring steps ^ ", g = " ^ PolyML.makestring g ^
723 ", g_n = " ^ PolyML.makestring g_n ^ ", mult = " ^ PolyML.makestring mult);
724 if m > M then if g_n %%|%% a andalso g_n %%|%% b then
725 (writeln ("HENSEL-poly ------------------------------> " ^ PolyML.makestring g_n);
729 try_new_r a b n M r r_list steps
732 val r = find_new_r a b r;
733 val r_list = r_list @ [r];
734 val g_r = gcd_poly' (order (eval a (n - 2) r))
735 (order (eval b (n - 2) r)) (n - 1) 0
737 if lex_ord (lmonom g) (lmonom g_r)
739 (writeln ("recurs.call 1 HENSEL-poly");
740 HENSEL_lifting a b n M 1 r r_list steps g g_n mult
742 else if (lexp g) = (lexp g_r)
744 let val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
746 if (nth steps (length steps - 1)) = (cero_poly (n - 1))
748 (writeln ("recurs.call 2 HENSEL-poly");
749 HENSEL_lifting a b n M (M + 1) r r_list steps g_r g_n new
752 (writeln ("recurs.call 3 HENSEL-poly");
753 HENSEL_lifting a b n M (m + 1) r r_list steps g_r g_n new
756 else (* \<not> lex_ord (lmonom g) (lmonom g_r) *)
757 (writeln ("recurs.call 4 HENSEL-poly");
758 HENSEL_lifting a b n M (m + 1) r r_list steps g g_n mult
761 )(*end function body*)
763 (* gcd_poly :: poly \<Rightarrow> poly \<Rightarrow> poly
764 gcd_poly a b = ((a', b'), c)
765 assumes a \<noteq> [] \<and> b \<noteq> [] \<and>
766 yields a = a' *\<^isub>p c \<and> b = b' *\<^isub>p c \<and>
767 (\<forall>c'. (c' dvd\<^isub>p a \<and> c' dvd\<^isub>p b) \<longrightarrow> c' \<le>\<^isub>p c) *)
768 fun gcd_poly (a as (_, es)::_ : poly) b =
769 let val c = gcd_poly' a b (length es) 0
770 val (a': poly, _) = a %%/%% c
771 val (b': poly, _) = b %%/%% c
775 section {* example from the last mail *}
776 ML {* (* test, see text below *)
777 val a = [(1,[2, 0]), (~1,[0, 2])];
778 val b = [(1,[2, 0]), (~1,[1, 1])] ;
779 val ((a', b'), c) = gcd_poly a b;
781 a' = [(~1, [1, 0]), (~1, [0, 1])];
783 c = [(~1, [1, 0]), (1, [0, 1])];
784 (* checking the postcondition: *)
785 a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)
786 (* \<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 *)
789 text {* last mail to Tobias Nipkow:
790 Eine erste Idee, wie die Integration der Diplomarbeit für einen Benutzer
791 von Isabelle aussehen könnte, wäre zum Beispiel im
794 assumes asm3: "x^2 - x*y \<noteq> 0" and asm4: "x \<noteq> 0"
795 shows "(x^2 - y^2) / (x^2 - x*y) = (x + y) / (x::real)"
796 apply (insert asm3 asm4)
802 asm1: "(x^2 - y^2) = (x + y) * (x - y)" and asm2: "x^2 - x*y = x * (x - y)"
804 im Hintergrund automatisch zu erzeugen (mit der Garantie, dass "(x - y)" der GCD ist)
805 und sie dem Simplifier (für die Rule nonzero_mult_divide_mult_cancel_right) zur
806 Verfügung zu stellen, sodass anstelle von "sorry" einfach "done" stehen kann.
807 Und weiters wäre eventuell asm3 zu "x - y \<noteq> 0" zu vereinfachen, eine
808 Rewriteorder zum Herstellen einer Normalform festzulegen, etc.