659 subsection {* gcd_poly algorithm, code for [1] p.95 *} |
659 subsection {* gcd_poly algorithm, code for [1] p.95 *} |
660 ML {* |
660 ML {* |
661 (* naming of n, M, m, r, ... follows Winkler p. 95 |
661 (* naming of n, M, m, r, ... follows Winkler p. 95 |
662 n: amount of variables |
662 n: amount of variables |
663 r: *) |
663 r: *) |
664 fun gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r = |
664 fun gcd_poly' [(c1, es1)] [(c2, es2)] _ _ = [(Integer.gcd c1 c2, map2 Integer.min es1 es2)] |
665 if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else |
665 | gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r = |
666 if n > 1 |
666 if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else |
667 then |
667 if n > 1 |
668 let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2)); |
668 then |
669 in try_new_r a b n M r [] [] end |
669 let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2)); |
670 else |
670 in try_new_r a b n M r [] [] end |
671 let val (a', b') = (multi_to_uni a, multi_to_uni b) |
671 else |
672 in |
672 let val (a', b') = (multi_to_uni a, multi_to_uni b) |
673 (* a b are made primitive and this is compensated by multiplying with gcd of divisors *) |
673 in |
674 uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %* |
674 (* a b are made primitive and this is compensated by multiplying with gcd of divisors *) |
675 (abs (Integer.gcd (gcd_coeff a') (gcd_coeff b')))) |
675 uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %* |
676 end |
676 (abs (Integer.gcd (gcd_coeff a') (gcd_coeff b')))) |
|
677 end |
677 |
678 |
678 (* fn: poly -> poly -> int -> |
679 (* fn: poly -> poly -> int -> |
679 int -> int -> int list -> poly list -> poly*) |
680 int -> int -> int list -> poly list -> poly*) |
680 and try_new_r a b n M r r_list steps = |
681 and try_new_r a b n M r r_list steps = |
681 let |
682 let |