1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Thu May 23 12:29:02 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 10:41:57 2013 +0200
1.3 @@ -390,7 +390,13 @@
1.4 let
1.5 val g = primitive_up (try_new_prime_up a b d M p g_p p)
1.6 in
1.7 - if (g %|% a) andalso (g %|% b) then g:unipoly else HENSEL_lifting_up a b d M p
1.8 + if (g %|% a) andalso (g %|% b) then
1.9 +(writeln ("HENSEL-univ M = " ^ PolyML.makestring M ^ ", p = " ^ PolyML.makestring p ^
1.10 +": " ^ PolyML.makestring g ^ " %|% " ^ PolyML.makestring a ^
1.11 +" /\\ " ^ PolyML.makestring g ^ " %|% " ^ PolyML.makestring b);
1.12 + g:unipoly
1.13 +)
1.14 + else HENSEL_lifting_up a b d M p
1.15 end
1.16 end
1.17
1.18 @@ -403,7 +409,9 @@
1.19 let val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))
1.20 in
1.21 if a = b then a else
1.22 +(writeln ("gcd_up calls HENSEL-univ");
1.23 HENSEL_lifting_up a b d (2 * d * LANDAU_MIGNOTTE_bound a b) 1
1.24 +)
1.25 end;
1.26 *}
1.27
1.28 @@ -672,7 +680,7 @@
1.29 subsection {* gcd_poly algorithm, code for [1] p.95 *}
1.30 ML {*
1.31 (* naming of n, M, m, r, ... follows Winkler p. 95 *)
1.32 - fun gcd_poly' (a as (c1, es1)::ps1 : poly) (b as (c2, es2)::ps2 : poly) n r =
1.33 + fun gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r =
1.34 if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else
1.35 if n - 1 = 0
1.36 then
1.37 @@ -685,7 +693,9 @@
1.38 else
1.39 let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2));
1.40 in try_new_r a b n M r [] [] end
1.41 -
1.42 +
1.43 + (* fn: poly -> poly -> int ->
1.44 + int -> int -> int list -> poly list -> poly*)
1.45 and try_new_r a b n M r r_list steps =
1.46 let
1.47 val r = find_new_r a b r;
1.48 @@ -693,10 +703,22 @@
1.49 val g_r = gcd_poly' (order (eval a (n - 2) r))
1.50 (order (eval b (n - 2) r)) (n - 1) 0
1.51 val steps = steps @ [g_r];
1.52 - in HENSEL_lifting a b n M 1 r r_list steps g_r ( cero_poly n ) [1] end
1.53 -
1.54 - and HENSEL_lifting a b n M m r r_list steps g g_n mult =
1.55 - if m > M then if g_n %%|%% a andalso g_n %%|%% b then g_n else
1.56 + in
1.57 +(writeln ("try_new_r calls HENSEL-poly");
1.58 + HENSEL_lifting a b n M 1 r r_list steps g_r ( cero_poly n ) [1]
1.59 +)
1.60 + end
1.61 + (* fn: poly -> poly -> int ->
1.62 + int -> int -> int -> int list -> poly list ->
1.63 + | poly -> poly -> int list -> poly*)
1.64 + and HENSEL_lifting a b n M m r r_list steps g g_n mult =
1.65 + if m > M then if g_n %%|%% a andalso g_n %%|%% b then
1.66 +(writeln ("HENSEL-poly M = " ^ PolyML.makestring M ^ "HENSEL-poly M = " ^ PolyML.makestring M ^
1.67 +": " ^ PolyML.makestring g_n ^ " %%|%% " ^ PolyML.makestring a ^
1.68 +" /\\ " ^ PolyML.makestring g_n ^ " %%|%% " ^ PolyML.makestring b);
1.69 + g_n
1.70 +)
1.71 + else
1.72 try_new_r a b n M r r_list steps
1.73 else
1.74 let
1.75 @@ -706,16 +728,28 @@
1.76 (order (eval b (n - 2) r)) (n - 1) 0
1.77 in
1.78 if lex_ord (lmonom g) (lmonom g_r)
1.79 - then HENSEL_lifting a b n M 1 r r_list steps g g_n mult
1.80 + then
1.81 +(writeln ("recurs.call 1 HENSEL-poly");
1.82 + HENSEL_lifting a b n M 1 r r_list steps g g_n mult
1.83 +)
1.84 else if (lexp g) = (lexp g_r)
1.85 - then let
1.86 - val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
1.87 - in
1.88 - if (nth steps (length steps - 1)) = (cero_poly (n - 1))
1.89 - then HENSEL_lifting a b n M (M + 1) r r_list steps g_r g_n new
1.90 - else HENSEL_lifting a b n M (m + 1) r r_list steps g_r g_n new
1.91 - end
1.92 - else HENSEL_lifting a b n M (m + 1) r r_list steps g g_n mult
1.93 + then
1.94 + let val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
1.95 + in
1.96 + if (nth steps (length steps - 1)) = (cero_poly (n - 1))
1.97 + then
1.98 +(writeln ("recurs.call 2 HENSEL-poly");
1.99 + HENSEL_lifting a b n M (M + 1) r r_list steps g_r g_n new
1.100 +)
1.101 + else
1.102 +(writeln ("recurs.call 3 HENSEL-poly");
1.103 + HENSEL_lifting a b n M (m + 1) r r_list steps g_r g_n new
1.104 +)
1.105 + end
1.106 + else (* \<not> lex_ord (lmonom g) (lmonom g_r) *)
1.107 +(writeln ("recurs.call 4 HENSEL-poly");
1.108 + HENSEL_lifting a b n M (m + 1) r r_list steps g g_n mult
1.109 +)
1.110 end
1.111
1.112 (* gcd_poly :: poly \<Rightarrow> poly \<Rightarrow> poly
1.113 @@ -736,10 +770,12 @@
1.114 val b = [(1,[2, 0]), (~1,[1, 1])] ;
1.115 val ((a', b'), c) = gcd_poly a b;
1.116
1.117 - a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)
1.118 - c = [(~1, [1, 0]), (1, [0, 1])];
1.119 a' = [(~1, [1, 0]), (~1, [0, 1])];
1.120 b' = [(~1, [1, 0])];
1.121 + c = [(~1, [1, 0]), (1, [0, 1])];
1.122 + (* checking the postcondition: *)
1.123 + a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)
1.124 + (* \<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 *)
1.125 *}
1.126
1.127 text {* last mail to Tobias Nipkow:
2.1 --- a/test/Tools/isac/Knowledge/gcd_poly.sml Thu May 23 12:29:02 2013 +0200
2.2 +++ b/test/Tools/isac/Knowledge/gcd_poly.sml Fri May 24 10:41:57 2013 +0200
2.3 @@ -696,24 +696,50 @@
2.4 "----------- fun gcd_poly ------------------------------";
2.5 "----------- fun gcd_poly ------------------------------";
2.6 if gcd_poly
2.7 -[(~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.8 -[(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.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 + (* ~3 x^2 + x^5 + 3 y - 6 x y - x^3 y + 2 x^4 y + x^3 y^2 - x y^3 + 2 x^2 y^3*)
2.11 + [(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.12 + (* 2 x^2 - 2 y + 4 x y - x^3 y + x y^2 - x^2 y^2 - y^3 + 2 x y^3*)
2.13 = (([(~3, [0, 0]), (1, [3, 0]), (1, [1, 2])], [(2, [0, 0]), (~1, [1, 1]), (1, [0, 2])]),
2.14 - [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])]) then () else error
2.15 -("gcd_poly\n" ^
2.16 -"[(~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.17 -"[(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.18 -"= (([(~3, [0, 0]), (1, [3, 0]), (1, [1, 2])], [(2, [0, 0]), (~1, [1, 1]), (1, [0, 2])]),"^
2.19 -" [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])]) changed");
2.20 -(* -xy +xy^2z+yz - 1*)(* xy +1*) (*=*) (*xy -1*)
2.21 + (*( ~3 + x^3 + x y^2 ,, 2 - x y + y^2 ),,*)
2.22 + [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])]) then () else error "gcd_poly ex1 changed";
2.23 + (* x^2 - y + 2 x y*)
2.24
2.25 if gcd_poly
2.26 -[(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
2.27 -[(1,[0,0,0]),(1,[1,1,0])]
2.28 -= (([(~1, [0, 0, 0]), (1, [0, 1, 1])], [(1, [0, 0, 0])]), [(1, [0, 0, 0]), (1, [1, 1, 0])]) then ()
2.29 -else error
2.30 -"gcd_poly [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] "
2.31 -"=(([(~1, [0, 0, 0]), (1, [0, 1, 1])], [(1, [0, 0, 0])]), [(1, [0, 0, 0]), (1, [1, 1, 0])])changed";
2.32 + [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
2.33 + (* ~1 + y z + x y^2 z - x y *)
2.34 + [(1,[0,0,0]),(1,[1,1,0])]
2.35 + (* 1 + x y *)
2.36 + = (([(~1, [0,0,0]), (1, [0,1,1])], [(1, [0,0,0])]), [(1, [0,0,0]), (1, [1,1,0])])
2.37 + (* ( ~1 + y z ,, 1 ),, 1 + x y *)
2.38 +then () else error "gcd_poly ex2 changed";
2.39 +
2.40 +(* example from unipoly: *)
2.41 +val a = uni_to_multi [~18, ~15, ~20, 12, 20, ~13, 2];
2.42 +val b = uni_to_multi [8, 28, 22, ~11, ~14, 1, 2];
2.43 +val ((a', b'), c) = gcd_poly a b;
2.44 +
2.45 +a' = [(9, [0]), (3, [1]), (13, [2]), (~11, [3]), (2, [4])];
2.46 +b' = [(~4, [0]), (~12, [1]), (~7, [2]), (3, [3]), (2, [4])];
2.47 +c = [(~2, [0]), (~1, [1]), (1, [2])];
2.48 +(* checking the postcondition: *)
2.49 +if a = (a' %%*%% c) andalso b = (b' %%*%% c) then () else error "gcd_poly ex-unipoly changed";
2.50 +(* \<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 *)
2.51 +
2.52 +val cu = gcd_up [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2];
2.53 +if cu = multi_to_uni c then () else error "gcd_up <> gcd_poly";
2.54 +
2.55 +(* example from mail to Tobias Nipkos *)
2.56 +val a = [(1,[2, 0]), (~1,[0, 2])]; (* x^2 - y^2 *)
2.57 +val b = [(1,[2, 0]), (~1,[1, 1])]; (* x^2 - x y *)
2.58 +val ((a', b'), c) = gcd_poly a b;
2.59 +
2.60 +a' = [(~1, [1, 0]), (~1, [0, 1])]; (* -x - y NOT NICE !!!*)
2.61 +b' = [(~1, [1, 0])]; (* -x NOT NICE !!!*)
2.62 +c = [(~1, [1, 0]), (1, [0, 1])]; (* -x + y *)
2.63 +(* checking the postcondition: *)
2.64 +if a = (a' %%*%% c) andalso b = (b' %%*%% c) then () else error "gcd_poly mail changed";
2.65 +(* \<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 *)
2.66
2.67 " ========== END ======================================= ";
2.68
3.1 --- a/test/Tools/isac/Test_Some2.thy Thu May 23 12:29:02 2013 +0200
3.2 +++ b/test/Tools/isac/Test_Some2.thy Fri May 24 10:41:57 2013 +0200
3.3 @@ -6,25 +6,35 @@
3.4 ML_file "/usr/local/isabisac/test/Tools/isac/Knowledge/gcd_poly.sml"
3.5
3.6 ML {*
3.7 -gcd_poly
3.8 -[(~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])]
3.9 -[(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])]
3.10 -
3.11 -*}
3.12 -ML {*
3.13 -gcd_poly
3.14 -[(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
3.15 -[(1,[0,0,0]),(1,[1,1,0])]
3.16 *}
3.17 ML {*
3.18 *}
3.19 ML {*
3.20 +if gcd_poly
3.21 + [(~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])]
3.22 + [(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])]
3.23 += (([(~3, [0, 0]), (1, [3, 0]), (1, [1, 2])], [(2, [0, 0]), (~1, [1, 1]), (1, [0, 2])]),
3.24 + [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])]) then () else error "gcd_poly ex1 changed";
3.25 *}
3.26 ML {*
3.27 +if gcd_poly
3.28 + [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
3.29 + [(1,[0,0,0]),(1,[1,1,0])]
3.30 + = (([(~1, [0,0,0]), (1, [0,1,1])], [(1, [0,0,0])]), [(1, [0,0,0]), (1, [1,1,0])])
3.31 +then () else error "gcd_poly ex2 changed";
3.32 *}
3.33 ML {*
3.34 +val a = uni_to_multi [~18, ~15, ~20, 12, 20, ~13, 2];
3.35 +val b = uni_to_multi [8, 28, 22, ~11, ~14, 1, 2];
3.36 +val ((a', b'), c) = gcd_poly a b;
3.37 *}
3.38 ML {*
3.39 +val a = [(1,[2, 0]), (~1,[0, 2])]; (* x^2 - y^2 *)
3.40 +val b = [(1,[2, 0]), (~1,[1, 1])]; (* x^2 - x y *)
3.41 +val ((a', b'), c) = gcd_poly a b;
3.42 +*}
3.43 +ML {*
3.44 +try_new_r
3.45 *}
3.46 ML {*
3.47 *}