1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Wed Jan 23 17:29:49 2013 +0100
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Thu Jan 24 17:17:03 2013 +0100
1.3 @@ -29,12 +29,13 @@
1.4 ML {* (* 5 div 2 = 2; ~5 div 2 = ~3; ~5 div2 2 = ~2;*)
1.5 infix div2
1.6 fun a div2 b =
1.7 - if (a div b) < 0 then abs a div (abs b) * ~1
1.8 + if a div b < 0 then (abs a) div (abs b) * ~1
1.9 else a div b;
1.10 *}
1.11
1.12 subsection {* modulo calculations for integers *}
1.13 ML {*
1.14 + (* ??? *)
1.15 fun mod_inv r m =
1.16 let
1.17 fun modi (r, rold, m, rinv) =
1.18 @@ -46,8 +47,10 @@
1.19 else 0
1.20 in modi (r, r, m, 1) end;
1.21
1.22 + (* ??? *)
1.23 fun mod_div a b m = a * (mod_inv b m) mod m;
1.24
1.25 + (* ??? *)
1.26 fun approx_root a =
1.27 let fun root' a n =
1.28 if n * n < a then root' a (n + 1) else n
1.29 @@ -100,49 +103,50 @@
1.30
1.31 subsection {* auxiliary functions for univariate polynomials *}
1.32 ML {*
1.33 + (* leading coefficient *)
1.34 fun lc (uvp: unipoly) =
1.35 if nth uvp (List.length uvp - 1) <> 0
1.36 then nth uvp (List.length uvp - 1)
1.37 else lc (nth_drop (length uvp - 1) uvp);
1.38
1.39 + (* degree *)
1.40 fun deg (uvp: unipoly) =
1.41 if nth uvp (length uvp - 1) <> 0
1.42 then length uvp - 1
1.43 else deg (nth_drop (length uvp - 1) uvp);
1.44
1.45 - fun lc_of_unipoly_not_0 [] = [](* and delete lc=0 WN: id for(multi)poly?*)
1.46 - | lc_of_unipoly_not_0 (uvp: unipoly) =
1.47 - if nth uvp (length uvp -1) =0
1.48 - then lc_of_unipoly_not_0 (nth_drop (length uvp-1) uvp)
1.49 - else uvp;
1.50 + (* drop leading coefficients equal 0 *)
1.51 + fun drop_lc0 [] = []
1.52 + | drop_lc0 (uvp: unipoly) =
1.53 + if nth uvp (length uvp - 1) = 0
1.54 + then drop_lc0 (nth_drop (length uvp - 1) uvp)
1.55 + else uvp;
1.56
1.57 - fun normirt_polynom (poly1: unipoly) (m: int) =
1.58 + (* ??? *)
1.59 + fun normalise (p: unipoly) (m: int) =
1.60 let
1.61 - val poly1 = lc_of_unipoly_not_0 poly1
1.62 - val lc_a = lc poly1;
1.63 - fun normirt poly1 b m lc_a i =
1.64 + val p = drop_lc0 p
1.65 + val lcp = lc p;
1.66 + fun norm nrm i =
1.67 if i = 0
1.68 - then [mod_div (nth poly1 i) lc_a m] @ b
1.69 - else normirt poly1 ([mod_div (nth poly1 i) lc_a m] @ b) m lc_a (i - 1) ;
1.70 + then [mod_div (nth p i) lcp m] @ nrm
1.71 + else norm ([mod_div (nth p i) lcp m] @ nrm) (i - 1) ;
1.72 in
1.73 - if length(poly1)=0
1.74 - then poly1
1.75 - else normirt poly1 [] m lc_a (length(poly1)-1)
1.76 - end
1.77 + if List.length p = 0 then p else norm [] (List.length p - 1)
1.78 + end;
1.79
1.80 + (* scalar multiplication *)
1.81 + infix %*
1.82 + fun (p: unipoly) %* n = map (Integer.mult n) p: unipoly
1.83 +
1.84 + (* scalar divison *)
1.85 + fun swapf f a b = f b a;
1.86 + infix %/
1.87 + fun (p: unipoly) %/ n = map (swapf (curry op div2) n) p: unipoly;
1.88
1.89 - infix %*
1.90 - fun (a: unipoly) %* b = map2 Integer.mult a (replicate (length (a)) b)
1.91 -
1.92 - infix %/
1.93 - fun (poly: unipoly) %/ b = (* =quotient*)
1.94 - let fun division poly b = poly div2 b;
1.95 - in map2 division poly (replicate (length (poly)) b) end
1.96 -
1.97 + (* minus of polys *)
1.98 infix %-%
1.99 - fun (a: unipoly) %-% (b: unipoly) =
1.100 - let fun minus a b = a-b;
1.101 - in map2 minus a b end
1.102 + fun (p1: unipoly) %-% (p2: unipoly) = map2 (curry op -) p1 p2
1.103
1.104 (* if poly2|poly1 *)
1.105 infix %|%
1.106 @@ -152,10 +156,10 @@
1.107 else false
1.108 | (poly2: unipoly) %|% (poly1: unipoly) =
1.109 let
1.110 - val b = (replicate (length poly1 - length(lc_of_unipoly_not_0 poly2)) 0) @
1.111 - lc_of_unipoly_not_0 poly2 ;
1.112 + val b = (replicate (length poly1 - length(drop_lc0 poly2)) 0) @
1.113 + drop_lc0 poly2 ;
1.114 val c = lc poly1 div2 lc b;
1.115 - val rest = lc_of_unipoly_not_0 (poly1 %-% (b %* c));
1.116 + val rest = drop_lc0 (poly1 %-% (b %* c));
1.117 in
1.118 if rest = []
1.119 then true
1.120 @@ -209,20 +213,20 @@
1.121
1.122 infix mod_poly
1.123 fun uvp mod_poly m =
1.124 - let fun mod_poly' uvp m n =
1.125 - if n < length (uvp)
1.126 - then mod_poly' ((nth_drop 0 uvp)@[(nth uvp 0) mod m]) m ( n + 1 )
1.127 - else uvp (*end of mod_poly'*)
1.128 - in
1.129 - mod_poly' uvp m 0 end
1.130 + let
1.131 + fun mod_pol uvp n =
1.132 + if n < List.length uvp
1.133 + then mod_pol ((List.tl uvp) @ [(List.hd uvp) mod m]) (n + 1 )
1.134 + else uvp
1.135 + in mod_pol uvp 0 end
1.136
1.137 fun mod_poly_gcd (upoly1: unipoly) (upoly2: unipoly) (m: int) =
1.138 let
1.139 val moda = upoly1 mod_poly m;
1.140 - val modb = (replicate (length upoly1 - length(lc_of_unipoly_not_0 (upoly2 mod_poly m))) 0)
1.141 - @ (lc_of_unipoly_not_0 (upoly2 mod_poly m)) ;
1.142 + val modb = (replicate (length upoly1 - length(drop_lc0 (upoly2 mod_poly m))) 0)
1.143 + @ (drop_lc0 (upoly2 mod_poly m)) ;
1.144 val c = mod_div (lc moda) (lc modb) m
1.145 - val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) mod_poly m)
1.146 + val rest = drop_lc0 ((moda %-% (modb %* c)) mod_poly m)
1.147 in if rest = []
1.148 then [upoly2, [0]]
1.149 else
1.150 @@ -233,7 +237,7 @@
1.151
1.152 fun mod_gcd_p (poly1: unipoly) (poly2: unipoly) (p: int) =
1.153 let val gcd_p = mod_poly_gcd poly1 poly2 p
1.154 - in normirt_polynom (nth gcd_p 0) p
1.155 + in normalise (nth gcd_p 0) p
1.156 end
1.157
1.158 fun pp [_: int] = [1]
1.159 @@ -249,7 +253,7 @@
1.160 if P > M then g else
1.161 let
1.162 val p = p next_prime_not_dividing d
1.163 - val g_p = poly_centr(((normirt_polynom ( mod_gcd_p a b p) p) %* (d mod p)) mod_poly p) p
1.164 + val g_p = poly_centr(((normalise ( mod_gcd_p a b p) p) %* (d mod p)) mod_poly p) p
1.165 in
1.166 if deg g_p < deg g
1.167 then if (deg g_p) = 0 then [1] else WHILE a b d M p g_p p
1.168 @@ -264,7 +268,7 @@
1.169 fun GOTO2 a b d M p =
1.170 let
1.171 val p = p next_prime_not_dividing d
1.172 - val g_p = poly_centr(((normirt_polynom ( mod_gcd_p a b p) p) %* (d mod p)) mod_poly p) p
1.173 + val g_p = poly_centr(((normalise ( mod_gcd_p a b p) p) %* (d mod p)) mod_poly p) p
1.174 in if (deg g_p) = 0 then [1]
1.175 else let val g = pp (WHILE a b d M p g_p p)
1.176 in
2.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Wed Jan 23 17:29:49 2013 +0100
2.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Thu Jan 24 17:17:03 2013 +0100
2.3 @@ -16,8 +16,8 @@
2.4
2.5 subsection {* a variant for div *}
2.6
2.7 -fun div2 :: "int => int => int" (infixl "div2" 70)
2.8 - where "a div2 b = (if a div b < 0 then (\<bar>a\<bar> div \<bar>b\<bar>) * -1 else a div b)"
2.9 +fun div2 :: "int => int => int" (infixl "div2" 70) where
2.10 +"a div2 b = (if a div b < 0 then (\<bar>a\<bar> div \<bar>b\<bar>) * -1 else a div b)"
2.11 text {* Meindl: infix div' ... div' didnt work, infix didnt work *}
2.12
2.13 text {*"----------- fun div2 -- --------------------------------"*};
2.14 @@ -34,19 +34,19 @@
2.15
2.16 subsection {* modulo calculations for integers *}
2.17
2.18 -function modi :: "int => int => int => int => int"
2.19 - where "modi r rold m rinv =
2.20 +function modi :: "int => int => int => int => int" where
2.21 + "modi r rold m rinv =
2.22 (if rinv < m
2.23 - then
2.24 - if r mod m = 1
2.25 - then rinv
2.26 - else modi (rold * (rinv + 1)) rold m (rinv + 1)
2.27 - else 0)"
2.28 + then
2.29 + if r mod m = 1
2.30 + then rinv
2.31 + else modi (rold * (rinv + 1)) rold m (rinv + 1)
2.32 + else 0)"
2.33 by auto
2.34 termination sorry
2.35 (* modi is just local for mod_inv *)
2.36 -fun mod_inv :: "int => int => int"
2.37 - where "mod_inv r m = modi r r m 1"
2.38 +fun mod_inv :: "int => int => int" where
2.39 + "mod_inv r m = modi r r m 1"
2.40
2.41 text {*"----------- fun mod_inv --------------------------------"*}
2.42 value "modi 5 5 7 1" -- "= 3"
2.43 @@ -57,19 +57,19 @@
2.44 value "mod_inv 3 7" -- "= 5"
2.45 value "mod_inv 4 339" -- "= 85"
2.46
2.47 -fun mod_div :: "int => int => int => int"
2.48 - where "mod_div a b m = a * (mod_inv b m) mod m"
2.49 +fun mod_div :: "int => int => int => int" where
2.50 + "mod_div a b m = a * (mod_inv b m) mod m"
2.51
2.52 -function root1 :: "int => int => int"
2.53 - where "root1 a n = (if n * n < a then root1 a (n + 1) else n)"
2.54 - by auto
2.55 - termination sorry
2.56 +function root1 :: "int => int => int" where
2.57 + "root1 a n = (if n * n < a then root1 a (n + 1) else n)"
2.58 +by auto
2.59 +termination sorry
2.60 (* root1 is only local to approx_root *)
2.61 -fun approx_root :: "int => int"
2.62 - where "approx_root a = root1 a 1"
2.63 +fun approx_root :: "int => int" where
2.64 + "approx_root a = root1 a 1"
2.65
2.66 -fun chinese_remainder :: "int => int => int => int => int"
2.67 - where "chinese_remainder r1 m1 r2 m2 =
2.68 +fun chinese_remainder :: "int => int => int => int => int" where
2.69 + "chinese_remainder r1 m1 r2 m2 =
2.70 (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1"
2.71
2.72 text {*"----------- fun chinese_remainder ----------------------"*}
2.73 @@ -79,8 +79,8 @@
2.74 subsection {* operations with primes *}
2.75
2.76 (* assumes: 'primes' contains all of first primes /\ n =< (max primes)^2 *)
2.77 -fun is_prime :: "int list \<Rightarrow> int \<Rightarrow> bool"
2.78 - where "is_prime primes n =
2.79 +fun is_prime :: "int list \<Rightarrow> int \<Rightarrow> bool" where
2.80 + "is_prime primes n =
2.81 (if List.length primes > 0
2.82 then
2.83 if (n mod (List.hd primes)) = 0
2.84 @@ -95,24 +95,24 @@
2.85 value "is_prime [2, 3] 6" -- "= False"
2.86 value "is_prime [2, 3] 25" -- "= True ... because 5 not in list"
2.87
2.88 -(* make_prims is local to make_primes *)
2.89 -function make_primes :: "int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list"
2.90 - where "make_primes primes last_p n =
2.91 +(* make_primes is local to primes_upto only *)
2.92 +function make_primes :: "int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list" where
2.93 + "make_primes primes last_p n =
2.94 (if (List.nth primes (List.length primes - 1)) < n
2.95 then
2.96 if is_prime primes (last_p + 2)
2.97 then make_primes (primes @ [(last_p + 2)]) (last_p + 2) n
2.98 else make_primes primes (last_p + 2) n
2.99 else primes)"
2.100 - by pat_completeness (simp del: is_prime.simps)
2.101 - termination sorry
2.102 +by pat_completeness (simp del: is_prime.simps)
2.103 +termination sorry
2.104
2.105 (* make_primes :: int \<Rightarrow> int list
2.106 make_primes n = [2, 3, 5, 7, .., p_k] -- sloppy formulations
2.107 assumes 0 < n
2.108 yields SOME k. ([2, 3, 5, 7, .., p_k] contains all primes \<le> p_n) \<and> n \<le> p_k *)
2.109 -fun primes_upto :: "int \<Rightarrow> int list"
2.110 - where "primes_upto n = (if n < 3 then [2] else make_primes [2,3] 3 n)"
2.111 +fun primes_upto :: "int \<Rightarrow> int list" where
2.112 + "primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"
2.113
2.114 value "primes_upto 2" -- "[2]"
2.115 value "primes_upto 3" -- "[2, 3]"
2.116 @@ -162,34 +162,103 @@
2.117
2.118 subsection {* auxiliary functions for univariate polynomials *}
2.119 type_synonym unipoly = "int list"
2.120 -value "[1,2,3,4,5] :: unipoly"
2.121 +value "[0, 1, 2, 3, 4, 5] :: unipoly"
2.122
2.123 -fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list"
2.124 - where "nth_drop n xs =
2.125 - List.take n xs @ List.drop (n + 1) xs"
2.126 +(* not in List.thy, copy from library.ML *)
2.127 +fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list" where
2.128 + "nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"
2.129
2.130 -function lc :: "unipoly \<Rightarrow> int"
2.131 - where "lc (up :: unipoly) =
2.132 - (if List.nth up (List.length up - 1) \<noteq> 0
2.133 - then List.nth up (List.length up - 1)
2.134 - else lc (nth_drop (List.length up - 1) up))"
2.135 - by auto
2.136 - termination sorry
2.137 +(* leading coefficient *)
2.138 +function lc :: "unipoly \<Rightarrow> int" where
2.139 + "lc uvp =
2.140 + (if List.nth uvp (List.length uvp - 1) \<noteq> 0
2.141 + then List.nth uvp (List.length uvp - 1)
2.142 + else lc (nth_drop (List.length uvp - 1) uvp))"
2.143 +by auto
2.144 +termination sorry
2.145
2.146 +text {*"----------- fun lc -------------------------------------"*}
2.147 value "lc [3, 4, 5, 6]" -- "= 6"
2.148 value "lc [3, 4, 5, 6, 0]" -- "= 6"
2.149
2.150 -function deg :: "unipoly \<Rightarrow> nat"
2.151 - where "deg (up :: unipoly) =
2.152 - (if nth up (length up - 1) \<noteq> 0
2.153 - then length up - 1
2.154 - else deg (nth_drop (length up - 1) up))"
2.155 - by auto
2.156 - termination sorry
2.157 +(* degree *)
2.158 +function deg :: "unipoly \<Rightarrow> nat" where
2.159 + "deg uvp =
2.160 + (if nth uvp (length uvp - 1) \<noteq> 0
2.161 + then length uvp - 1
2.162 + else deg (nth_drop (length uvp - 1) uvp))"
2.163 +by auto
2.164 +termination sorry
2.165 +(* TODO: ?let l = length uvp - 1? more efficient??? compare drop_lc0 !!! *)
2.166
2.167 +text {*"----------- fun deg ------------------------------------"*}
2.168 value "deg [3, 4, 5, 6]" -- "= 3"
2.169 value "deg [3, 4, 5, 6, 0]" -- "= 3"
2.170
2.171 +(* drop leading coefficients equal 0 *)
2.172 +fun drop_lc0 :: "unipoly \<Rightarrow> unipoly" where
2.173 + "drop_lc0 [] = []"
2.174 +| "drop_lc0 uvp =
2.175 + (let l = List.length uvp - 1
2.176 + in
2.177 + if nth uvp l = 0
2.178 + then drop_lc0 (nth_drop l uvp)
2.179 + else uvp)"
2.180 +
2.181 +text {*"----------- fun drop_lc0 -------------------------------";*}
2.182 +value "drop_lc0 [0, 1, 2, 3, 4, 5, 0, 0]" -- "= [0, 1, 2, 3, 4, 5]"
2.183 +value "drop_lc0 [0, 1, 2, 3, 4, 5]" -- "= [0, 1, 2, 3, 4, 5]"
2.184 +
2.185 +(* norm is local to normalise only *)
2.186 +fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> int (*nat?*) \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where
2.187 + "norm pol nrm m lcp i =
2.188 + (if i = 0
2.189 + then [mod_div (nth pol i) lcp m] @ nrm
2.190 + else norm pol ([mod_div (nth pol i) lcp m] @ nrm) m lcp (i - 1))"
2.191 +(* ??? *)
2.192 +fun normalise :: "unipoly \<Rightarrow> int (*nat?*) \<Rightarrow> unipoly" where
2.193 + "normalise p m =
2.194 + (let
2.195 + p = drop_lc0 p;
2.196 + lcp = lc p
2.197 + in
2.198 + if List.length p = 0 then p else norm p [] m lcp (List.length p - 1))"
2.199 +
2.200 +text {*"----------- fun normalise ------------------------------";*}
2.201 +value "normalise [-18, -15, -20, 12, 20, -13, 2] 5" -- "= [1, 0, 0, 1, 0, 1, 1]"
2.202 +
2.203 +(* scalar multiplication *)
2.204 +fun imult :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "imult" 70) where
2.205 + "p imult m = map (op * m) p"
2.206 +
2.207 +text {*"----------- fun %* ------------------------------";*}
2.208 +value "[5, 4, 7, 8, 1] imult 5" -- "= [25, 20, 35, 40, 5]"
2.209 +value "[5, 4, -7, 8, -1] imult 5" -- "= [25, 20, -35, 40, -5]"
2.210 +
2.211 +(* scalar divison *)
2.212 +fun swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where
2.213 + "swapf f a b = f b a"
2.214 +fun idiv :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "idiv" 70) where
2.215 + "p idiv m = map (swapf op div2 m) p"
2.216 +
2.217 +text {*"----------- fun %/ -------------------------------";*}
2.218 +value "[4, 3, 2, 5, 6] idiv 3" -- "= [1, 1, 0, 1, 2]"
2.219 +value "[4, 3, 2, 0] idiv 3" -- "= [1, 1, 0, 0]"
2.220 +
2.221 +(* not in List.thy, copy from library.ML *)
2.222 +fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where
2.223 + "map2 _ [] [] = []"
2.224 +| "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys"
2.225 +| "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
2.226 +
2.227 +(* minus of polys *)
2.228 +fun pminus :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "pminus" 70) where
2.229 + "p1 pminus p2 = map2 (op -) p1 p2"
2.230 +
2.231 +text {*"----------- fun %-% ------------------------";*}
2.232 +value "[8, -7, 0, 1] pminus [-2, 2, 3, 0]" -- "= [10, -9, -3, 1]"
2.233 +value "[8, 7, 6, 5, 4] pminus [2, 2, 3, 1, 1]" -- "= [6, 5, 3, 4, 3]"
2.234 +
2.235
2.236
2.237
3.1 --- a/test/Tools/isac/Knowledge/gcd_poly.sml Wed Jan 23 17:29:49 2013 +0100
3.2 +++ b/test/Tools/isac/Knowledge/gcd_poly.sml Thu Jan 24 17:17:03 2013 +0100
3.3 @@ -17,16 +17,16 @@
3.4 "----------- fun next_prime_not_dividing ----------------";
3.5 "----------- fun lc -------------------------------------";
3.6 "----------- fun deg ------------------------------------";
3.7 -"----------- fun lc_of_unipoly_not_0 --------------------";
3.8 -"----------- fun normirt_polynom ------------------------";
3.9 +"----------- fun drop_lc0 -------------------------------";
3.10 +"----------- fun normalise ------------------------------";
3.11 (* order until here: (GCD_Poly_FP =) GCD_Poly = gcd_poly *)
3.12 "----------- fun %* -------------------------------------";
3.13 "----------- fun %/ -------------------------------------";
3.14 "----------- fun %-% ------------------------------------";
3.15 +"----------- fun %|% ------------------------------------";
3.16 (*"----------- fun %+% ------------------------------------";*)
3.17 +"----------- fun chinese_remainder_poly -----------------";
3.18 "----------- fun mod_poly -------------------------------";
3.19 -"----------- fun chinese_remainder_poly -----------------";
3.20 -"----------- fun %|% ------------------------------------";
3.21 "----------- fun mod_poly_gcd ---------------------------";
3.22 (*"----------- fun poly_centr -----------------------------";*)
3.23 "----------- fun mod_gcd_p ------------------------------";
3.24 @@ -139,51 +139,43 @@
3.25 if deg [3,4,5,6] = 3 then () else error "lc (3,4,5,6) = 3 changed" ;
3.26 if deg [3,4,5,6,0] = 3 then () else error "lc (3,4,5,6) = 6 changed";
3.27
3.28 -"----------- fun lc_of_unipoly_not_0 --------------------";
3.29 -"----------- fun lc_of_unipoly_not_0 --------------------";
3.30 -"----------- fun lc_of_unipoly_not_0 --------------------";
3.31 -if lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] then ()
3.32 -else error "lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] changed";
3.33 -if lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0,1, 2, 3, 4, 5] then ()
3.34 -else error "lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0, 1, 2, 3, 4, 5] changed";
3.35 +"----------- fun drop_lc0 -------------------------------";
3.36 +"----------- fun drop_lc0 -------------------------------";
3.37 +"----------- fun drop_lc0 -------------------------------";
3.38 +if drop_lc0 [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5] then ()
3.39 +else error "drop_lc0 [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5] changed";
3.40 +if drop_lc0 [0, 1, 2, 3, 4, 5] = [0, 1, 2, 3, 4, 5] then ()
3.41 +else error "drop_lc0 [0, 1, 2, 3, 4, 5]=[0, 1, 2, 3, 4, 5] changed";
3.42
3.43 -"----------- fun normirt_polynom ------------------------";
3.44 -"----------- fun normirt_polynom ------------------------";
3.45 -"----------- fun normirt_polynom ------------------------";
3.46 -if (normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] then ()
3.47 -else error "(normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] changed"
3.48 -
3.49 -"----------- fun mod_poly -------------------------------";
3.50 -"----------- fun mod_poly -------------------------------";
3.51 -"----------- fun mod_poly -------------------------------";
3.52 -if ([5,4,7,8,1] mod_poly 5) = [0, 4, 2, 3, 1] then ()
3.53 -else error "[5,4,7,8,1] mod_poly 5 = [0, 4, 2, 3, 1] changed" ;
3.54 -if ([5,4,~7,8,~1] mod_poly 5) = [0, 4, 3, 3, 4] then ()
3.55 -else error "[5,4,~7,8,~1] mod_poly 5 = [0, 4, 3, 3, 4] changed" ;
3.56 +"----------- fun normalise ------------------------------";
3.57 +"----------- fun normalise ------------------------------";
3.58 +"----------- fun normalise ------------------------------";
3.59 +if (normalise [~18, ~15, ~20, 12, 20, ~13, 2] 5) = [1, 0, 0, 1, 0, 1, 1 ] then ()
3.60 +else error "(normalise [~18, ~15, ~20, 12, 20, ~13, 2] 5) = [1, 0, 0, 1, 0, 1, 1 ] changed"
3.61
3.62 "----------- fun %* ------------------------------";
3.63 "----------- fun %* ------------------------------";
3.64 "----------- fun %* ------------------------------";
3.65 -if ([5,4,7,8,1] %* 5) = [25, 20, 35, 40, 5] then ()
3.66 -else error "[5,4,7,8,1] %* 5 = [25, 20, 35, 40, 5] changed";
3.67 -if ([5,4,~7,8,~1] %* 5) = [25, 20, ~35, 40, ~5] then ()
3.68 -else error "[5,4,~7,8,~1] %* 5 = [25, 20, ~35, 40, ~5] changed";
3.69 +if ([5, 4, 7, 8, 1] %* 5) = [25, 20, 35, 40, 5] then ()
3.70 +else error "[5, 4, 7, 8, 1] %* 5 = [25, 20, 35, 40, 5] changed";
3.71 +if ([5, 4, ~7, 8, ~1] %* 5) = [25, 20, ~35, 40, ~5] then ()
3.72 +else error "[5, 4, ~7, 8, ~1] %* 5 = [25, 20, ~35, 40, ~5] changed";
3.73
3.74 "----------- fun %/ -------------------------------";
3.75 "----------- fun %/ -------------------------------";
3.76 "----------- fun %/ -------------------------------";
3.77 -if ([4,3,2,5,6] %/ 3) =[1, 1, 0, 1, 2] then ()
3.78 -else error "%/ [4,3,2,5,6] 3 =[1, 1, 0, 1, 2] changed";
3.79 -if ([4,3,2,0] %/ 3) =[1, 1, 0, 0] then ()
3.80 -else error "%/ [4,3,2,0] 3 =[1, 1, 0, 0] changed";
3.81 +if ([4, 3, 2, 5, 6] %/ 3) = [1, 1, 0, 1, 2] then ()
3.82 +else error "%/ [4, 3, 2, 5, 6] 3 = [1, 1, 0, 1, 2] changed";
3.83 +if ([4, 3, 2, 0] %/ 3) = [1, 1, 0, 0] then ()
3.84 +else error "%/ [4, 3, 2, 0] 3 = [1, 1, 0, 0] changed";
3.85
3.86 -"----------- fun %-% ------------------------";
3.87 -"----------- fun %-% ------------------------";
3.88 -"----------- fun %-% ------------------------";
3.89 -if ([8,~7,0,1] %-% [~2,2,3,0])=[10,~9,~3,1] then ()
3.90 -else error "[8,~7,0,1] %-% [~2,2,3,0]=[10,~9,~3,1] changed";
3.91 -if ([8,7,6,5,4] %-% [2,2,3,1,1])=[6,5,3,4,3] then ()
3.92 -else error "[8,7,6,5,4] %-% [2,2,3,1,1]=[6,5,3,4,3] changed";
3.93 +"----------- fun %-% ------------------------------";
3.94 +"----------- fun %-% ------------------------------";
3.95 +"----------- fun %-% -----------------------------";
3.96 +if ([8, ~7, 0, 1] %-% [~2, 2, 3, 0]) = [10, ~9, ~3, 1] then ()
3.97 +else error "([8, ~7, 0, 1] %-% [~2, 2, 3, 0]) = [10, ~9, ~3, 1] changed";
3.98 +if ([8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1]) = [6, 5, 3, 4, 3] then ()
3.99 +else error "([8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1]) = [6, 5, 3, 4, 3] changed";
3.100
3.101 (*"----------- fun %+% ------------------------------";
3.102 "----------- fun %+% ------------------------------";
3.103 @@ -193,34 +185,42 @@
3.104 if ([8,7,6,5,4] %+% [2,2,3,1,1])=[10,9,9,6,5] then ()
3.105 else error "[8,7,6,5,4] %+% [2,2,3,1,1]=[10,9,9,6,5] changed";*)
3.106
3.107 -"----------- fun chinese_remainder_poly -----------------";
3.108 -"----------- fun chinese_remainder_poly -----------------";
3.109 -"----------- fun chinese_remainder_poly -----------------";
3.110 -if (chinese_remainder_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()
3.111 -else error "chinese_remainder_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";
3.112 -
3.113 "----------- fun %|% -------------------------------";
3.114 "----------- fun %|% -------------------------------";
3.115 "----------- fun %|% -------------------------------";
3.116 "~~~~~ fun %|% , args:"; val (poly1, poly2) = ([9,12,4],[3,2]);
3.117 val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
3.118 val c = lc poly1 div lc b;
3.119 -val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
3.120 +val rest = drop_lc0 ( poly1 %-% (b %* c ));
3.121 rest = [];(*=false*)
3.122 c<>0 andalso length rest >= length poly2; (*=true*)
3.123 "~~~~~ fun %|% , args:"; val (poly1, poly2) = (rest, poly2);
3.124 val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
3.125 val c = lc poly1 div lc b;
3.126 -val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
3.127 +val rest = drop_lc0 ( poly1 %-% (b %* c ));
3.128 rest = [];(*=true*)
3.129 "~~~~~ to return val:"; val (divide) = (true);
3.130
3.131 if ([4] %|% [6]) = false then () else error "[4] %|% [6] = false changed";
3.132 -if ( [8] %|%[16,0]) = true then () else error "[8] %|%[16,0] = true changed";
3.133 +if ( [8] %|%[16,0]) = true then () else error "[8] %|% [16,0] = true changed";
3.134 if ([3,2] %|%[0,0,9,12,4] ) = true then ()
3.135 else error "[3,2] %|%[0,0,9,12,4] = true changed";
3.136 if ([8,0] %|% [16]) = true then () else error "[8,0] %|% [16] = true changed";
3.137
3.138 +"----------- fun chinese_remainder_poly -----------------";
3.139 +"----------- fun chinese_remainder_poly -----------------";
3.140 +"----------- fun chinese_remainder_poly -----------------";
3.141 +if (chinese_remainder_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()
3.142 +else error "chinese_remainder_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";
3.143 +
3.144 +"----------- fun mod_poly -------------------------------";
3.145 +"----------- fun mod_poly -------------------------------";
3.146 +"----------- fun mod_poly -------------------------------";
3.147 +if ([5,4,7,8,1] mod_poly 5) = [0, 4, 2, 3, 1] then ()
3.148 +else error "[5,4,7,8,1] mod_poly 5 = [0, 4, 2, 3, 1] changed" ;
3.149 +if ([5,4,~7,8,~1] mod_poly 5) = [0, 4, 3, 3, 4] then ()
3.150 +else error "[5,4,~7,8,~1] mod_poly 5 = [0, 4, 3, 3, 4] changed" ;
3.151 +
3.152 "----------- fun mod_poly_gcd ------------------------";
3.153 "----------- fun mod_poly_gcd ------------------------";
3.154 "----------- fun mod_poly_gcd ------------------------";
3.155 @@ -235,33 +235,33 @@
3.156 "~~~~~ fun mod_poly_gcd , args:";
3.157 val (a,b,m) = ([ ~13, 2,12],[ ~14, 2],5);
3.158 val moda = a mod_poly m;
3.159 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b mod_poly m) ;
3.160 +val modb = (replicate (length a - length(drop_lc0 b)) 0) @ (drop_lc0 b mod_poly m) ;
3.161 val c = mod_div (lc moda) (lc modb) m;
3.162 -val rest = lc_of_unipoly_not_0 (moda %-% (modb %* c) mod_poly m) mod_poly m;
3.163 +val rest = drop_lc0 (moda %-% (modb %* c) mod_poly m) mod_poly m;
3.164 rest = []; (*=false*)
3.165 length rest < length b; (*=false*)
3.166 "~~~~~ fun mod_poly_gcd , args:";
3.167 val (a,b,m,d) = (rest, b , m , [c]);
3.168 val moda = a mod_poly m
3.169 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b mod_poly m) ;
3.170 +val modb = (replicate (length a - length(drop_lc0 b)) 0) @ (drop_lc0 b mod_poly m) ;
3.171 val c = mod_div (lc moda) (lc modb) m
3.172 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) mod_poly m);
3.173 +val rest = drop_lc0 ((moda %-% (modb %* c)) mod_poly m);
3.174 rest = [];(*=flase*)
3.175 length rest < length b; (*=true*)
3.176 "~~~~~ fun mod_poly_gcd , args:";
3.177 val (a,b,m,d) = (b, rest, m, [c] @ d);
3.178 val moda = a mod_poly m
3.179 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b mod_poly m) ;
3.180 +val modb = (replicate (length a - length(drop_lc0 b)) 0) @ (drop_lc0 b mod_poly m) ;
3.181 val c = mod_div (lc moda) (lc modb) m
3.182 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) mod_poly m);
3.183 +val rest = drop_lc0 ((moda %-% (modb %* c)) mod_poly m);
3.184 rest = [];(*=flase*)
3.185 length rest < length b; (*=false*)
3.186 "~~~~~ fun mod_poly_gcd , args:";
3.187 val (a,b,m,d) = (b, rest, m, [c] @ d);
3.188 val moda = a mod_poly m
3.189 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b mod_poly m) ;
3.190 +val modb = (replicate (length a - length(drop_lc0 b)) 0) @ (drop_lc0 b mod_poly m) ;
3.191 val c = mod_div (lc moda) (lc modb) m
3.192 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) mod_poly m);
3.193 +val rest = drop_lc0 ((moda %-% (modb %* c)) mod_poly m);
3.194 rest = [];(*=true*)
3.195 "~~~~~ to return val:"; val ([b, rest, lsit])=[b, [0], [c] @ d];
3.196
4.1 --- a/test/Tools/isac/Test_Isac.thy Wed Jan 23 17:29:49 2013 +0100
4.2 +++ b/test/Tools/isac/Test_Isac.thy Thu Jan 24 17:17:03 2013 +0100
4.3 @@ -155,6 +155,7 @@
4.4 use "Knowledge/atools.sml"
4.5 use "Knowledge/simplify.sml"
4.6 use "Knowledge/poly.sml"
4.7 + use "Knowledge/gcd_poly_winkler.sml"
4.8 (*use "Knowledge/rational.sml" WN120317.TODO postponed to joint work with dmeindl *)
4.9 use "Knowledge/equation.sml"
4.10 use "Knowledge/root.sml"
5.1 --- a/test/Tools/isac/Test_Some2.thy Wed Jan 23 17:29:49 2013 +0100
5.2 +++ b/test/Tools/isac/Test_Some2.thy Thu Jan 24 17:17:03 2013 +0100
5.3 @@ -6,79 +6,72 @@
5.4 use"~~/test/Tools/isac/Knowledge/gcd_poly.sml"
5.5
5.6
5.7 -ML {* " ~~~ 1 ~~~";
5.8 -val thy = @{theory GCD_Poly};
5.9 -val ctxt = Proof_Context.init_global thy;
5.10 +ML {*
5.11 print_depth 7;
5.12 *}
5.13 -ML {*
5.14 -"----------- fun normirt_polynom --------------------";
5.15 -val a = [~18, ~15, ~20, 12, 20, ~13, 2];
5.16 -normirt_polynom a 3;
5.17 -val b = [8, 28, 22, ~11, ~14, 1, 2];
5.18 -val c = [0, 0, 3, 0, 0];
5.19 -*}
5.20
5.21 -
5.22 -ML {*
5.23 -PolyML.makestring
5.24 -*}
5.25 -
5.26 -ML {*
5.27 -fun GCD_MOD a =
5.28 - (writeln (PolyML.makestring a); GOTO2 a - 1)
5.29 -and GOTO2 b =
5.30 - if b > 0 then GCD_MOD (b - 1) else b
5.31 -*}
5.32 -
5.33 -ML {*
5.34 -GCD_MOD 5
5.35 -*}
5.36 -
5.37 -ML{*
5.38 -
5.39 -*}
5.40 -
5.41 - ML{*
5.42 +ML{*
5.43 if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5) = [1, 1, 1, 1] then ()
5.44 else error "mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5 = [1, 1, 1, 1] changed";
5.45 if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7) = [5, 1, 0, 5, 1] then ()
5.46 else error "mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [5, 1, 0, 5, 1] changed";
5.47 -
5.48 *}
5.49
5.50 -ML {*
5.51 +ML {*
5.52 + fun mod_pol uvp m n =
5.53 + (if n < length uvp
5.54 + then mod_pol ((List.tl uvp) @ [(List.hd uvp) mod m]) m (n + 1 )
5.55 + else uvp)
5.56 +
5.57 + infix mod_poly
5.58 + fun (uvp: unipoly) mod_poly m = mod_pol uvp m 0;
5.59 +
5.60 +"----------- fun mod_poly -------------------------------";
5.61 +if ([5,4,7,8,1] mod_poly 5) = [0, 4, 2, 3, 1] then ()
5.62 +else error "[5,4,7,8,1] mod_poly 5 = [0, 4, 2, 3, 1] changed" ;
5.63 +"======================================================================= later"
5.64 *}
5.65
5.66
5.67 -ML{*
5.68 -*}
5.69 -
5.70 ML {*
5.71 -*}
5.72 -
5.73 -ML {*
5.74 -*}
5.75 -
5.76 -ML {*
5.77 -*}
5.78 -
5.79 -ML {*
5.80 +*}
5.81 +ML {*
5.82 + infix %|%
5.83 + fun [b: int] %|% [a: int] =
5.84 + if abs b<= abs a andalso a mod b = 0
5.85 + then true
5.86 + else false
5.87 + | (poly2: unipoly) %|% (poly1: unipoly) =
5.88 + let
5.89 + val b = (replicate (length poly1 - length(drop_lc0 poly2)) 0) @
5.90 + drop_lc0 poly2 ;
5.91 + val c = lc poly1 div2 lc b;
5.92 + val rest = drop_lc0 (poly1 %-% (b %* c));
5.93 + in
5.94 + if rest = []
5.95 + then true
5.96 + else
5.97 + if c<>0 andalso length rest >= length poly2
5.98 + then poly2 %|% rest
5.99 + else false
5.100 + end
5.101 +*}
5.102 +ML {*
5.103 +"----------- fun %|% -------------------------------";
5.104 +if ([4] %|% [6]) = false then () else error "[4] %|% [6] = false changed";
5.105 +if ([8] %|% [16,0]) = true then () else error "[8] %|% [16,0] = true changed";
5.106 +if ([8,0] %|% [16]) = true then () else error "[8,0] %|% [16] = true changed";
5.107 +if ([3, 2] %|% [0, 0 ,9, 12, 4]) = true then ()
5.108 + else error "[3, 2] %|% [0, 0 ,9, 12, 4] = true changed";
5.109 +[ 2, 3] %|% [ 4, 12, 9 ];
5.110 +[0, 2, 3] %|% [0, 4, 12, 9, 0, 0];
5.111 +[0, 2, 3] %|% [0, 4, 12, 9 ];
5.112 +*}
5.113 +ML {*
5.114 *}
5.115
5.116 -ML{*
5.117 -*}
5.118
5.119 - ML{*
5.120 -*}
5.121 -
5.122 -ML{*
5.123 -*}
5.124 -
5.125 -ML{*
5.126 -*}
5.127 -
5.128 -ML{* (*==================*)
5.129 +ML {* (*==================*)
5.130 *} ML{* "~~~~~ fun , args:"; val () = ();
5.131 "~~~~~ to return val:"; val () = ();
5.132