1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 14:25:34 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 16:50:31 2013 +0200
1.3 @@ -97,6 +97,9 @@
1.4 fun drop_hd_zeros (0 :: ps) = drop_hd_zeros ps
1.5 | drop_hd_zeros p = p
1.6 fun drop_tl_zeros p = (rev o drop_hd_zeros o rev) p;
1.7 +
1.8 + (* swaps the arguments of a function *)
1.9 + fun swapf f a b = f b a
1.10 *}
1.11
1.12 subsection {* modulo calculations for integers *}
1.13 @@ -225,9 +228,7 @@
1.14
1.15 (* scalar divison *)
1.16 infix %/
1.17 - fun (p: unipoly) %/ n =
1.18 - let fun swapf f a b = f b a (* swaps the arguments of a function *)
1.19 - in map (swapf (curry op div2) n) p : unipoly end;
1.20 + fun (p: unipoly) %/ n = map (swapf (curry op div2) n) p : unipoly
1.21
1.22 (* minus of polys *)
1.23 infix %-%
1.24 @@ -432,10 +433,10 @@
1.25 fun lexp (p: poly) = (snd o hd o rev) p
1.26 fun lmonom (p: poly) = (hd o rev) p
1.27
1.28 - fun cero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;
1.29 + fun zero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;
1.30
1.31 (* drop leading coefficients equal 0 TODO: signature?*)
1.32 - fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then cero_poly (length es) else p
1.33 + fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then zero_poly (length es) else p
1.34 | drop_lc0 p = if lcoeff p = 0 then drop_lc0 (drop_last p) else p
1.35
1.36 (* gcd of coefficients WN find right location in file *)
1.37 @@ -443,7 +444,7 @@
1.38
1.39 (* given p = 3x+4z = [(3,[1,0]),(4,[0,1])] with x<z
1.40 add_variable p 2 => [(3,[1,0,0]),(4,[0,0,1])] with x<y<z *)
1.41 - fun add_variable (p: poly) var = map (apsnd (nth_ins var 0)) p : poly
1.42 + fun add_variable pos (p: poly) = map (apsnd (nth_ins pos 0)) p : poly
1.43 *}
1.44
1.45 subsection {* ordering and other auxiliary funcions *}
1.46 @@ -464,7 +465,7 @@
1.47 (* add monomials in ordered polynomal *)
1.48 fun add_monoms (p: poly) =
1.49 let
1.50 - fun add [(0, es)] [] = cero_poly (length es)
1.51 + fun add [(0, es)] [] = zero_poly (length es)
1.52 | add [(0, _)] (new: poly) = new
1.53 | add [m] (new: poly) = new @ [m]
1.54 | add ((c1, es1) :: (c2, es2) :: ms) (new: poly) =
1.55 @@ -486,14 +487,11 @@
1.56 subsection {* evaluation, translation uni--multivariate, etc *}
1.57 ML {*
1.58 (* evaluate a polynomial in a certain variable and remove this var from the exp-list *)
1.59 - fun eval (p: poly) var value = (* TODO use "map" instead*)
1.60 - let
1.61 - fun eval [] _ _ new = order new
1.62 - | eval ((c, es)::ps) var value new =
1.63 - eval ps var value (new @ [((Integer.pow (nth es var) value) * c, nth_drop var es)]);
1.64 - in eval p var value [] end
1.65 + fun eval_monom var value ((c, es): monom) =
1.66 + (c * (Integer.pow (nth es var) value), nth_drop var es) : monom
1.67 + fun eval (p: poly) var value = order (map (eval_monom var value) p)
1.68
1.69 - (* transform a multivariate to a univariate polynomial *)
1.70 + (* transform a multivariate to a univariate polynomial TODO map?*)
1.71 fun multi_to_uni (p as (_, es)::_ : poly) =
1.72 if length es = 1
1.73 then
1.74 @@ -505,7 +503,7 @@
1.75 in transfer (order p) [] end
1.76 else error "Polynom has more than one variable!";
1.77
1.78 - (* transform a univariate to a multivariate polynomial *)
1.79 + (* transform a univariate to a multivariate polynomial TODO map *)
1.80 fun uni_to_multi (up: unipoly) =
1.81 let fun transfer ([]: unipoly) (mp: poly) (_: int) = add_monoms mp
1.82 | transfer up mp var =
1.83 @@ -513,7 +511,7 @@
1.84 in transfer up [] 0 end
1.85
1.86 (* multiply a poly with a variable represented by a unipoly (at position var)
1.87 - e.g new z: (3x + 2y)*(z - 4) TODO fold / map ? *)
1.88 + e.g new z: (3x + 2y)*(z - 4) TODO mp %%*%% (uni_to_multi up) !TODO equal es ! *)
1.89 fun mult_with_new_var ([]: poly) (_: unipoly) (_: int) = []
1.90 | mult_with_new_var mp up var =
1.91 let
1.92 @@ -581,19 +579,6 @@
1.93 fun p1 %%-%% (p2: poly) = p1 %%+%% (map (apfst (Integer.mult ~1)) p2)
1.94
1.95 infix %%*%%
1.96 - (* the same algorithms as done by hand
1.97 - fun (p1: poly) %%*%% (p2: poly) =
1.98 - let
1.99 - fun mult ([]: poly) (_: poly) (_: poly) (new: poly) = order new
1.100 - | mult p1 [] regular_p2 new = mult (tl p1) regular_p2 regular_p2 new
1.101 - | mult (p1 as (c1, es1)::_) ((c2, es2)::ps2) regular_p2 (new: poly) =
1.102 - mult p1 ps2 regular_p2 (new @ [(c1 * c2, map2 Integer.add es1 es2)])
1.103 - in
1.104 - if length p1 > length p2
1.105 - then mult p1 p2 p2 []
1.106 - else mult p2 p1 p1 []
1.107 - end
1.108 - *)
1.109 fun mult_monom ((c1, es1): monom) ((c2, es2): monom) =
1.110 (c1 * c2, map2 Integer.add es1 es2): monom;
1.111 fun mult_monoms (ms: poly) m = map (mult_monom m) ms : poly
1.112 @@ -623,7 +608,7 @@
1.113 let
1.114 val p1 as (c, es)::_ = drop_lc0 (order p1);
1.115 val p2 = drop_lc0 (order p2);
1.116 - val [c] = cero_poly (length es);
1.117 + val [c] = zero_poly (length es);
1.118 val d = (replicate (List.length p1 - List.length p2) c) @ p2 ;
1.119 val quot = [(lcoeff p1 div2 lcoeff d, lexp p1 %-% lexp p2)];
1.120 val rest = drop_lc0 (p1 %%-%% (d %%*%% quot));
1.121 @@ -637,12 +622,12 @@
1.122
1.123 subsection {* Newton interpolation *}
1.124 ML{* (* first step *)
1.125 - fun newton_first (x: int list) (f: poly list) (ord: int) =
1.126 + fun newton_first (x: int list) (f: poly list) (pos: int) =
1.127 let
1.128 val new_vals = [(hd x) * ~1, 1];
1.129 - val p = (add_variable (hd f) ord) %%+%%
1.130 + val p = (add_variable pos (hd f)) %%+%%
1.131 ((mult_with_new_var
1.132 - (((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x))) new_vals ord);
1.133 + (((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x))) new_vals pos);
1.134 val steps = [((nth f 1) %%-%% (hd f)) %%/ ((nth x 1) - (hd x))]
1.135 in (p, new_vals, steps) end
1.136
1.137 @@ -653,10 +638,10 @@
1.138 [((last new_steps) %%-%% (hd steps)) %%/ ((last x) - (nth x (length x - 3)))])
1.139 (nth_drop (length x -2) x)
1.140
1.141 - fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) ord =
1.142 + fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) pos =
1.143 if length x = 2
1.144 then
1.145 - let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] ord
1.146 + let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] pos
1.147 in (p', new_vals, steps) end
1.148 else
1.149 let
1.150 @@ -665,7 +650,7 @@
1.151 val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/
1.152 ((last x) - (nth x (length x - 2)))];
1.153 val steps = next_step steps new_steps x;
1.154 - val p' = p %%+%% (mult_with_new_var (last steps) new_vals ord);
1.155 + val p' = p %%+%% (mult_with_new_var (last steps) new_vals pos);
1.156 in (order p', new_vals, steps) end
1.157 *}
1.158
1.159 @@ -697,7 +682,7 @@
1.160 val g_r = gcd_poly' (order (eval a (n - 2) r))
1.161 (order (eval b (n - 2) r)) (n - 1) 0
1.162 val steps = steps @ [g_r];
1.163 - in HENSEL_lifting a b n M 1 r r_list steps g_r ( cero_poly n ) [1] end
1.164 + in HENSEL_lifting a b n M 1 r r_list steps g_r ( zero_poly n ) [1] end
1.165 (* fn: poly -> poly -> int ->
1.166 int -> int -> int -> int list -> poly list ->
1.167 | poly -> poly -> int list -> poly*)
1.168 @@ -717,7 +702,7 @@
1.169 then
1.170 let val (g_n, new, steps) = NEWTON r_list [g, g_r] steps mult g_n (n - 2)
1.171 in
1.172 - if (nth steps (length steps - 1)) = (cero_poly (n - 1))
1.173 + if (nth steps (length steps - 1)) = (zero_poly (n - 1))
1.174 then HENSEL_lifting a b n M (M + 1) r r_list steps g_r g_n new
1.175 else HENSEL_lifting a b n M (m + 1) r r_list steps g_r g_n new
1.176 end
2.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri May 24 14:25:34 2013 +0200
2.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri May 24 16:50:31 2013 +0200
2.3 @@ -728,7 +728,6 @@
2.4 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
2.5
2.6 argumentns "a b d M P g p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
2.7 -(*declare [[show_types]]*)
2.8 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
2.9 where
2.10 "try_new_prime_up a b d M P g p =
2.11 @@ -747,8 +746,7 @@
2.12 let
2.13 P = P * p;
2.14 g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
2.15 - in
2.16 - try_new_prime_up a b d M P g p)"
2.17 + in try_new_prime_up a b d M P g p)"
2.18 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
2.19 termination try_new_prime_up (*by lexicographic_order LOOPS +PROOF:? / by size_change LOOPS*)
2.20 (*apply (relation "measures
3.1 --- a/test/Tools/isac/Knowledge/gcd_poly.sml Fri May 24 14:25:34 2013 +0200
3.2 +++ b/test/Tools/isac/Knowledge/gcd_poly.sml Fri May 24 16:50:31 2013 +0200
3.3 @@ -39,7 +39,7 @@
3.4 "----------- fun lcoeff ---------------------------------";
3.5 "----------- fun lexp -----------------------------------";
3.6 "----------- fun add_variable ---------------------------";
3.7 -"----------- fun cero_poly ------------------------------";
3.8 +"----------- fun zero_poly ------------------------------";
3.9 "----------- drop_lc0 -----------------------------------";
3.10 "----------- fun add_monoms -----------------------------";
3.11 "----------- fun lex_ord --------------------------------";
3.12 @@ -403,23 +403,23 @@
3.13 "----------- fun add_variable ---------------------------";
3.14 "----------- fun add_variable ---------------------------";
3.15 "----------- fun add_variable ---------------------------";
3.16 -if add_variable
3.17 - [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] then () else
3.18 +if add_variable 0
3.19 + [(1,[1,2]),(2,[2,3])] = [(1, [0, 1, 2]), (2, [0, 2, 3])] then () else
3.20 error "add_variable [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] changed";
3.21 -if add_variable
3.22 - [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3])] then () else
3.23 +if add_variable 1
3.24 + [(1,[1,2]),(2,[2,3])] = [(1, [1, 0, 2]), (2, [2, 0, 3])] then () else
3.25 error "add_variable [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3]) changed";
3.26 -if add_variable
3.27 - [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] then () else
3.28 +if add_variable 2
3.29 + [(1,[1,2]),(2,[2,3])] = [(1, [1, 2, 0]), (2, [2, 3, 0])] then () else
3.30 error "add_variable [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] changed";
3.31
3.32 -"----------- fun cero_poly ------------------------------";
3.33 -"----------- fun cero_poly ------------------------------";
3.34 -"----------- fun cero_poly ------------------------------";
3.35 -if cero_poly 1 = [(0, [0])] then () else
3.36 -error "cero_poly 1 = [(0, [0])] changed";
3.37 -if cero_poly 5 = [(0, [0, 0, 0, 0, 0])] then () else
3.38 -error "cero_poly 1 = [(0, [0, 0, 0, 0, 0])] changed";
3.39 +"----------- fun zero_poly ------------------------------";
3.40 +"----------- fun zero_poly ------------------------------";
3.41 +"----------- fun zero_poly ------------------------------";
3.42 +if zero_poly 1 = [(0, [0])] then () else
3.43 +error "zero_poly 1 = [(0, [0])] changed";
3.44 +if zero_poly 5 = [(0, [0, 0, 0, 0, 0])] then () else
3.45 +error "zero_poly 1 = [(0, [0, 0, 0, 0, 0])] changed";
3.46
3.47 "----------- drop_lc0 -----------------------------------";
3.48 "----------- drop_lc0 -----------------------------------";
3.49 @@ -617,10 +617,10 @@
3.50 "----------- fun NEWTON --------------------------------";
3.51 "----------- fun NEWTON --------------------------------";
3.52 "----------- fun NEWTON --------------------------------";
3.53 -if NEWTON [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_poly 2) 1 =
3.54 +if NEWTON [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (zero_poly 2) 1 =
3.55 ([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])
3.56 then () else error
3.57 -"NEWTON [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_poly 2) 1 = "
3.58 +"NEWTON [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (zero_poly 2) 1 = "
3.59 "([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";
3.60 if NEWTON [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 =
3.61 ([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]])
4.1 --- a/test/Tools/isac/Test_Some2.thy Fri May 24 14:25:34 2013 +0200
4.2 +++ b/test/Tools/isac/Test_Some2.thy Fri May 24 16:50:31 2013 +0200
4.3 @@ -10,26 +10,23 @@
4.4 ML {*
4.5 *}
4.6 ML {*
4.7 -gcd_poly
4.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])]
4.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])];
4.10 *}
4.11 ML {*
4.12 -gcd_poly
4.13 - [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])]
4.14 - [(1,[0,0,0]),(1,[1,1,0])];
4.15 -*}
4.16 -ML {*
4.17 -gcd_poly (uni_to_multi [~18, ~15, ~20, 12, 20, ~13, 2]) (uni_to_multi [8, 28, 22, ~11, ~14, 1, 2]);
4.18 -*}
4.19 -ML {*
4.20 - (* x^2 - y^2 *) (* x^2 - x y *)
4.21 -gcd_poly [(1,[2, 0]), (~1,[0, 2])] [(1,[2, 0]), (~1,[1, 1])];
4.22 -*}
4.23 +*}
4.24 ML {*
4.25 *}
4.26 ML {*
4.27 *}
4.28 +ML {*
4.29 +*}
4.30 +ML {*
4.31 +*}
4.32 +ML {*
4.33 +*}
4.34 +ML {*
4.35 +*}
4.36 +ML {*
4.37 +*}
4.38 (*------------------------------------------------------------------------------------------
4.39 "~~~~~ fun , args:"; val () = ();
4.40 "~~~~~ to return val:"; val () = ();