1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 10:41:57 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 12:40:15 2013 +0200
1.3 @@ -382,7 +382,10 @@
1.4
1.5 argumentns "a b d M p" named according to [1] p.93 *)
1.6 fun HENSEL_lifting_up a b d M p =
1.7 - let
1.8 + let
1.9 +val _ = writeln ("HENSEL-univ d = " ^ PolyML.makestring d ^
1.10 +", M = " ^ PolyML.makestring (2 * d * LANDAU_MIGNOTTE_bound a b) ^
1.11 +", p = " ^ PolyML.makestring p)
1.12 val p = p next_prime_not_dvd d
1.13 val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
1.14 in
1.15 @@ -391,9 +394,7 @@
1.16 val g = primitive_up (try_new_prime_up a b d M p g_p p)
1.17 in
1.18 if (g %|% a) andalso (g %|% b) then
1.19 -(writeln ("HENSEL-univ M = " ^ PolyML.makestring M ^ ", p = " ^ PolyML.makestring p ^
1.20 -": " ^ PolyML.makestring g ^ " %|% " ^ PolyML.makestring a ^
1.21 -" /\\ " ^ PolyML.makestring g ^ " %|% " ^ PolyML.makestring b);
1.22 +(writeln ("HENSEL-univ -------------------> " ^ PolyML.makestring g);
1.23 g:unipoly
1.24 )
1.25 else HENSEL_lifting_up a b d M p
1.26 @@ -493,8 +494,8 @@
1.27
1.28 subsection {* evaluation, translation uni--multivariate, etc *}
1.29 ML {*
1.30 - (* evaluate a polynomial in a certain variable and remove this from the exp-list *)
1.31 - fun eval (p: poly) var value =
1.32 + (* evaluate a polynomial in a certain variable and remove this var from the exp-list *)
1.33 + fun eval (p: poly) var value = (* TODO use "map" instead*)
1.34 let
1.35 fun eval [] _ _ new = order new
1.36 | eval ((c, es)::ps) var value new =
1.37 @@ -679,20 +680,22 @@
1.38
1.39 subsection {* gcd_poly algorithm, code for [1] p.95 *}
1.40 ML {*
1.41 - (* naming of n, M, m, r, ... follows Winkler p. 95 *)
1.42 + (* naming of n, M, m, r, ... follows Winkler p. 95
1.43 + n: amount of variables
1.44 + r: *)
1.45 fun gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r =
1.46 if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else
1.47 - if n - 1 = 0
1.48 + if n > 1
1.49 then
1.50 + let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
1.51 + in try_new_r a b n M r [] [] end
1.52 + else
1.53 let val (a', b') = (multi_to_uni a, multi_to_uni b)
1.54 in
1.55 (* a b are made primitive and this is compensated by multiplying with gcd of divisors *)
1.56 uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %*
1.57 (abs (Integer.gcd (gcd_coeff a') (gcd_coeff b'))))
1.58 end
1.59 - else
1.60 - let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
1.61 - in try_new_r a b n M r [] [] end
1.62
1.63 (* fn: poly -> poly -> int ->
1.64 int -> int -> int list -> poly list -> poly*)
1.65 @@ -712,10 +715,14 @@
1.66 int -> int -> int -> int list -> poly list ->
1.67 | poly -> poly -> int list -> poly*)
1.68 and HENSEL_lifting a b n M m r r_list steps g g_n mult =
1.69 +((*begin function body*)
1.70 +writeln ("HENSEL-poly n = " ^ PolyML.makestring n ^
1.71 +", M = " ^ PolyML.makestring M ^ ", m = " ^ PolyML.makestring m ^
1.72 +", r = " ^ PolyML.makestring r ^ ", r_list = " ^ PolyML.makestring r_list^
1.73 +", steps = " ^ PolyML.makestring steps ^ ", g = " ^ PolyML.makestring g ^
1.74 +", g_n = " ^ PolyML.makestring g_n ^ ", mult = " ^ PolyML.makestring mult);
1.75 if m > M then if g_n %%|%% a andalso g_n %%|%% b then
1.76 -(writeln ("HENSEL-poly M = " ^ PolyML.makestring M ^ "HENSEL-poly M = " ^ PolyML.makestring M ^
1.77 -": " ^ PolyML.makestring g_n ^ " %%|%% " ^ PolyML.makestring a ^
1.78 -" /\\ " ^ PolyML.makestring g_n ^ " %%|%% " ^ PolyML.makestring b);
1.79 +(writeln ("HENSEL-poly ------------------------------> " ^ PolyML.makestring g_n);
1.80 g_n
1.81 )
1.82 else
1.83 @@ -750,7 +757,8 @@
1.84 (writeln ("recurs.call 4 HENSEL-poly");
1.85 HENSEL_lifting a b n M (m + 1) r r_list steps g g_n mult
1.86 )
1.87 - end
1.88 + end
1.89 +)(*end function body*)
1.90
1.91 (* gcd_poly :: poly \<Rightarrow> poly \<Rightarrow> poly
1.92 gcd_poly a b = ((a', b'), c)
2.1 --- a/test/Tools/isac/Test_Some2.thy Fri May 24 10:41:57 2013 +0200
2.2 +++ b/test/Tools/isac/Test_Some2.thy Fri May 24 12:40:15 2013 +0200
2.3 @@ -10,31 +10,23 @@
2.4 ML {*
2.5 *}
2.6 ML {*
2.7 -if gcd_poly
2.8 +gcd_poly
2.9 [(~3,[2,0]),(1,[5,0]),(3,[0,1]),(~6,[1,1]),(~1,[3,1]),(2,[4,1]),(1,[3,2]),(~1,[1,3]),(2,[2,3])]
2.10 - [(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])]
2.11 -= (([(~3, [0, 0]), (1, [3, 0]), (1, [1, 2])], [(2, [0, 0]), (~1, [1, 1]), (1, [0, 2])]),
2.12 - [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])]) then () else error "gcd_poly ex1 changed";
2.13 + [(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])];
2.14 *}
2.15 ML {*
2.16 -if gcd_poly
2.17 +gcd_poly
2.18 [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
2.19 - [(1,[0,0,0]),(1,[1,1,0])]
2.20 - = (([(~1, [0,0,0]), (1, [0,1,1])], [(1, [0,0,0])]), [(1, [0,0,0]), (1, [1,1,0])])
2.21 -then () else error "gcd_poly ex2 changed";
2.22 + [(1,[0,0,0]),(1,[1,1,0])];
2.23 *}
2.24 ML {*
2.25 -val a = uni_to_multi [~18, ~15, ~20, 12, 20, ~13, 2];
2.26 -val b = uni_to_multi [8, 28, 22, ~11, ~14, 1, 2];
2.27 -val ((a', b'), c) = gcd_poly a b;
2.28 +gcd_poly (uni_to_multi [~18, ~15, ~20, 12, 20, ~13, 2]) (uni_to_multi [8, 28, 22, ~11, ~14, 1, 2]);
2.29 *}
2.30 ML {*
2.31 -val a = [(1,[2, 0]), (~1,[0, 2])]; (* x^2 - y^2 *)
2.32 -val b = [(1,[2, 0]), (~1,[1, 1])]; (* x^2 - x y *)
2.33 -val ((a', b'), c) = gcd_poly a b;
2.34 + (* x^2 - y^2 *) (* x^2 - x y *)
2.35 +gcd_poly [(1,[2, 0]), (~1,[0, 2])] [(1,[2, 0]), (~1,[1, 1])];
2.36 *}
2.37 ML {*
2.38 -try_new_r
2.39 *}
2.40 ML {*
2.41 *}