1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri Aug 06 12:09:06 2021 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri Aug 06 18:27:05 2021 +0200
1.3 @@ -5,11 +5,11 @@
1.4 begin
1.5
1.6 text \<open>
1.7 - This code has been translated from GCD_Poly.thy by Diana Meindl,
1.8 + This code has been translated from GCD_Poly_OLD.thy by Diana Meindl,
1.9 who follows Franz Winkler, Polynomial algorithyms in computer algebra, Springer 1996.
1.10 - Winkler's original identifiers are in test/./gcd_poly_winkler.sml;
1.11 - test/../gcd_poly.sml documents the changes towards GCD_Poly.thy;
1.12 - the style of GCD_Poly.thy has been adapted to the function package.
1.13 + Winkler's original identifiers are in test/./gcd_poly_old.sml;
1.14 + test/../gcd_poly_ml.sml documents the changes towards GCD_Poly_ML.thy;
1.15 + the style of GCD_Poly_OLD.thy has been adapted to the function package.
1.16 \<close>
1.17
1.18 section \<open>Isabelle's predefined polynomials\<close>
2.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_ML.thy Fri Aug 06 12:09:06 2021 +0200
2.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_ML.thy Fri Aug 06 18:27:05 2021 +0200
2.3 @@ -6,7 +6,8 @@
2.4 10 20 30 40 50 60 70 80 90 100
2.5 *)
2.6
2.7 -theory GCD_Poly_ML imports Base_Tools begin
2.8 +theory GCD_Poly_ML imports Base_Tools (** )GCD_Poly_OLD ( *huge loading time*)
2.9 +begin
2.10
2.11 text \<open>
2.12 Below we reference
3.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_OLD.thy Fri Aug 06 12:09:06 2021 +0200
3.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_OLD.thy Fri Aug 06 18:27:05 2021 +0200
3.3 @@ -6,6 +6,8 @@
3.4 (*"~~/src/HOL/Number_Theory/Primes" ??? rm ???, thys for code_gen missing*)
3.5 (*"~~/src/HOL/Library/Code_Target_Numeral" ...note (4)*)
3.6 begin
3.7 +ML_file \<open>gcd_poly_old.sml\<close> (* this ML code was derived into the Isar code below *)
3.8 +
3.9 value "xxx" \<comment> \<open>TODO\<close>
3.10
3.11 subsection \<open>Euclidean Algorithm with generalised division\<close>
3.12 @@ -31,7 +33,7 @@
3.13 Term: \<lambda>a p n. coeff p n div a
3.14 Type: ?'b \<Rightarrow> ?'b poly \<Rightarrow> nat \<Rightarrow> ?'b
3.15
3.16 -Reason: The type of the term cannot be instancied to "'a \<Rightarrow> 'a poly \<Rightarrow> nat \<Rightarrow> 'a".
3.17 +Reason: The type of the term cannot be instantiated to "'a \<Rightarrow> 'a poly \<Rightarrow> nat \<Rightarrow> 'a".
3.18 *)
3.19 (*/------------------------------------------------------------------------------\
3.20 this is an odd bypass for the above problem:
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/src/Tools/isac/Knowledge/gcd_poly_old.sml Fri Aug 06 18:27:05 2021 +0200
4.3 @@ -0,0 +1,1579 @@
4.4 +(* Title: src/../gcd_poly_old.sml
4.5 + Author: Diana Meindl
4.6 + Copyright (c) Diana Meindl 2011
4.7 +
4.8 +Programcode according to [1] for later lookup for mathematicians.
4.9 +The tests below remain according to [1],
4.10 +while "$ISABELLE_ISAC_TEST/Tools/isac/Knowledge/gcd_poly_ml.sml" start exactly from the same state
4.11 +and in time follows development in "$ISABELLE_ISAC/Knowledge/GCD_Poly_ML.thy".
4.12 +
4.13 +This file includes *** tests ***.
4.14 +
4.15 +[1] Franz Winkler. Polynomial Algorithms in Computer Algebra.
4.16 +Springer-Verlag/Wien 1996.
4.17 +*)
4.18 +
4.19 +(*fun nth _ [] = error "nth _ []" (*Isabelle2002, still saved the hours of update*)
4.20 + | nth 1 (x::_) = x
4.21 + | nth n (_::xs) = nth (n- 1) xs;*)
4.22 + fun nth xs i = List.nth (xs, i);
4.23 +"--------------------------------------------------------";
4.24 +"table of contents --------------------------------------";
4.25 +"--------------------------------------------------------";
4.26 +"----------- auxiliary functions univariate -------------";
4.27 +"=========== code for [1] p.93: univariate ==============";
4.28 +"----------- auxiliary functions multivariate -----------";
4.29 +"=========== code for [1] p.95: multivariate ============";
4.30 +"--------------------------------------------------------";
4.31 +"--- begin of univariate polynomials --------------------";
4.32 +"----------- fun mod_inv --------------------------------";
4.33 +"----------- fun CRA_2 ----------------------------------";
4.34 +"----------- fun lc -------------------------------------";
4.35 +"----------- fun deg ------------------------------------";
4.36 +"----------- fun primeList ------------------------------";
4.37 +"----------- fun find_next_prime_not_divide -------------";
4.38 +"----------- fun poly_mod -------------------------------";
4.39 +"----------- fun %* -------------------------------------";
4.40 +"----------- fun %/ -------------------------------------";
4.41 +"----------- fun %-% ------------------------------------";
4.42 +"----------- fun %+% ------------------------------------";
4.43 +"----------- fun CRA_2_poly -----------------------------";
4.44 +"----------- fun mod_div --------------------------------";
4.45 +"----------- fun lc_of_unipoly_not_0 --------------------";
4.46 +"----------- fun %|% ------------------------------------";
4.47 +"----------- fun mod_poly_gcd ---------------------------";
4.48 +"----------- fun normirt_polynom ------------------------";
4.49 +"----------- fun poly_centr -----------------------------";
4.50 +"----------- fun mod_gcd_p ------------------------------";
4.51 +"----------- fun gcd_n ----------------------------------";
4.52 +"----------- fun pp -------------------------------------";
4.53 +"----------- fun GCD_MOD---------------------------------";
4.54 +"--- begin of multivariate polynomials ------------------";
4.55 +"----------- fun get_coef --------------------------- ---";
4.56 +"----------- fun get_varlist ----------------------------";
4.57 +"----------- fun get_coeflist ---------------------------";
4.58 +"----------- fun add_var_to_multipoly -------------------";
4.59 +"----------- fun cero_multipoly -------------------------";
4.60 +"----------- fun short_mv -------------------------------";
4.61 +"----------- fun greater_var ----------------------------";
4.62 +"----------- fun order_multipoly ------------------------";
4.63 +"----------- fun lc_multipoly ---------------------------";
4.64 +"----------- fun deg_multipoly --------------------------";
4.65 +"----------- fun max_deg_var ----------------------------";
4.66 +"----------- infix %%/ ----------------------------------";
4.67 +"----------- fun cont -----------------------------------";
4.68 +"----------- fun cont_multipoly -------------------------";
4.69 +"----------- fun evaluate_mv ----------------------------";
4.70 +"----------- fun mutli_to_uni ---------------------------";
4.71 +"----------- fun find_new_r ----------------------------";
4.72 +"----------- fun newton_mv -----------------------------";
4.73 +"----------- fun mult_with_new_var ---------------------";
4.74 +"----------- fun greater_deg_multipoly -----------------";
4.75 +"----------- infix %%+%% -------------------------------";
4.76 +"----------- infix %%-%% -------------------------------";
4.77 +"----------- infix %%*%% -------------------------------";
4.78 +"----------- fun newton_first --------------------------";
4.79 +"----------- fun newton_mv -----------------------------";
4.80 +"----------- fun listgreater ---------------------------";
4.81 +"----------- finfix %%|%% -----------------------------";
4.82 +"----------- fun GCD_MODm ------------------------------";
4.83 +
4.84 +(*text {*
4.85 + Below we reference
4.86 + F. Winkler, Polynomial Algorithms. ...
4.87 + by page 93 and 95. The identifiers used in this book are re-used in this thesis
4.88 + in order to support reference (although some of these identifiers
4.89 + doe not conform with the Isabelle coding standards)
4.90 +*}*)
4.91 +
4.92 + (*default_print_depth 3; 20*)
4.93 + type unipoly = int list;
4.94 +
4.95 +"----------- auxiliary functions univariate -------------";
4.96 +"----------- auxiliary functions univariate -------------";
4.97 +"----------- auxiliary functions univariate -------------";
4.98 +
4.99 +(*subsection {* calculations for integers *}*)
4.100 +
4.101 + infix div'
4.102 + fun a div' b =
4.103 + if a < 0 then abs a div b * ~1
4.104 + else a div b;
4.105 +
4.106 +(*subsection {* modulo calculations for integers *}*)
4.107 + fun mod_inv r m =
4.108 + let
4.109 + fun modi (r, rold, m, rinv) =
4.110 + if rinv < m
4.111 + then
4.112 + if r mod m = 1
4.113 + then rinv
4.114 + else modi (rold * (rinv + 1), rold, m, rinv + 1)
4.115 + else 0
4.116 + in modi (r, r, m, 1) end;
4.117 +
4.118 + fun mod_div a b m =
4.119 + a * (mod_inv b m) mod m;
4.120 +
4.121 + fun aprox_root a =
4.122 + let fun root' a n =
4.123 + if n*n < a then root' a (n+1) else n
4.124 + in root' a 1
4.125 + end
4.126 +
4.127 + (* r = r1 mod m1 and r= r2 mod m2 *)
4.128 + fun CRA_2 (r1: int, m1: int, r2: int, m2: int ) =
4.129 + (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1
4.130 +
4.131 +(*subsection {* prime opertations *}*)
4.132 +
4.133 + fun is_prime primelist number =
4.134 + if length primelist >0
4.135 + then
4.136 + if (number mod (nth primelist 0))=0
4.137 + then false
4.138 + else is_prime (nth_drop 0 primelist) number
4.139 + else true
4.140 +
4.141 + fun primeList number =
4.142 + let
4.143 + fun make_primelist list last number =
4.144 + if (nth list (length list - 1)) < number
4.145 + then
4.146 + if ( is_prime list (last + 2))
4.147 + then make_primelist (list @ [(last + 2)]) (last + 2) number
4.148 + else make_primelist list (last + 2) number
4.149 + else list
4.150 + in
4.151 + if number < 3
4.152 + then [2]
4.153 + else make_primelist [2,3] 3 number end
4.154 +
4.155 + (*find a prime greater p not dividing the number a*)
4.156 + fun find_next_prime_not_divide a p =
4.157 + let
4.158 + val next = nth (primeList (p + 1)) (length (primeList (p + 1)) - 1) ;
4.159 + in
4.160 + if a mod next <> 0
4.161 + then next
4.162 + else find_next_prime_not_divide a next
4.163 + end
4.164 +
4.165 +(*subsection {* calculations for univariate polynomials *}*)
4.166 +
4.167 + fun lc (uvp: unipoly) =
4.168 + if nth uvp (length uvp- 1) <>0
4.169 + then nth uvp (length uvp- 1)
4.170 + else lc (nth_drop (length uvp- 1) uvp);
4.171 +
4.172 + fun deg (uvp: unipoly) =
4.173 + if nth uvp (length uvp- 1) <>0
4.174 + then length uvp - 1
4.175 + else deg (nth_drop (length uvp- 1) uvp)
4.176 +
4.177 + fun lc_of_unipoly_not_0 [] = [](* and delete lc=0*)
4.178 + | lc_of_unipoly_not_0 (uvp: unipoly) =
4.179 + if nth uvp (length uvp - 1) =0
4.180 + then lc_of_unipoly_not_0 (nth_drop (length uvp- 1) uvp)
4.181 + else uvp;
4.182 +
4.183 + fun normirt_polynom (poly1: unipoly) (m: int) =
4.184 + let
4.185 + val poly1 = lc_of_unipoly_not_0 poly1
4.186 + val lc_a=lc poly1;
4.187 + fun normirt poly1 b m lc_a i =
4.188 + if i=0
4.189 + then [mod_div (nth poly1 i) lc_a m]@b
4.190 + else normirt poly1 ( [mod_div(nth poly1 i) lc_a m]@b) m lc_a (i- 1) ;
4.191 + in
4.192 + if length(poly1)=0
4.193 + then poly1
4.194 + else normirt poly1 [] m lc_a (length(poly1)- 1)
4.195 + end
4.196 +
4.197 + infix %*
4.198 + fun (a: unipoly) %* b = map2 Integer.mult a (replicate (length (a)) b)
4.199 +
4.200 + infix %/
4.201 + fun (poly: unipoly) %/ b = (* =quotient*)
4.202 + let fun division poly b = poly div' b;
4.203 + in map2 division poly (replicate (length (poly)) b) end
4.204 +
4.205 + infix %-
4.206 + fun (poly: unipoly) %- b =
4.207 + let fun minus poly b = poly - b;
4.208 + in map2 minus poly (replicate (length (poly)) b) end
4.209 +
4.210 + infix %+%
4.211 + fun (a: unipoly) %+% (b: unipoly) = map2 Integer.add a b
4.212 +
4.213 + infix %-%
4.214 + fun (a: unipoly) %-% (b: unipoly) =
4.215 + let fun minus a b = a-b;
4.216 + in map2 minus a b end
4.217 +
4.218 + (* if poly2|poly1 *)
4.219 + infix %|%
4.220 + fun [b: int] %|% [a: int] =
4.221 + if abs b<= abs a andalso a mod b = 0
4.222 + then true
4.223 + else false
4.224 + | (poly2: unipoly) %|% (poly1: unipoly) =
4.225 + let
4.226 + val b = (replicate (length poly1 - length(lc_of_unipoly_not_0 poly2)) 0) @ lc_of_unipoly_not_0 poly2 ;
4.227 + val c = lc poly1 div' lc b;
4.228 + val rest = lc_of_unipoly_not_0 (poly1 %-% (b %* c));
4.229 + in
4.230 + if rest = []
4.231 + then true
4.232 + else
4.233 + if c<>0 andalso length rest >= length poly2
4.234 + then poly2 %|% rest
4.235 + else false
4.236 + end
4.237 +
4.238 + fun poly_centr (poly: unipoly) (m: int) =
4.239 + let
4.240 + fun midle m =
4.241 + let val mid = m div' 2
4.242 + in
4.243 + if m mod mid = 0
4.244 + then mid
4.245 + else mid+1
4.246 + end
4.247 + fun centr a m mid =
4.248 + if a > mid
4.249 + then a - m
4.250 + else a
4.251 + fun polyCentr poly poly' m mid counter =
4.252 + if length(poly) > counter
4.253 + then polyCentr poly (poly' @ [centr (nth poly counter) m mid]) m mid (counter+1)
4.254 + else (poly': unipoly)
4.255 + in polyCentr poly [] m (midle m) 0
4.256 + end
4.257 +
4.258 + fun sum_lmb (a: unipoly) exp =
4.259 + let fun sum' [a] _ = a
4.260 + | sum' (a: int list) exp =
4.261 + sum' ([((nth a 0)) + ( Integer.pow exp (nth a 1))] @ nth_drop 0 (nth_drop 0 a)) exp
4.262 + in sum' ([Integer.pow exp (nth a 0)] @ nth_drop 0 a) exp
4.263 + end;
4.264 +Integer.min ;
4.265 + fun landau_mignotte_bound a b =
4.266 + (Integer.pow (Integer.min (deg a) (deg b)) 2) * (abs (Integer.gcd (lc a) (lc b)))*
4.267 + (Int.min (abs((aprox_root (sum_lmb a 2)) div ~(lc a)), abs(( aprox_root (sum_lmb b 2)) div ~(lc b))));
4.268 +
4.269 +(*subsection {* modulo calculations for polynomials *}*)
4.270 +
4.271 + fun CRA_2_poly (ma, mb) (r1, r2) =
4.272 + let
4.273 + fun CRA_2' (ma, mb) (r1, r2) = CRA_2 (r1, ma,r2, mb)
4.274 + in map (CRA_2' (ma, mb)) (r1 ~~ r2) end
4.275 +
4.276 + infix poly_mod
4.277 + fun uvp poly_mod m =
4.278 + let fun poly_mod' uvp m n =
4.279 + if n < length (uvp)
4.280 + then poly_mod' ((nth_drop 0 uvp)@[(nth uvp 0) mod m]) m ( n + 1 )
4.281 + else uvp (*end of poly_mod'*)
4.282 + in
4.283 + poly_mod' uvp m 0 end
4.284 +
4.285 + fun mod_poly_gcd (upoly1: unipoly) (upoly2: unipoly) (m: int) =
4.286 + let
4.287 + val moda = upoly1 poly_mod m;
4.288 + val modb = (replicate (length upoly1 - length(lc_of_unipoly_not_0 (upoly2 poly_mod m))) 0)
4.289 + @ (lc_of_unipoly_not_0 (upoly2 poly_mod m)) ;
4.290 + val c = mod_div (lc moda) (lc modb) m
4.291 + val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m)
4.292 + in if rest = []
4.293 + then [upoly2, [0]]
4.294 + else
4.295 + if length rest < length upoly2
4.296 + then mod_poly_gcd upoly2 rest m
4.297 + else mod_poly_gcd rest upoly2 m
4.298 + end
4.299 +
4.300 + fun mod_gcd_p (poly1: unipoly) (poly2: unipoly) (p: int) =
4.301 + let val gcd_p = mod_poly_gcd poly1 poly2 p
4.302 + in normirt_polynom (nth gcd_p 0) p
4.303 + end
4.304 +
4.305 + fun gcd_n (up: unipoly) =
4.306 + if length up = 2
4.307 + then abs (Integer.gcd (nth up 0)(nth up 1))
4.308 + else gcd_n ([abs (Integer.gcd (nth up 0)(nth up 1))]@(nth_drop 0 (nth_drop 0 up)))
4.309 +
4.310 +fun gcd_n up = abs (Integer.gcds up);
4.311 +
4.312 + fun pp [_: int] = [1]
4.313 + | pp (poly1: unipoly) =
4.314 + if (poly1 poly_mod (gcd_n poly1)) = (replicate (length (poly1)) 0)
4.315 + then poly1 %/ (gcd_n poly1)
4.316 + else poly1;
4.317 +
4.318 +"=========== code for [1] p.93: univariate ==============";
4.319 +"=========== code for [1] p.93: univariate ==============";
4.320 +"=========== code for [1] p.93: univariate ==============";
4.321 +(*subsection {* GCD_MOD Algorgithmus, code for [1] p.93: univariate *}*)
4.322 +
4.323 + fun GCD_MOD (a: unipoly) (b: unipoly) =
4.324 + if a = [0] orelse b = [0] then [1] else
4.325 + if a = b then a else
4.326 + let
4.327 +(*1*)val d = abs (Integer.gcd (lc a) (lc b));
4.328 + val M = 2*d*landau_mignotte_bound a b;
4.329 + fun GOTO2 a b d M p = (*==============================*)
4.330 + let
4.331 +(*2*) val p = find_next_prime_not_divide d p
4.332 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
4.333 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
4.334 + fun GOTO3 a b d M p g_p = (*~~~~~~~~~~~~~~~~~~~~~~*)
4.335 +(*3*) if (deg g_p) = 0
4.336 + then [1]
4.337 + else
4.338 + let
4.339 + val P = p
4.340 + val g = g_p
4.341 + fun WHILE a b d M P g p = (*------------------*)
4.342 +(*4*) if P > M
4.343 + then g
4.344 + else
4.345 + let
4.346 + val p = find_next_prime_not_divide d p
4.347 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
4.348 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
4.349 + in
4.350 + if deg g_p < deg g
4.351 + then GOTO3 a b d M p g_p
4.352 + else
4.353 + if deg g_p = deg g
4.354 + then
4.355 + let
4.356 + val g = CRA_2_poly (P,p)(g,g_p)
4.357 + val P = P*p
4.358 + val g = poly_centr(g poly_mod P)P
4.359 + in WHILE a b d M P g p end
4.360 + else WHILE a b d M P g p
4.361 + end (*----------------------------------*)
4.362 +
4.363 +
4.364 + val g = WHILE a b d M P g p (* << 1.Mal -----*)
4.365 + in g end (*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
4.366 +
4.367 + val g = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)
4.368 +(*5*) val g = pp g
4.369 + in
4.370 + if (g %|% a) andalso (g %|% b)
4.371 + then g
4.372 + else GOTO2 a b d M p
4.373 + end (*==============================================*)
4.374 +
4.375 +
4.376 + val g = GOTO2 a b d M (*p*)1(* << 1. Mal =============*)
4.377 + in g end;
4.378 +
4.379 +"----------- auxiliary functions multivariate -----------";
4.380 +"----------- auxiliary functions multivariate -----------";
4.381 +"----------- auxiliary functions multivariate -----------";
4.382 +(*subsection {* GCD_MODm Algorgithmus, auxiliary functions multivariate *}*)
4.383 +
4.384 +type monom = (int * int list);
4.385 +type multipoly = monom list;
4.386 +
4.387 +(*subsection {* calculations for multivariate polynomials *}*)
4.388 +
4.389 + fun listgreater a b =
4.390 + let fun gr a b = a>=b;
4.391 + fun all [] = true
4.392 + | all list =
4.393 + if nth list 0 then all (nth_drop 0 list) else false
4.394 + in all (map2 gr a b)
4.395 + end
4.396 +
4.397 + fun get_coef (mvp: multipoly) (place: int) =
4.398 + let fun coef ((coef,_): monom) = coef;
4.399 + in coef (nth mvp place) end
4.400 +
4.401 + fun get_varlist (mvp: multipoly) (place: int) =
4.402 + let fun varlist ((_,var): monom) = var;
4.403 + in varlist (nth mvp place) end
4.404 +
4.405 + fun get_coeflist (poly: multipoly) =
4.406 + let
4.407 + fun get_coeflist' (poly: multipoly) list =
4.408 + if poly = []
4.409 + then
4.410 + list
4.411 + else
4.412 + get_coeflist' (nth_drop 0 poly) (list @ [get_coef poly 0])
4.413 + in
4.414 + get_coeflist' poly []
4.415 + end
4.416 +
4.417 + fun add_var_to_multipoly (mvp: multipoly) (order: int) =
4.418 + let fun add ([]: multipoly) (_: int) (new: multipoly) = new
4.419 + | add (mvp: multipoly) (order: int) (new: multipoly) =
4.420 + let val (first, last) = chop order (get_varlist mvp 0)
4.421 + in add (nth_drop 0 mvp) (order) (new @ [((get_coef mvp 0),(first @ [0] @ last))]) end
4.422 + in add mvp order [] end
4.423 +
4.424 + fun cero_multipoly 1 = [(0,[0])]
4.425 + | cero_multipoly (diferent_var: int) =
4.426 + add_var_to_multipoly (cero_multipoly (diferent_var- 1)) 0;
4.427 +
4.428 + fun short_mv (mvp: multipoly) =
4.429 + let fun short (mvp: multipoly) (new: multipoly) =
4.430 + if length mvp =1
4.431 + then if (get_coef mvp 0) = 0
4.432 + then if new = [] then cero_multipoly (length (get_varlist mvp 0)) else new
4.433 + else new @ mvp
4.434 + else if get_varlist mvp 0 = get_varlist mvp 1
4.435 + then short ( [(get_coef mvp 0 + get_coef mvp 1,get_varlist mvp 0)]
4.436 + @ nth_drop 0 (nth_drop 0 mvp)) new
4.437 + else if (get_coef mvp 0) = 0 then short (nth_drop 0 mvp) new
4.438 + else short (nth_drop 0 mvp) (new @ [nth mvp 0])
4.439 + in short mvp [] end
4.440 +
4.441 + (* if a is greater than b *)
4.442 + fun greater_var [a: int] [b: int] =
4.443 + a > b
4.444 + | greater_var (a: int list) (b: int list) =
4.445 + if (nth a (length a - 1))= (nth b (length b - 1))
4.446 + then greater_var (nth_drop (length a- 1) a) (nth_drop (length b- 1) b)
4.447 + else (nth a (length a - 1)) > (nth b (length b - 1))
4.448 +
4.449 + fun order_multipoly (a: multipoly)=
4.450 + let
4.451 + val ordered = [nth a 0];
4.452 + fun order_mp [] [] = cero_multipoly (length (get_varlist a 0))
4.453 + | order_mp [] (ordered: multipoly) = short_mv ordered
4.454 + | order_mp a (ordered: multipoly) =
4.455 + if greater_var (get_varlist a 0) (get_varlist ordered (length ordered - 1))
4.456 + then order_mp (nth_drop 0 a)(ordered @ [nth a 0])
4.457 + else let
4.458 + val rest = [nth ordered (length ordered - 1)];
4.459 + fun order_mp' [] (new: multipoly) (rest: multipoly) = new @ rest
4.460 + | order_mp' (ordered': multipoly) (new: multipoly) (rest: multipoly) =
4.461 + if greater_var (get_varlist new 0) (get_varlist ordered' (length ordered' - 1))
4.462 + then ordered' @ new @ rest
4.463 + else order_mp' (nth_drop (length ordered' - 1) ordered') new
4.464 + ([nth ordered' (length ordered' - 1)]@ rest)
4.465 + in order_mp (nth_drop 0 a)
4.466 + (order_mp' (nth_drop (length ordered - 1) ordered) [nth a 0] rest)
4.467 + end
4.468 + in order_mp (nth_drop 0 a) ordered
4.469 + end
4.470 +
4.471 + fun lc_multipoly (mpoly: multipoly) =
4.472 + get_coef (order_multipoly mpoly) (length (order_multipoly mpoly) - 1);
4.473 +
4.474 + (* greatest variablegroup *)
4.475 + fun deg_multipoly (mvp: multipoly) =
4.476 + get_varlist (order_multipoly mvp) (length (order_multipoly mvp) - 1)
4.477 +
4.478 + fun max_deg_var [m: monom] (x: int) = nth (get_varlist [m] 0) x |
4.479 + max_deg_var (mvp: multipoly) (x: int) =
4.480 + let
4.481 + fun max_monom (m1: monom) (m2: monom) (x: int)=
4.482 + if nth (get_varlist [m1] 0) x < nth (get_varlist [m2] 0) x then m2 else m1;
4.483 + in
4.484 + max_deg_var ([max_monom (nth(mvp) 0)( nth(mvp) 1) x] @ nth_drop 0 (nth_drop 0 mvp)) x
4.485 + end
4.486 +
4.487 + fun greater_deg (monom1: monom) (monom2: monom)=
4.488 + if (max_deg_var [monom1] (length(get_varlist [monom1] 0) - 1)) >
4.489 + (max_deg_var [monom2] (length (get_varlist [monom2] 0) - 1))
4.490 + then 1
4.491 + else if (max_deg_var [monom1] (length(get_varlist [monom1] 0) - 1)) <
4.492 + (max_deg_var [monom2] (length (get_varlist [monom2] 0) - 1))
4.493 + then 2
4.494 + else if length (get_varlist [monom1] 0) >1
4.495 + then
4.496 + greater_deg (1,(nth_drop 0 (get_varlist [monom1] 0))) (1,(nth_drop 0 (get_varlist [monom2] 0)))
4.497 + else 0
4.498 +
4.499 + infix %%/
4.500 + fun (poly: multipoly) %%/ b = (* =quotient*)
4.501 + let fun division monom b = (((get_coef [monom] 0) div' b),get_varlist [monom] 0);
4.502 + in order_multipoly(map2 division poly (replicate (length(poly)) b)) end
4.503 +
4.504 + fun cont [a: int] = a
4.505 + | cont (poly1: unipoly) = gcd_n poly1
4.506 + fun cont_multipoly [(a,_): monom] = a
4.507 + | cont_multipoly (poly: multipoly) = gcd_n (get_coeflist poly)
4.508 +
4.509 + fun evaluate_mv (mvp: multipoly) (var: int) (value: int) =
4.510 + let fun eval ([]: multipoly) (_: int) (_: int) (new: multipoly) = order_multipoly new |
4.511 + eval (mvp: multipoly) (var: int) (value: int) (new: multipoly) =
4.512 + eval (nth_drop 0 mvp) var value
4.513 + (new @ [((Integer.pow (nth (get_varlist mvp 0) var)value) * get_coef mvp 0,
4.514 + nth_drop var (get_varlist mvp 0))]);
4.515 + in eval mvp var value [] end
4.516 +
4.517 + fun multi_to_uni (mvp: multipoly) =
4.518 + if length (get_varlist mvp 0) = 1
4.519 + then let fun mtu ([]: multipoly) (uvp: unipoly) = uvp |
4.520 + mtu (mvp: multipoly) (uvp: unipoly) =
4.521 + if length uvp = (nth(get_varlist mvp 0) 0)
4.522 + then mtu (nth_drop 0 mvp) (uvp @ [get_coef mvp 0])
4.523 + else mtu mvp (uvp @ [0])
4.524 + in mtu (order_multipoly mvp) [] end
4.525 + else error "Polynom has more than one variable!";
4.526 +
4.527 + fun uni_to_multi (uvp: unipoly) =
4.528 + let fun utm ([]: unipoly) (mvp: multipoly) (_: int)= short_mv mvp
4.529 + | utm (uvp: unipoly) (mvp: multipoly) (counter: int) =
4.530 + utm (nth_drop 0 uvp) (mvp @ [((nth uvp 0),[counter])]) (counter+1)
4.531 + in utm uvp [] 0 end
4.532 +
4.533 + fun find_new_r (mvp1: multipoly) (mvp2: multipoly) (old: int) =
4.534 + let val poly1 = evaluate_mv mvp1 (length (get_varlist mvp1 0) - 2) (old + 1)
4.535 + val poly2 = evaluate_mv mvp2 (length (get_varlist mvp2 0) - 2) (old + 1);
4.536 + in
4.537 + if max_deg_var poly1 (length (get_varlist poly1 0) - 1)= max_deg_var mvp1 (length (get_varlist mvp1 0) - 1) orelse
4.538 + max_deg_var poly2 (length (get_varlist poly2 0) - 1)= max_deg_var mvp2 (length (get_varlist mvp2 0) - 1)
4.539 + then old + 1
4.540 + else find_new_r (mvp1) (mvp2) (old + 1)
4.541 + end
4.542 +
4.543 + fun mult_with_new_var ([]: multipoly) (_: unipoly) (_: int) = []
4.544 + | mult_with_new_var (mvp: multipoly) (uvp: unipoly) (order: int) =
4.545 + let fun mult ([]: multipoly) (_: unipoly) (_: int) (new: multipoly) (_: int) = short_mv new
4.546 + | mult (mvp: multipoly) (uvp: unipoly) (order: int) (new: multipoly) (_: int) =
4.547 + let fun mult' (_: multipoly) ([]: unipoly) (_: int) (new': multipoly) (_) = order_multipoly new'
4.548 + | mult' (mvp': multipoly) (uvp': unipoly) (order: int) (new': multipoly) (c2': int) =
4.549 + let val (first, last) = chop order (get_varlist mvp' 0)
4.550 + in mult' mvp' (nth_drop 0 uvp') order
4.551 + (new' @ [(((get_coef mvp' 0) * (nth uvp' 0)),(first @ [c2'] @ last))]) (c2'+1) end
4.552 + in mult (nth_drop 0 mvp) uvp order (new @ (mult' mvp uvp order [] 0)) (0) end
4.553 + in mult mvp uvp order [] 0 end
4.554 +
4.555 + fun greater_deg_multipoly (var1: int list) (var2: int list) =
4.556 + if var1 = [] andalso var2 =[] then 0
4.557 + else if (nth var1 (length var1 - 1)) = (nth var2 (length var1 - 1) )
4.558 + then greater_deg_multipoly (nth_drop (length var1 - 1) var1) (nth_drop (length var1 - 1) var2)
4.559 + else if (nth var1 (length var1 - 1)) > (nth var2 (length var1 - 1))
4.560 + then 1 else 2 ;
4.561 +
4.562 + infix %%+%%
4.563 + fun ([]: multipoly) %%+%% (mvp2: multipoly) = mvp2
4.564 + | (mvp1: multipoly) %%+%% ([]: multipoly) = mvp1
4.565 + | (mvp1: multipoly) %%+%% (mvp2: multipoly) =
4.566 + let fun plus ([]: multipoly) (mvp2: multipoly) (new: multipoly) = order_multipoly (new @ mvp2)
4.567 + | plus (mvp1: multipoly) ([]: multipoly) (new: multipoly) = order_multipoly (new @ mvp1)
4.568 + | plus (mvp1: multipoly) (mvp2: multipoly) (new: multipoly) =
4.569 + if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 0
4.570 + then plus (nth_drop 0 mvp1) (nth_drop 0 mvp2)
4.571 + (new @ [(((get_coef mvp1 0) + (get_coef mvp2 0)), get_varlist mvp1 0)])
4.572 + else if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 1
4.573 + then plus mvp1 (nth_drop 0 mvp2) (new @ [nth mvp2 0])
4.574 + else plus (nth_drop 0 mvp1) mvp2 (new @ [nth mvp1 0])
4.575 + in plus mvp1 mvp2 [] end
4.576 +
4.577 + infix %%-%%
4.578 + fun (mvp1: multipoly) %%-%% (mvp2: multipoly) =
4.579 + let fun neg ([]: multipoly) (new: multipoly) = order_multipoly new
4.580 + | neg (mvp: multipoly) (new: multipoly) =
4.581 + neg (nth_drop 0 mvp) (new @ [(((get_coef mvp 0) * ~1), get_varlist mvp 0)])
4.582 + val neg_mvp2 = neg mvp2
4.583 + in mvp1 %%+%% (neg_mvp2 []) end
4.584 +
4.585 + infix %%*%%
4.586 + fun (mvp1: multipoly) %%*%% (mvp2: multipoly) =
4.587 + let fun mult ([]: multipoly) (_: multipoly) (_: multipoly) (new: multipoly) = order_multipoly new
4.588 + | mult (mvp1: multipoly) ([]: multipoly) (regular_mvp2: multipoly) (new: multipoly) = mult (nth_drop 0 mvp1) regular_mvp2 regular_mvp2 new
4.589 + | mult (mvp1: multipoly) (mvp2: multipoly) (regular_mvp2: multipoly) (new: multipoly) =
4.590 + mult mvp1 (nth_drop 0 mvp2) regular_mvp2 (new @ [(((get_coef mvp1 0) * (get_coef mvp2 0)), ((get_varlist mvp1 0) %+% (get_varlist mvp2 0)))])
4.591 + in if (length mvp1) > (length mvp2) then mult mvp1 mvp2 mvp2 [] else mult mvp2 mvp1 mvp1 [] end
4.592 +
4.593 + (* if poly1|poly2 *)
4.594 + infix %%|%%
4.595 + fun [(coef2, var2): monom] %%|%% [(coef1, var1): monom] =
4.596 + ( (listgreater var1 var2)
4.597 + andalso (coef1 mod coef2) = 0)
4.598 + | (_: multipoly) %%|%% [(0,_)]= true
4.599 + | (poly2: multipoly) %%|%% (poly1: multipoly) =
4.600 + if [nth poly2 (length poly2 - 1)] %%|%% [nth poly1 (length poly1 - 1)]
4.601 + then poly2 %%|%% (poly1 %%-%%
4.602 + ([((get_coef poly1 (length poly1 - 1)) div (get_coef poly2 (length poly2 - 1)),
4.603 + (get_varlist poly1 (length poly1 - 1)) %-%(get_varlist poly2 (length poly2 - 1)))] %%*%%
4.604 + poly2))
4.605 + else false
4.606 +
4.607 +(*subsection {* Newtoninterpolation *}*)
4.608 + (* first step *)
4.609 + fun newton_first (x: int list) (f: multipoly list) (order: int) =
4.610 + let val polynom =(add_var_to_multipoly (nth f 0) order) %%+%%
4.611 + ((mult_with_new_var (((nth f 1)%%-%% (nth f 0))%%/
4.612 + (nth x 1) - (nth x 0))) [(nth x 0) * ~1, 1] order);
4.613 + val new_value_poly = [(nth x 0) * ~1, 1];
4.614 + val steps = [((nth f 1) %%-%%(nth f 0))%%/ ((nth x 1) - (nth x 0))]
4.615 + in (polynom, new_value_poly, steps) end
4.616 +
4.617 + fun newton_mv (x: int list) (f: multipoly list) (steps: multipoly list) (t: unipoly) (p: multipoly) (order: int) =
4.618 + if length x = 2
4.619 + then let val (polynom, new_value_poly, steps) = newton_first x [(nth f 0), (nth f 1)] order
4.620 + in (polynom, new_value_poly, steps) end
4.621 + else let val new_value_poly = multi_to_uni((uni_to_multi t) %%*%% (uni_to_multi [(nth x (length x - 2) )* ~1, 1]));
4.622 + val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))];
4.623 + fun next_step ([]: multipoly list) (new_steps: multipoly list) (_: int list) = new_steps
4.624 + | next_step (steps: multipoly list) (new_steps: multipoly list) (x': int list) =
4.625 + next_step (nth_drop 0 steps)
4.626 + (new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0)) %%/
4.627 + ((nth x' (length x' - 1)) - (nth x' (length x' - 3))))])
4.628 + ( nth_drop (length x' - 2) x')
4.629 + val steps = next_step steps new_steps x;
4.630 + val polynom' = (p %%+%% (mult_with_new_var (nth steps (length steps - 1)) new_value_poly order));
4.631 + in (order_multipoly(polynom'), new_value_poly, steps) end;
4.632 +
4.633 +"=========== code for [1] p.95: multivariate ============";
4.634 +"=========== code for [1] p.95: multivariate ============";
4.635 +"=========== code for [1] p.95: multivariate ============";
4.636 +(*subsection {* GCD_MODm algorithm, code for [1] p.95: multivariate *}*)
4.637 +
4.638 +fun GCD_MODm a b n s r=
4.639 + if greater_var (deg_multipoly b) (deg_multipoly a) then GCD_MODm b a n s r else
4.640 +(*0*) if s = 0
4.641 + then uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %*
4.642 + (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))
4.643 + else
4.644 + let
4.645 +(*1*) val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
4.646 +(*2*) fun GOTO2 a b n s M r_list steps = (*==============================*)
4.647 + let
4.648 + val r = find_new_r a b r;
4.649 + val r_list = r_list @ [r];
4.650 + val g_r = GCD_MODm (order_multipoly (evaluate_mv a (s- 1) r)) ( order_multipoly (evaluate_mv b (s- 1) r)) (n- 1) (s- 1) 0
4.651 + val steps = steps @ [g_r];
4.652 +(*3*) fun GOTO3 a b n s M g_r r r_list steps = (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)
4.653 + let
4.654 + val m = 1;
4.655 + fun WHILE a b n s M m r r_list steps g g_n mult = (* ----------------------- *)
4.656 + if m > M
4.657 + then
4.658 + if g_n %%|%% a andalso g_n %%|%% b
4.659 + then g_n
4.660 + else GOTO2 a b n s M r_list steps
4.661 + else
4.662 + let
4.663 + val r = find_new_r a b r;
4.664 + val r_list = r_list @ [r];
4.665 + val g_r = GCD_MODm (order_multipoly (evaluate_mv a (s- 1) r))
4.666 + (order_multipoly (evaluate_mv b (s- 1) r)) (n- 1) (s- 1) 0
4.667 + in
4.668 + if greater_var (deg_multipoly g) (deg_multipoly g_r)
4.669 + then GOTO3 a b n s M g_r r r_list steps
4.670 + else if (deg_multipoly g)= (deg_multipoly g_r)
4.671 + then
4.672 + let
4.673 + val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1)
4.674 + in if (nth steps (length steps - 1)) = (cero_multipoly s)
4.675 + then WHILE a b n s M (M+1) r r_list steps g_r g_n new
4.676 + else WHILE a b n s M (m+1) r r_list steps g_r g_n new
4.677 + end
4.678 + else WHILE a b n s M (m+1) r r_list steps g g_n mult
4.679 + end (* WHILE *)
4.680 + in (* GOTO3*) (* ----------------------- *)
4.681 + WHILE a b n s M m r r_list steps g_r ( cero_multipoly (s+1)) [1]
4.682 + end (*GOTO3*)
4.683 + in (* GOTO2*) (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)
4.684 + GOTO3 a b n s M g_r r r_list steps
4.685 + end (*GOTO2*)
4.686 + in (*==============================*)
4.687 + GOTO2 a b n s M [] []
4.688 + end (*end 0*);
4.689 +
4.690 +
4.691 +(******************************************** tests ********************************************)
4.692 +(******************************************** tests ********************************************)
4.693 +(******************************************** tests ********************************************)
4.694 +
4.695 +"----------- fun mod_inv --------------------------------";
4.696 +"----------- fun mod_inv --------------------------------";
4.697 +"----------- fun mod_inv--- -----------------------------";
4.698 +"~~~~~ fun mod_inv, args:"; val (r, m) = (5,7);
4.699 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (r, r, m, 1);
4.700 +rinv < m; (*=true*)
4.701 +r mod m = 1; (*=false*)
4.702 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
4.703 +rinv < m;(*=true*)
4.704 +r mod m = 1; (*=false*)
4.705 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
4.706 +rinv < m;(*=true*)
4.707 +r mod m = 1; (*=true*)
4.708 +"~~~~~ to mod_inv return val:"; val (rinv) = (rinv);
4.709 +
4.710 +if mod_inv 5 7 = 3 then () else error "mod_inv 5 7 = 3 changed";
4.711 +if mod_inv 3 7 = 5 then () else error "mod_inv 3 7 = 5 changed";
4.712 +if mod_inv 4 339 = 85 then () else error "mod_inv 4 339 = 85 changed";
4.713 +
4.714 +"----------- fun CRA_2 ----------------------------------";
4.715 +"----------- fun CRA_2 ----------------------------------";
4.716 +"----------- fun CRA_2 ----------------------------------";
4.717 +
4.718 +if CRA_2(17,9,3,4) = 35 then () else error "CRA_2(17,9,3,4)=35 changed";
4.719 +if CRA_2(7,2,6,11) = 17 then () else error "CRA_2(7,2,6,11) = 17 changed";
4.720 +
4.721 +"----------- fun lc -------------------------------------";
4.722 +"----------- fun lc -------------------------------------";
4.723 +"----------- fun lc -------------------------------------";
4.724 +if lc [3,4,5,6] = 6 then () else error "lc (3,4,5,6) = 6 changed" ;
4.725 +if lc [3,4,5,6,0] = 6 then () else error "lc (3,4,5,6,0) = 6 changed" ;
4.726 +
4.727 +"----------- fun deg ------------------------------------";
4.728 +"----------- fun deg ------------------------------------";
4.729 +"----------- fun deg ------------------------------------";
4.730 +if deg [3,4,5,6] = 3 then () else error "lc (3,4,5,6) = 3 changed" ;
4.731 +if deg [3,4,5,6,0] = 3 then () else error "lc (3,4,5,6) = 6 changed";
4.732 +
4.733 +"----------- fun primeList ------------------------------";
4.734 +"----------- fun primeList ------------------------------";
4.735 +"----------- fun primeList ------------------------------";
4.736 +if primeList 1 = [2] then () else error " primeList 1 = [2] changed";
4.737 +if primeList 3 = [2,3] then () else error " primeList 3 = [2,3] changed";
4.738 +if primeList 6 = [2,3,5,7] then () else error " primeList 6 = [2,3,5,7] changed";
4.739 +if primeList 7 = [2,3,5,7] then () else error " primeList 7 = [2,3,5,7] changed";
4.740 +
4.741 +"----------- fun find_next_prime_not_divide -----------------";
4.742 +"----------- fun find_next_prime_not_divide -----------------";
4.743 +"----------- fun find_next_prime_not_divide -----------------";
4.744 +if find_next_prime_not_divide 15 2 = 7 then () else error "findNextPrimeNotdivide 15 2 = 7 changed";
4.745 +if find_next_prime_not_divide 90 1 = 7 then () else error "findNextPrimeNotdivide 90 1 = 7 changed";
4.746 +
4.747 +"----------- fun poly_mod -------------------------------";
4.748 +"----------- fun poly_mod -------------------------------";
4.749 +"----------- fun poly_mod -------------------------------";
4.750 +if ([5,4,7,8,1] poly_mod 5) = [0, 4, 2, 3, 1] then ()
4.751 +else error "[5,4,7,8,1] poly_mod 5 = [0, 4, 2, 3, 1] changed" ;
4.752 +if ([5,4,~7,8,~1] poly_mod 5) = [0, 4, 3, 3, 4] then ()
4.753 +else error "[5,4,~7,8,~1] poly_mod 5 = [0, 4, 3, 3, 4] changed" ;
4.754 +
4.755 +"----------- fun %* ------------------------------";
4.756 +"----------- fun %* ------------------------------";
4.757 +"----------- fun %* ------------------------------";
4.758 +if ([5,4,7,8,1] %* 5) = [25, 20, 35, 40, 5] then ()
4.759 +else error "[5,4,7,8,1] %* 5 = [25, 20, 35, 40, 5] changed";
4.760 +if ([5,4,~7,8,~1] %* 5) = [25, 20, ~35, 40, ~5] then ()
4.761 +else error "[5,4,~7,8,~1] %* 5 = [25, 20, ~35, 40, ~5] changed";
4.762 +
4.763 +"----------- fun %/ -------------------------------";
4.764 +"----------- fun %/ -------------------------------";
4.765 +"----------- fun %/ -------------------------------";
4.766 +if ([4,3,2,5,6] %/ 3) =[1, 1, 0, 1, 2] then ()
4.767 +else error "%/ [4,3,2,5,6] 3 =[1, 1, 0, 1, 2] changed";
4.768 +if ([4,3,2,0] %/ 3) =[1, 1, 0, 0] then ()
4.769 +else error "%/ [4,3,2,0] 3 =[1, 1, 0, 0] changed";
4.770 +
4.771 +"----------- fun %-% ------------------------";
4.772 +"----------- fun %-% ------------------------";
4.773 +"----------- fun %-% ------------------------";
4.774 +if ([8,~7,0,1] %-% [~2,2,3,0])=[10,~9,~3,1] then ()
4.775 +else error "[8,~7,0,1] %-% [~2,2,3,0]=[10,~9,~3,1] changed";
4.776 +if ([8,7,6,5,4] %-% [2,2,3,1,1])=[6,5,3,4,3] then ()
4.777 +else error "[8,7,6,5,4] %-% [2,2,3,1,1]=[6,5,3,4,3] changed";
4.778 +
4.779 +"----------- fun %+% ------------------------------";
4.780 +"----------- fun %+% ------------------------------";
4.781 +"----------- fun %+% ------------------------------";
4.782 +if ([8,~7,0,1] %+% [~2,2,3,0])=[6,~5,3,1] then ()
4.783 +else error "[8,~7,0,1] %+% [~2,2,3,0]=[6,~5,3,1] changed";
4.784 +if ([8,7,6,5,4] %+% [2,2,3,1,1])=[10,9,9,6,5] then ()
4.785 +else error "[8,7,6,5,4] %+% [2,2,3,1,1]=[10,9,9,6,5] changed";
4.786 +
4.787 +"----------- fun CRA_2_poly -----------------------------";
4.788 +"----------- fun CRA_2_poly -----------------------------";
4.789 +"----------- fun CRA_2_poly -----------------------------";
4.790 +if (CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()
4.791 +else error "CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";
4.792 +
4.793 +"----------- fun mod_div --------------------------------";
4.794 +"----------- fun mod_div --------------------------------";
4.795 +"----------- fun mod_div --------------------------------";
4.796 +if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";
4.797 +if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";
4.798 +if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";
4.799 +
4.800 +"----------- fun lc_of_unipoly_not_0 --------------------";
4.801 +"----------- fun lc_of_unipoly_not_0 --------------------";
4.802 +"----------- fun lc_of_unipoly_not_0 --------------------";
4.803 +if lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] then ()
4.804 +else error "lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] changed";
4.805 +if lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0,1, 2, 3, 4, 5] then ()
4.806 +else error "lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0, 1, 2, 3, 4, 5] changed";
4.807 +
4.808 +"----------- fun %|% -------------------------------";
4.809 +"----------- fun %|% -------------------------------";
4.810 +"----------- fun %|% -------------------------------";
4.811 +"~~~~~ fun %|% , args:"; val (poly1, poly2) = ([9,12,4],[3,2]);
4.812 +val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
4.813 +val c = lc poly1 div lc b;
4.814 +val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
4.815 +rest = [];(*=false*)
4.816 +c<>0 andalso length rest >= length poly2; (*=true*)
4.817 +"~~~~~ fun %|% , args:"; val (poly1, poly2) = (rest, poly2);
4.818 +val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
4.819 +val c = lc poly1 div lc b;
4.820 +val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
4.821 +rest = [];(*=true*)
4.822 +"~~~~~ to return val:"; val (divide) = (true);
4.823 +
4.824 +if ([4] %|% [6]) = false then () else error "[4] %|% [6] = false changed";
4.825 +if ( [8] %|%[16,0]) = true then () else error "[8] %|%[16,0] = true changed";
4.826 +if ([3,2] %|%[0,0,9,12,4] ) = true then ()
4.827 +else error "[3,2] %|%[0,0,9,12,4] = true changed";
4.828 +if ([8,0] %|% [16]) = true then () else error "[8,0] %|% [16] = true changed";
4.829 +
4.830 +"----------- fun mod_poly_gcd ------------------------";
4.831 +"----------- fun mod_poly_gcd ------------------------";
4.832 +"----------- fun mod_poly_gcd ------------------------";
4.833 +
4.834 +if ( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7)=[[2, 6, 0, 2, 6], [0]]
4.835 +then () else error "( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [ [5, 1, 0, 5, 1], [0]] changed";
4.836 +if ( mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7)=[[2, 6, 0, 2, 6], [0]] then ()
4.837 +else error "mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7=[[1, 3, 0, 1, 3], [0]] changed";
4.838 +if mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] then ()
4.839 +else error " mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] changed";
4.840 +
4.841 +"~~~~~ fun mod_poly_gcd , args:";
4.842 +val (a,b,m) = ([ ~13, 2,12],[ ~14, 2],5);
4.843 +val moda = a poly_mod m;
4.844 +val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
4.845 +val c = mod_div (lc moda) (lc modb) m;
4.846 +val rest = lc_of_unipoly_not_0 (moda %-% (modb %* c) poly_mod m) poly_mod m;
4.847 +rest = []; (*=false*)
4.848 +length rest < length b; (*=false*)
4.849 +"~~~~~ fun mod_poly_gcd , args:";
4.850 +val (a,b,m,d) = (rest, b , m , [c]);
4.851 +val moda = a poly_mod m
4.852 +val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
4.853 +val c = mod_div (lc moda) (lc modb) m
4.854 +val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
4.855 +rest = [];(*=flase*)
4.856 +length rest < length b; (*=true*)
4.857 +"~~~~~ fun mod_poly_gcd , args:";
4.858 +val (a,b,m,d) = (b, rest, m, [c] @ d);
4.859 +val moda = a poly_mod m
4.860 +val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
4.861 +val c = mod_div (lc moda) (lc modb) m
4.862 +val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
4.863 +rest = [];(*=flase*)
4.864 +length rest < length b; (*=false*)
4.865 +"~~~~~ fun mod_poly_gcd , args:";
4.866 +val (a,b,m,d) = (b, rest, m, [c] @ d);
4.867 +val moda = a poly_mod m
4.868 +val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
4.869 +val c = mod_div (lc moda) (lc modb) m
4.870 +val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
4.871 +rest = [];(*=true*)
4.872 +"~~~~~ to return val:"; val ([b, rest, lsit])=[b, [0], [c] @ d];
4.873 +
4.874 +
4.875 +"----------- fun normirt_polynom ------------------------";
4.876 +"----------- fun normirt_polynom ------------------------";
4.877 +"----------- fun normirt_polynom ------------------------";
4.878 +if (normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] then ()
4.879 +else error "(normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] changed"
4.880 +
4.881 +"----------- fun poly_centr -----------------------------";
4.882 +"----------- fun poly_centr -----------------------------";
4.883 +"----------- fun poly_centr -----------------------------";
4.884 +if (poly_centr [7,3,5,8,1,3] 10) = [~3, 3, 5, ~2, 1, 3] then ()
4.885 +else error "poly_centr [7,3,5,8,1,3] 10 = [~3, 3, 5, ~2, 1, 3] cahnged";
4.886 +
4.887 +"----------- fun mod_gcd_p -------------------------------";
4.888 +"----------- fun mod_gcd_p -------------------------------";
4.889 +"----------- fun mod_gcd_p -------------------------------";
4.890 +if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5) = [1, 1, 1, 1] then ()
4.891 +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";
4.892 +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 ()
4.893 +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";
4.894 +
4.895 +
4.896 +"----------- fun gcd_n ----------------------------";
4.897 +"----------- fun gcd_n ----------------------------";
4.898 +"----------- fun gcd_n ----------------------------";
4.899 +if gcd_n [3,9,12,21] = 3 then () else error "gcd_n [3,9,12,21] = 3 changed";
4.900 +if gcd_n [12,16,32,44] = 4 then () else error "gcd_n [12,16,32,44] = 4 changed";
4.901 +if gcd_n [4,5,12] = 1 then () else error "gcd_n [4,5,12] = 1 changed";
4.902 +
4.903 +"----------- fun pp -------------------------------------";
4.904 +"----------- fun pp -------------------------------------";
4.905 +"----------- fun pp -------------------------------------";
4.906 +if pp [12,16,32,44] = [3, 4, 8, 11] then () else error "pp [12,16,32,44] = [3, 4, 8, 11] changed";
4.907 +if pp [4,5,12] = [4, 5, 12] then () else error "pp [4,5,12] = [4, 5, 12] changed";
4.908 +
4.909 +"----------- fun GCD_MOD --------------------------------";
4.910 +"----------- fun GCD_MOD --------------------------------";
4.911 +"----------- fun GCD_MOD --------------------------------";
4.912 +
4.913 +"~~~~~ fun GCD_MOD a b , args:"; val (a, b) = ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2]);
4.914 +val d = abs (Integer.gcd (lc a) (lc b));
4.915 +val M =2*d* landau_mignotte_bound a b;
4.916 +(*val g = GOTO2 a b d (*p*)1 M*)
4.917 + "~~~~~ fun GOTO2 a b d M p , args:"; val (a, b, d, p, M) = (a,b,d,3,M);
4.918 + val p = find_next_prime_not_divide d p
4.919 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
4.920 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.921 + (*val (p, g) = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)*)
4.922 + "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
4.923 + (deg g_p) = 0; (*=false*)
4.924 + val P = p
4.925 + val g = g_p;
4.926 + (*(p, g) = WHILE a b d M P g (* << 1.Mal -----*)*)
4.927 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,p,g);
4.928 + P > M ;(*false*)
4.929 + val p = find_next_prime_not_divide d p
4.930 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
4.931 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.932 + deg g_p < deg g; (*=fasle*)
4.933 + deg g_p = deg g; (*false*);
4.934 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.935 + P > M ;(*false*)
4.936 + val p = find_next_prime_not_divide d p
4.937 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
4.938 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.939 + deg g_p < deg g; (*=false*)
4.940 + deg g_p = deg g; (*=true*)
4.941 + val g = CRA_2_poly (P,p)(g,g_p)
4.942 + val P = P*p;
4.943 + val g = poly_centr(g poly_mod P)P;
4.944 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.945 + P > M ;(*false*)
4.946 + val p = find_next_prime_not_divide d p
4.947 + val c_p = normirt_polynom (mod_gcd_p a b p) p
4.948 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.949 + deg g_p < deg g; (*=true*)
4.950 + "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
4.951 + (deg g_p) = 0; (*false*)
4.952 + val P = p
4.953 + val g = g_p;
4.954 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.955 + P > M ;(*false*)
4.956 + val p = find_next_prime_not_divide d p
4.957 + val c_p = normirt_polynom (mod_gcd_p a b p) p
4.958 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.959 + deg g_p < deg g; (*=fasle*)
4.960 + deg g_p = deg g;(*true*)
4.961 + val g = CRA_2_poly (P,p)(g,g_p)
4.962 + val P = P*p;
4.963 + val g = poly_centr(g poly_mod P)P;
4.964 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.965 + P > M ;(*false*)
4.966 + val p = find_next_prime_not_divide d p
4.967 + val c_p = normirt_polynom (mod_gcd_p a b p) p
4.968 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.969 + deg g_p < deg g; (*=false*)
4.970 + deg g_p = deg g;(*true*)
4.971 + val g = CRA_2_poly (P,p)(g,g_p)
4.972 + val P = P*p;
4.973 + val g = poly_centr(g poly_mod P)P;
4.974 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.975 + P > M ;(*false*)
4.976 + val p = find_next_prime_not_divide d p
4.977 + val c_p = normirt_polynom (mod_gcd_p a b p) p
4.978 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
4.979 + deg g_p < deg g; (*=false*)
4.980 + deg g_p = deg g;(*true*)
4.981 + val g = CRA_2_poly (P,p)(g,g_p)
4.982 + val P = P*p;
4.983 + val g = poly_centr(g poly_mod P)P;
4.984 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
4.985 + P > M ;(*true*)
4.986 + "~~~~~ to return WHILE val:"; val (g,p) = (g,p);
4.987 + "~~~~~ to return GOTO3 val:"; val (g,p) = (g,p);
4.988 + val g = pp g;
4.989 + (g %|% a) andalso (g %|% b);(*=true*)
4.990 +"~~~~~ to return GOTO2 val:"; val (g) = (g);
4.991 +"~~~~~ to return GCD_MOD val:"; val (g) = (g);
4.992 +
4.993 +"----------- fun GCD_MOD---------------------------------";
4.994 +"----------- fun GCD_MOD --------------------------------";
4.995 +"----------- fun GCD_MOD --------------------------------";
4.996 +if GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] then ()
4.997 +else error "GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] changed";
4.998 +
4.999 +"------------------------------------------------------------";
4.1000 +"--------------------- Multivariate Case --------------------";
4.1001 +"------------------------------------------------------------";
4.1002 +
4.1003 +"----------- fun get_coef --------------------------- ---";
4.1004 +"----------- fun get_coef -------------------------------";
4.1005 +"----------- fun get_coef -------------------------------";
4.1006 +if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 then () else
4.1007 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 changed";
4.1008 +if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 then () else
4.1009 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 changed";
4.1010 +if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 then () else
4.1011 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 changed";
4.1012 +if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 then () else
4.1013 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 changed";
4.1014 +if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 then () else
4.1015 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 changed";
4.1016 +
4.1017 +
4.1018 +"----------- fun get_varlist ----------------------------";
4.1019 +"----------- fun get_varlist ----------------------------";
4.1020 +"----------- fun get_varlist ----------------------------";
4.1021 +if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] then () else
4.1022 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] changed";
4.1023 +if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] then () else
4.1024 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] changed";
4.1025 +if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] then () else
4.1026 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] changed";
4.1027 +if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] then () else
4.1028 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] changed";
4.1029 +if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] then () else
4.1030 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] changed";
4.1031 +
4.1032 +"----------- fun get_coeflist ---------------------------";
4.1033 +"----------- fun get_coeflist ---------------------------";
4.1034 +"----------- fun get_coeflist ---------------------------";
4.1035 +if get_coeflist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] = [1,2,3,4,5] then () else
4.1036 +error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] = [1,2,3,4,5] changed";
4.1037 +
4.1038 +"----------- fun add_var_to_multipoly -------------------";
4.1039 +"----------- fun add_var_to_multipoly -------------------";
4.1040 +"----------- fun add_var_to_multipoly -------------------";
4.1041 +if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] then () else
4.1042 +error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] changed";
4.1043 +if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3])] then () else
4.1044 +error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3]) changed";
4.1045 +if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] then () else
4.1046 +error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] changed";
4.1047 +
4.1048 +"----------- fun cero_multipoly -------------------------";
4.1049 +"----------- fun cero_multipoly -------------------------";
4.1050 +"----------- fun cero_multipoly -------------------------";
4.1051 +if cero_multipoly 1 = [(0, [0])] then () else
4.1052 +error "cero_multipoly 1 = [(0, [0])] changed";
4.1053 +if cero_multipoly 5 = [(0, [0, 0, 0, 0, 0])] then () else
4.1054 +error "cero_multipoly 1 = [(0, [0, 0, 0, 0, 0])] changed";
4.1055 +
4.1056 +"----------- fun short_mv -------------------------------";
4.1057 +"----------- fun short_mv -------------------------------";
4.1058 +"----------- fun short_mv -------------------------------";
4.1059 +"~~~~~ fun short_mv, args:"; val (mvp) = ([(~12, [1]), (2, [1]), (4, [0])]);
4.1060 +"~~~~~ fun short , args:"; val (mvp, new) = (mvp,([]:multipoly));
4.1061 +length mvp =1(*false*);
4.1062 + get_varlist mvp 0 = get_varlist mvp 1; (* true*)
4.1063 +"~~~~~ fun short , args:"; val (mvp, new) =
4.1064 +([(get_coef mvp 0 + get_coef mvp 1,get_varlist mvp 0)] @ nth_drop 0 (nth_drop 0 mvp), new );
4.1065 +length mvp =1(*false*);
4.1066 +get_varlist mvp 0 = get_varlist mvp 1;(*false*)
4.1067 +"~~~~~ fun short , args:"; val (mvp, new) = ((nth_drop 0 mvp),new @ [nth mvp 0]);
4.1068 +length mvp =1(*true*);
4.1069 +"~~~~~ to return val:"; val (shortform) = (new @ mvp);
4.1070 +
4.1071 +if short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])]
4.1072 +then () else error "short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])] changed"
4.1073 +
4.1074 +"----------- fun greater_var ----------------------------";
4.1075 +"----------- fun greater_var ----------------------------";
4.1076 +"----------- fun greater_var ----------------------------";
4.1077 +if greater_var [3] [4]
4.1078 +then error " greater_var [3] [4] changed = false" else ();
4.1079 +if greater_var [1,2] [0,3]
4.1080 +then error " greater_var [1,2] [0,3] changed = false" else ();
4.1081 +if greater_var [1,2] [3,0]
4.1082 +then () else error " greater_var [1,2] [3,0] = true changed";
4.1083 +if greater_var [1,2] [1,2]
4.1084 +then error " greater_var [1,2] [1,2] changed = false" else ();
4.1085 +
4.1086 +"----------- fun order_multipoly ------------------------";
4.1087 +"----------- fun order_multipoly ------------------------";
4.1088 +"----------- fun order_multipoly ------------------------";
4.1089 +if order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])]
4.1090 +then () else error "order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])] changed"
4.1091 +
4.1092 +"----------- fun lc_multipoly ---------------------------";
4.1093 +"----------- fun lc_multipoly ---------------------------";
4.1094 +"----------- fun lc_multipoly ---------------------------";
4.1095 +if lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1
4.1096 +then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1 changed";
4.1097 +if lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2
4.1098 +then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2 changed";
4.1099 +
4.1100 +"----------- fun deg_multipoly -------------------------";
4.1101 +"----------- fun deg_multipoly -------------------------";
4.1102 +"----------- fun deg_multipoly -------------------------";
4.1103 +if deg_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [1,2]
4.1104 +then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1 changed";
4.1105 +if deg_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = [1,2]
4.1106 +then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2 changed";
4.1107 +
4.1108 +
4.1109 +
4.1110 +"----------- fun max_deg_var ----------------------------";
4.1111 +"----------- fun max_deg_var ----------------------------";
4.1112 +"----------- fun max_deg_var ----------------------------";
4.1113 +if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 then ()
4.1114 +else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 changed";
4.1115 +if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 then ()
4.1116 +else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 changed";
4.1117 +
4.1118 +"----------- infix %%/ ----------------------------------";
4.1119 +"----------- infix %%/ ----------------------------------";
4.1120 +"----------- infix %%/ ----------------------------------";
4.1121 +if ([(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2) = [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] then ()
4.1122 +else error "[(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2 = [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] changed";
4.1123 +
4.1124 +"----------- fun cont -----------------------------------";
4.1125 +"----------- fun cont -----------------------------------";
4.1126 +"----------- fun cont -----------------------------------";
4.1127 +if cont [4,8,12,~2,0,~4] = 2 then () else error " cont [4,8,12,~2,0,~4] = 2 changed";
4.1128 +if cont [3,6,9,19] = 1 then () else error " cont [3,6,9,19] = 1 changed";
4.1129 +
4.1130 +"----------- fun cont_multipoly -------------------------";
4.1131 +"----------- fun cont_multipoly -------------------------";
4.1132 +"----------- fun cont_multipoly -------------------------";
4.1133 +if cont_multipoly [(3,[1,2,3,1])] = 3 then () else error " cont_multipoly [(3,[1,2,3,1])] = 3 changed";
4.1134 +if cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 then ()
4.1135 +else error "cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 changed";
4.1136 +if cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 then ()
4.1137 +else error "cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 changed";
4.1138 +
4.1139 +"----------- fun evaluate_mv ----------------------------";
4.1140 +"----------- fun evaluate_mv ----------------------------";
4.1141 +"----------- fun evaluate_mv ----------------------------";
4.1142 +if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])]
4.1143 +then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])] changed";
4.1144 +if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 = [(2, [0]), (~2, [2])]
4.1145 +then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 =[(2, [0]), (~2, [2])] changed";
4.1146 +if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])]
4.1147 +then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])] changed";
4.1148 +if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])]
4.1149 +then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])] changed";
4.1150 +
4.1151 +"----------- fun mutli_to_uni ---------------------------";
4.1152 +"----------- fun mutli_to_uni ---------------------------";
4.1153 +"----------- fun mutli_to_uni ---------------------------";
4.1154 +if multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1]
4.1155 +then () else error " multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1] changed";
4.1156 +
4.1157 +"~~~~~ fun multi_to_uni , args:"; val (mvp) = ([(~3, [0]), (1, [0]), (3, [1]), (~6, [1]),(2,[3])]);
4.1158 +val mvp = order_multipoly mvp;
4.1159 +"~~~~~ fun short, args:"; val (mvp, uvp) = ((short_mv mvp), ([]:unipoly));
4.1160 +length uvp = (nth(get_varlist mvp 0) 0);(*true*)
4.1161 +"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
4.1162 +length uvp = (nth(get_varlist mvp 0) 0);(*true*)
4.1163 +"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
4.1164 +length uvp = (nth(get_varlist mvp 0) 0);(*false*)
4.1165 +"~~~~~ fun short, args:"; val (mvp, uvp) = (mvp, (uvp @ [0]));
4.1166 +length uvp = (nth(get_varlist mvp 0) 0);(*true*)
4.1167 +"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
4.1168 +"~~~~~ to return val:"; val (unipoly) = (uvp);
4.1169 +
4.1170 +
4.1171 +
4.1172 +"----------- fun find_new_r ----------------------------";
4.1173 +"----------- fun find_new_r ----------------------------";
4.1174 +"----------- fun find_new_r ----------------------------";
4.1175 +if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 then ()
4.1176 +else error " find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 changed";
4.1177 +if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 then ()
4.1178 +else error "find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 changed";
4.1179 +if find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 then ()
4.1180 +else error "find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 changed";
4.1181 +
4.1182 +"----------- fun newton_mv -----------------------------";
4.1183 +"----------- fun newton_mv -----------------------------";
4.1184 +"----------- fun newton_mv -----------------------------";
4.1185 +if uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] then ()
4.1186 +else error "uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] changed";
4.1187 +if uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] then ()
4.1188 +else error "uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] changed";
4.1189 +
4.1190 +"----------- fun mult_with_new_var ---------------------";
4.1191 +"----------- fun mult_with_new_var ---------------------";
4.1192 +"----------- fun mult_with_new_var ---------------------";
4.1193 +if mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] then ()
4.1194 +else error "mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] changed";
4.1195 +
4.1196 +"----------- fun greater_deg_multipoly -----------------";
4.1197 +"----------- fun greater_deg_multipoly -----------------";
4.1198 +"----------- fun greater_deg_multipoly -----------------";
4.1199 +greater_deg_multipoly [1,2,3] [1,2,4] =2;
4.1200 +
4.1201 +"----------- infix %%+%% -------------------------------";
4.1202 +"----------- infix %%+%% -------------------------------";
4.1203 +"----------- infix %%+%% -------------------------------";
4.1204 +if ([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])]
4.1205 +then () else error "([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] changed";
4.1206 +
4.1207 +"----------- infix %%-%% -------------------------------";
4.1208 +"----------- infix %%-%% -------------------------------";
4.1209 +"----------- infix %%-%% -------------------------------";
4.1210 +if ([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])]
4.1211 +then () else error "([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])] changed";
4.1212 +
4.1213 +"----------- infix %%*%% -------------------------------";
4.1214 +"----------- infix %%*%% -------------------------------";
4.1215 +"----------- infix %%*%% -------------------------------";
4.1216 +if ([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])]
4.1217 +then () else error "([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])] changed";
4.1218 +if ([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] then ()
4.1219 +else error "([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] changed";
4.1220 +if ([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0
4.1221 +then () else error
4.1222 +"([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 changed";
4.1223 +
4.1224 +"----------- fun newton_first --------------------------";
4.1225 +"----------- fun newton_first --------------------------";
4.1226 +"----------- fun newton_first --------------------------";
4.1227 +if newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 =
4.1228 +([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])
4.1229 +then () else error
4.1230 +"newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 = "
4.1231 +"([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";
4.1232 +
4.1233 +"----------- fun newton_mv -----------------------------";
4.1234 +"----------- fun newton_mv -----------------------------";
4.1235 +"----------- fun newton_mv -----------------------------";
4.1236 +if newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 =
4.1237 +([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])
4.1238 +then () else error
4.1239 +"newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 = "
4.1240 +"([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";
4.1241 +if newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 =
4.1242 +([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]])
4.1243 +then () else error
4.1244 +"newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 ="
4.1245 +" ([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]]) changed";
4.1246 +
4.1247 +
4.1248 +"~~~~~ fun newton_mv, args:"; val (x,f,steps,t,p,order) = ([1, 2, 3, 4],[[(9, [0]), (5, [1])], [(16, [0]), (7, [1])]], [[(5, [0]), (2, [1])], [(1, [0])]], [2, ~3, 1], [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])], 0 );
4.1249 +length x = 2; (* false *)
4.1250 +val new_value_poly = multi_to_uni((uni_to_multi t) %%*%% (uni_to_multi [(nth x (length x - 2) )* ~1, 1]));
4.1251 +val new_steps = [((nth f (length f - 1)) %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))) %%-%% ((nth f (length f - 2)))];
4.1252 +
4.1253 +"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (steps, new_steps, x);
4.1254 +steps = []; (*false*)
4.1255 +
4.1256 +"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0))) %%/
4.1257 + ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' - 2) x');steps = []; (*false*)
4.1258 +
4.1259 +"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0))) %%/
4.1260 + ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' - 2) x');
4.1261 +steps = []; (*true*)
4.1262 +val steps = new_steps;
4.1263 +val polynom' = p %%+%% (mult_with_new_var (nth steps (length steps - 1)) new_value_poly order);
4.1264 +
4.1265 +"----------- fun listgreater ---------------------------";
4.1266 +"----------- fun listgreater ---------------------------";
4.1267 +"----------- fun listgreater ---------------------------";
4.1268 +if listgreater [1,2,3,4,5] [1,2,3,4,5] = true then ()
4.1269 +else error " listgreater [1,2,3,4,5] [1,2,3,4,5] = true changed";
4.1270 +if listgreater [1,2,3,4] [1,2,3,5] = false then ()
4.1271 +else error "listgreater [1,2,3,4] [1,2,3,5] = false changed" ;
4.1272 +if listgreater [1,4,5,4] [0,3,4,5] = false then ()
4.1273 +else error "listgreater [1,2,3,4] [0,3,4,5] = false changed ";
4.1274 +
4.1275 +"----------- fun greater_deg_multipoly -----------------";
4.1276 +"----------- fun greater_deg_multipoly -----------------";
4.1277 +"----------- fun greater_deg_multipoly -----------------";
4.1278 +if greater_deg_multipoly [1,2,3,4,5] [1,2,3,4,5] = 0 then ()
4.1279 +else error " greater_deg_multipoly [1,2,3,4,5] [1,2,3,4,5] = 0 changed";
4.1280 +if greater_deg_multipoly [1,2,3,4] [5,2,8,1] = 1 then ()
4.1281 +else error "greater_deg_multipoly [1,2,3,4] [1,2,3,5] = 1 changed" ;
4.1282 +if greater_deg_multipoly [1,4,5,4] [0,3,4,5] = 2 then ()
4.1283 +else error "greater_deg_multipoly [1,2,3,4] [0,3,4,5] = 2 changed ";
4.1284 +
4.1285 +"----------- finfix %%|%% ------------------------------";
4.1286 +"----------- finfix %%|%% ------------------------------";
4.1287 +"----------- finfix %%|%% ------------------------------";
4.1288 +if [(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] then ()
4.1289 +else error "[(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] = true changed";
4.1290 +if [(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] then ()
4.1291 +else error "[(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = true changed";
4.1292 +if [(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])]
4.1293 +then error "[(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = false changed" else ();
4.1294 +
4.1295 +"----------- fun GCD_MODm ------------------------------";
4.1296 +"----------- fun GCD_MODm ------------------------------";
4.1297 +"----------- fun GCD_MODm ------------------------------";
4.1298 +
4.1299 +if GCD_MODm
4.1300 + [(~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.1301 + [(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 1 0
4.1302 + = [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] then () else error
4.1303 +"GCD_MODm [(~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.1304 +"[(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 1 0 "
4.1305 +"= [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] changed";
4.1306 +(* -xy +xy^2z+yz - 1*)(* xy +1*) (*=*) (*xy - 1*)
4.1307 +if GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0
4.1308 += [(1, [0, 0, 0]), (1, [1, 1, 0])] then () else error
4.1309 +"GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0 "
4.1310 +"= [(1, [0, 0, 0]), (1, [1, 1, 0])] changed";
4.1311 +
4.1312 +" ~~~ 1 ~~~";
4.1313 +"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ([(1,[1,2,1])], [(1,[1,2,1])], 3, 2, 0);
4.1314 +s = 0; (*false*)
4.1315 +val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
4.1316 + " ~~~ 1 ~~~";
4.1317 +"~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);
4.1318 +val r = find_new_r a b r;
4.1319 +val r_list = r_list @ [r];
4.1320 +val (a',b',n',s',r',r_list',steps',M') = ( a,b,n,s,r,r_list,steps,M);
4.1321 +(*g_1=*)
4.1322 + " ~~~ 1.1 ~~~";
4.1323 +"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1324 +s = 0; (*false*)
4.1325 +
4.1326 +val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
4.1327 + " ~~~ 1.1 ~~~";
4.1328 +"~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);
4.1329 +val r = find_new_r a b r;
4.1330 +val r_list = r_list @ [r];
4.1331 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1332 +(*g_2=*)
4.1333 + " ~~~ 1.1.1 ~~~";
4.1334 +"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1335 +s = 0; (*true*)
4.1336 +val g_2 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))
4.1337 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1338 +val g_r = g_2;
4.1339 +val steps = steps @ [g_r];
4.1340 + " ~~~ 1.1 ~~~";
4.1341 +"~~~~~ fun GOTO3 , args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
4.1342 +(a, b, n, s, M, g_r, r, r_list, steps);
4.1343 +val m = 1;
4.1344 + " ~~~ 1.1.W1 ~~~";
4.1345 +"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
4.1346 +(a, b, n, s, M, m,r, r_list, steps, g_r, ( cero_multipoly (s+1)), [1]);
4.1347 +; (*false*)
4.1348 +val r = find_new_r a b r;
4.1349 +val r_list = r_list @ [r];
4.1350 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1351 +(*g_3=*)
4.1352 + " ~~~ 1.1.W1.1 ~~~";
4.1353 +"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1354 +s = 0; (*true*)
4.1355 +val g_3 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1356 + " ~~~ 1.1.W1 ~~~";
4.1357 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1358 +val g_r = g_3;
4.1359 +greater_var (deg_multipoly g) (deg_multipoly g_r) (*false*);
4.1360 +(deg_multipoly g)= (deg_multipoly g_r) (*true*);
4.1361 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1362 +if g_n = [(1,[1,1])] then () else error "g_n is false" ;
4.1363 +(nth steps (length steps - 1)) = (cero_multipoly s) (* false*);
4.1364 + " ~~~ 1.1.W2 ~~~";
4.1365 +"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
4.1366 +(a, b, n, s, M, (m+1), r, r_list, steps, g_r, g_n, new);
4.1367 +m > M; (*false*)
4.1368 +val r = find_new_r a b r;
4.1369 +val r_list = r_list @ [r];
4.1370 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1371 +(*g_4=*)
4.1372 +" ~~~ 1.1.W2.1 ~~~";
4.1373 +"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1374 +s = 0; (*true*)
4.1375 +val g_r =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1376 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1377 +" ~~~ 1.1.W2 ~~~";
4.1378 +greater_var (deg_multipoly g) (deg_multipoly g_r) (*false*);
4.1379 +(deg_multipoly g)= (deg_multipoly g_r) (*true*);
4.1380 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1381 +if g_n = [(1,[1,1])] then () else error "g_n is false" ;
4.1382 +(nth steps (length steps - 1)) = (cero_multipoly s) (* true*);
4.1383 + " ~~~ 1.1.W3 ~~~";
4.1384 +"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
4.1385 +(a, b, n, s, M, (M+1), r, r_list, steps, g_r, g_n, new);
4.1386 +m > M; (*true*)
4.1387 +g_n %%|%% a andalso g_n %%|%% b(*true*);
4.1388 + " ~~~ 1 ~~~";
4.1389 +val ( a,b,n,s,r,r_list,steps,M) = (a',b',n',s',r',r_list',steps',M');
4.1390 +val g_1=g_n;
4.1391 +val g_r = g_1;
4.1392 +val steps = steps @ [g_r];
4.1393 + " ~~~ 1 ~~~";
4.1394 +"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
4.1395 +(a, b, n, s, M, g_r, r, r_list, steps);
4.1396 +val m = 1;
4.1397 + " ~~~ 1.W1 ~~~";
4.1398 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1399 +(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
4.1400 +val g' = g;
4.1401 +m > M; (*false*)
4.1402 +val r = find_new_r a b r;
4.1403 +val r_list = r_list @ [r];
4.1404 +val (a',b',n',s',r',r_list',steps',m') = ( a,b,n,s,r,r_list,steps,m);
4.1405 +(* g_2 = GCD_MODm *)
4.1406 + " ~~~ 1.W1.1 ~~~";
4.1407 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1408 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1409 +s = 0; (*false*)
4.1410 +val M = 1 + Int.min(max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
4.1411 + " ~~~ 1.W1.1 ~~~";
4.1412 +"~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =
4.1413 +(a, b, n, s, M, [], []);
4.1414 +val r = find_new_r a b r;
4.1415 +val r_list = r_list @ [r];
4.1416 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1417 + " ~~~ 1.W1.1.1 ~~~";
4.1418 + (*g_r=*)
4.1419 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1420 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1421 +s = 0; (*true*)
4.1422 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1423 + " ~~~ 1.W1.1 ~~~";
4.1424 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1425 +val steps = steps @ [g_r];
4.1426 + " ~~~ 1.W1.1 ~~~";
4.1427 +"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
4.1428 +(a, b, n, s, M, g_r, r, r_list, steps);
4.1429 +val m = 1;
4.1430 +val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);
4.1431 + " ~~~ 1.W1.1.W1 ~~~";
4.1432 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1433 +(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
4.1434 +
4.1435 +m > M; (*false*)
4.1436 +val r = find_new_r a b r;
4.1437 +val r_list = r_list @ [r];
4.1438 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1439 +(* g_r = GCD_MODm *)
4.1440 + " ~~~ 1.W1.1.W1.1 ~~~";
4.1441 + (*g_r=*)
4.1442 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1443 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1444 +s = 0; (*true*)
4.1445 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1446 + " ~~~ 1.W1.1.W1 ~~~";
4.1447 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1448 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1449 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1450 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1451 +(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
4.1452 + " ~~~ 1.W1.1.W2 ~~~";
4.1453 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1454 +(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
4.1455 +m > M; (*false*)
4.1456 +val r = find_new_r a b r;
4.1457 +val r_list = r_list @ [r];
4.1458 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1459 +(* g_r = GCD_MODm *)
4.1460 + " ~~~ 1.W1.1.W2.1 ~~~";
4.1461 + (*g_r=*)
4.1462 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1463 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1464 +s = 0; (*true*)
4.1465 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1466 + " ~~~ 1.W1.1.W2 ~~~";
4.1467 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1468 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1469 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1470 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1471 +(nth steps (length steps - 1)) = (cero_multipoly s); (*true*)
4.1472 + " ~~~ 1.W1.1.W3 ~~~";
4.1473 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1474 +(a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);
4.1475 +m>M; (*true*)
4.1476 +g_n %%|%% a andalso g_n %%|%% b; (*true*)
4.1477 + " ~~~ 1.W1.1 ~~~";
4.1478 +val ( a,b,n,s,r,r_list,g,M) = (a',b',n',s',r',r_list',g',M');
4.1479 +val g_r = g_n;
4.1480 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1481 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1482 +val (g_n,steps,m,new) = (g_n'',steps',m'',new'');
4.1483 + " ~~~ 1.W2 ~~~";
4.1484 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps new g_n (s- 1);
4.1485 +(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
4.1486 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1487 +(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
4.1488 +m > M; (*false*)
4.1489 +val r = find_new_r a b r;
4.1490 +val r_list = r_list @ [r];
4.1491 +val (a',b',n',s',r',r_list',steps',m',M',g') = ( a,b,n,s,r,r_list,steps,m,M,g);
4.1492 +(* g_2 = GCD_MODm *)
4.1493 + " ~~~ 1.W2.1 ~~~";
4.1494 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1495 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1496 +s = 0; (*false*)
4.1497 +val M = 1 + Int.min(max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
4.1498 + " ~~~ 1.W2.1 ~~~";
4.1499 +"~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =
4.1500 +(a, b, n, s, M, [], []);
4.1501 +val r = find_new_r a b r;
4.1502 +val r_list = r_list @ [r];
4.1503 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1504 + " ~~~ 1.W2.1.1 ~~~";
4.1505 + (*g_r=*)
4.1506 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1507 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1508 +s = 0; (*true*)
4.1509 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1510 + " ~~~ 1.W2.1 ~~~";
4.1511 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1512 +val steps = steps @ [g_r];
4.1513 + " ~~~ 1.W2.1 ~~~";
4.1514 +"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
4.1515 +(a, b, n, s, M, g_r, r, r_list, steps);
4.1516 +val m = 1;
4.1517 +val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);
4.1518 + " ~~~ 1.W2.1.W1 ~~~";
4.1519 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1520 +(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
4.1521 +m > M; (*false*)
4.1522 +val r = find_new_r a b r;
4.1523 +val r_list = r_list @ [r];
4.1524 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1525 +(* g_r = GCD_MODm *)
4.1526 + " ~~~ 1.W2.1.W1.1 ~~~";
4.1527 + (*g_r=*)
4.1528 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1529 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1530 +s = 0; (*true*)
4.1531 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1532 + " ~~~ 1.W2.1.W1 ~~~";
4.1533 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1534 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1535 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1536 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1537 +(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
4.1538 + " ~~~ 1.W2.1.W2 ~~~";
4.1539 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1540 +(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
4.1541 +m > M; (*false*)
4.1542 +val r = find_new_r a b r;
4.1543 +val r_list = r_list @ [r];
4.1544 +val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
4.1545 +(* g_r = GCD_MODm *)
4.1546 + " ~~~ 1.W2.1.W2.1 ~~~";
4.1547 + (*g_r=*)
4.1548 +"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
4.1549 +((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
4.1550 +s = 0; (*true*)
4.1551 +val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
4.1552 + " ~~~ 1.W2.1.W2 ~~~";
4.1553 +val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
4.1554 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1555 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1556 +val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1557 +(nth steps (length steps - 1)) = (cero_multipoly s); (*true*)
4.1558 + " ~~~ 1.W2.1.W3 ~~~";
4.1559 +"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
4.1560 +(a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);
4.1561 +m>M; (*true*)
4.1562 +g_n %%|%% a andalso g_n %%|%% b; (*true*)
4.1563 + " ~~~ 1.W2.1 ~~~";
4.1564 +val ( a,b,n,s,r,r_list,g) = (a',b',n',s',r',r_list',g');
4.1565 +val g_r = g_n;
4.1566 +greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
4.1567 +(deg_multipoly g)= (deg_multipoly g_r); (*true*)
4.1568 +val (g_n,steps,m,mult) = (g_n'',steps',m'',new'');
4.1569 + " ~~~ 1.W2 ~~~";
4.1570 + val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
4.1571 +"~~~~~ to return val:"; val (g) = (g_n);
4.1572 +
4.1573 +
4.1574 +
4.1575 +if g = [(1, [1, 2, 1])] then () else error "GCD_MODm [(1, [1, 2, 1])] [(1, [1, 2, 1])] has changes.";
4.1576 +
4.1577 +" ========================== END ========================== ";
4.1578 + fun nth _ [] = error "nth _ []" (*Isabelle2002, still saved the hours of update*)
4.1579 + | nth 1 (x::_) = x
4.1580 + | nth n (_::xs) = nth (n- 1) xs;
4.1581 +(*fun nth xs i = List.nth (xs, i); recent Isabelle: TODO update all isac code *)
4.1582 +
5.1 --- a/test/Tools/isac/Knowledge/gcd_poly.thy Fri Aug 06 12:09:06 2021 +0200
5.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
5.3 @@ -1,7 +0,0 @@
5.4 -theory gcd_poly
5.5 -imports (*"$ISABELLE_ISAC/Knowledge/GCD_Poly"*) Isac.Isac_Knowledge
5.6 -begin
5.7 -
5.8 -(*here come the tests from GCD_Poly.thy*)
5.9 -
5.10 -end
6.1 --- a/test/Tools/isac/Knowledge/gcd_poly_ml.sml Fri Aug 06 12:09:06 2021 +0200
6.2 +++ b/test/Tools/isac/Knowledge/gcd_poly_ml.sml Fri Aug 06 18:27:05 2021 +0200
6.3 @@ -23,7 +23,7 @@
6.4 "----------- fun deg_up ---------------------------------";
6.5 "----------- fun drop_lc0_up ----------------------------";
6.6 "----------- fun normalise ------------------------------";
6.7 -(* order until here: (GCD_Poly_FP =) GCD_Poly = gcd_poly *)
6.8 +(* order until here: (GCD_Poly_FP =) GCD_Poly_OLD = gcd_poly *)
6.9 "----------- fun %* -------------------------------------";
6.10 "----------- fun %/ -------------------------------------";
6.11 "----------- fun %-% ------------------------------------";
6.12 @@ -998,7 +998,7 @@
6.13 r (* = [317, 295]*)) = p1 %*/% p2;
6.14 if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed 1";
6.15
6.16 -(* demo in GCD_Poly *)
6.17 +(* demo in GCD_Poly_OLD *)
6.18 val p1 = [2, 2, 2, 2, 2];
6.19 val p2 = [1, 2, 3];
6.20 val (n (* = 27*),
7.1 --- a/test/Tools/isac/Knowledge/gcd_poly_winkler.sml Fri Aug 06 12:09:06 2021 +0200
7.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
7.3 @@ -1,1578 +0,0 @@
7.4 -(* Title: test/../gcd_poly_winkler
7.5 - Author: Diana Meindl
7.6 - Copyright (c) Diana Meindl 2011
7.7 -12345678901234567890123456789012345678901234567890123456789012345678901234567890
7.8 - 10 20 30 40 50 60 70 80
7.9 -
7.10 -Programcode according to [1] for later lookup for mathematicians.
7.11 -The tests below remain according to [1],
7.12 -while "$ISABELLE_ISAC_TEST/Tools/isac/Knowledge/gcd_poly.sml" start exactly from the same state
7.13 -and in time follows development in "$ISABELLE_ISAC/Knowledge/GCD_Poly.thy".
7.14 -
7.15 -[1] Franz Winkler. Polynomial Algorithms in Computer Algebra.
7.16 -Springer-Verlag/Wien 1996.
7.17 -*)
7.18 -
7.19 -(*fun nth _ [] = error "nth _ []" (*Isabelle2002, still saved the hours of update*)
7.20 - | nth 1 (x::_) = x
7.21 - | nth n (_::xs) = nth (n- 1) xs;*)
7.22 - fun nth xs i = List.nth (xs, i); (*recent Isabelle: TODO update all isac code *)
7.23 -"--------------------------------------------------------";
7.24 -"table of contents --------------------------------------";
7.25 -"--------------------------------------------------------";
7.26 -"----------- auxiliary functions univariate -------------";
7.27 -"=========== code for [1] p.93: univariate ==============";
7.28 -"----------- auxiliary functions multivariate -----------";
7.29 -"=========== code for [1] p.95: multivariate ============";
7.30 -"--------------------------------------------------------";
7.31 -"--- begin of univariate polynomials --------------------";
7.32 -"----------- fun mod_inv --------------------------------";
7.33 -"----------- fun CRA_2 ----------------------------------";
7.34 -"----------- fun lc -------------------------------------";
7.35 -"----------- fun deg ------------------------------------";
7.36 -"----------- fun primeList ------------------------------";
7.37 -"----------- fun find_next_prime_not_divide -------------";
7.38 -"----------- fun poly_mod -------------------------------";
7.39 -"----------- fun %* -------------------------------------";
7.40 -"----------- fun %/ -------------------------------------";
7.41 -"----------- fun %-% ------------------------------------";
7.42 -"----------- fun %+% ------------------------------------";
7.43 -"----------- fun CRA_2_poly -----------------------------";
7.44 -"----------- fun mod_div --------------------------------";
7.45 -"----------- fun lc_of_unipoly_not_0 --------------------";
7.46 -"----------- fun %|% ------------------------------------";
7.47 -"----------- fun mod_poly_gcd ---------------------------";
7.48 -"----------- fun normirt_polynom ------------------------";
7.49 -"----------- fun poly_centr -----------------------------";
7.50 -"----------- fun mod_gcd_p ------------------------------";
7.51 -"----------- fun gcd_n ----------------------------------";
7.52 -"----------- fun pp -------------------------------------";
7.53 -"----------- fun GCD_MOD---------------------------------";
7.54 -"--- begin of multivariate polynomials ------------------";
7.55 -"----------- fun get_coef --------------------------- ---";
7.56 -"----------- fun get_varlist ----------------------------";
7.57 -"----------- fun get_coeflist ---------------------------";
7.58 -"----------- fun add_var_to_multipoly -------------------";
7.59 -"----------- fun cero_multipoly -------------------------";
7.60 -"----------- fun short_mv -------------------------------";
7.61 -"----------- fun greater_var ----------------------------";
7.62 -"----------- fun order_multipoly ------------------------";
7.63 -"----------- fun lc_multipoly ---------------------------";
7.64 -"----------- fun deg_multipoly --------------------------";
7.65 -"----------- fun max_deg_var ----------------------------";
7.66 -"----------- infix %%/ ----------------------------------";
7.67 -"----------- fun cont -----------------------------------";
7.68 -"----------- fun cont_multipoly -------------------------";
7.69 -"----------- fun evaluate_mv ----------------------------";
7.70 -"----------- fun mutli_to_uni ---------------------------";
7.71 -"----------- fun find_new_r ----------------------------";
7.72 -"----------- fun newton_mv -----------------------------";
7.73 -"----------- fun mult_with_new_var ---------------------";
7.74 -"----------- fun greater_deg_multipoly -----------------";
7.75 -"----------- infix %%+%% -------------------------------";
7.76 -"----------- infix %%-%% -------------------------------";
7.77 -"----------- infix %%*%% -------------------------------";
7.78 -"----------- fun newton_first --------------------------";
7.79 -"----------- fun newton_mv -----------------------------";
7.80 -"----------- fun listgreater ---------------------------";
7.81 -"----------- finfix %%|%% -----------------------------";
7.82 -"----------- fun GCD_MODm ------------------------------";
7.83 -
7.84 -(*text {*
7.85 - Below we reference
7.86 - F. Winkler, Polynomial Algorithms. ...
7.87 - by page 93 and 95. The identifiers used in this book are re-used in this thesis
7.88 - in order to support reference (although some of these identifiers
7.89 - doe not conform with the Isabelle coding standards)
7.90 -*}*)
7.91 -
7.92 - (*default_print_depth 3; 20*)
7.93 - type unipoly = int list;
7.94 -
7.95 -"----------- auxiliary functions univariate -------------";
7.96 -"----------- auxiliary functions univariate -------------";
7.97 -"----------- auxiliary functions univariate -------------";
7.98 -
7.99 -(*subsection {* calculations for integers *}*)
7.100 -
7.101 - infix div'
7.102 - fun a div' b =
7.103 - if a < 0 then abs a div b * ~1
7.104 - else a div b;
7.105 -
7.106 -(*subsection {* modulo calculations for integers *}*)
7.107 - fun mod_inv r m =
7.108 - let
7.109 - fun modi (r, rold, m, rinv) =
7.110 - if rinv < m
7.111 - then
7.112 - if r mod m = 1
7.113 - then rinv
7.114 - else modi (rold * (rinv + 1), rold, m, rinv + 1)
7.115 - else 0
7.116 - in modi (r, r, m, 1) end;
7.117 -
7.118 - fun mod_div a b m =
7.119 - a * (mod_inv b m) mod m;
7.120 -
7.121 - fun aprox_root a =
7.122 - let fun root' a n =
7.123 - if n*n < a then root' a (n+1) else n
7.124 - in root' a 1
7.125 - end
7.126 -
7.127 - (* r = r1 mod m1 and r= r2 mod m2 *)
7.128 - fun CRA_2 (r1: int, m1: int, r2: int, m2: int ) =
7.129 - (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1
7.130 -
7.131 -(*subsection {* prime opertations *}*)
7.132 -
7.133 - fun is_prime primelist number =
7.134 - if length primelist >0
7.135 - then
7.136 - if (number mod (nth primelist 0))=0
7.137 - then false
7.138 - else is_prime (nth_drop 0 primelist) number
7.139 - else true
7.140 -
7.141 - fun primeList number =
7.142 - let
7.143 - fun make_primelist list last number =
7.144 - if (nth list (length list - 1)) < number
7.145 - then
7.146 - if ( is_prime list (last + 2))
7.147 - then make_primelist (list @ [(last + 2)]) (last + 2) number
7.148 - else make_primelist list (last + 2) number
7.149 - else list
7.150 - in
7.151 - if number < 3
7.152 - then [2]
7.153 - else make_primelist [2,3] 3 number end
7.154 -
7.155 - (*find a prime greater p not dividing the number a*)
7.156 - fun find_next_prime_not_divide a p =
7.157 - let
7.158 - val next = nth (primeList (p + 1)) (length (primeList (p + 1)) - 1) ;
7.159 - in
7.160 - if a mod next <> 0
7.161 - then next
7.162 - else find_next_prime_not_divide a next
7.163 - end
7.164 -
7.165 -(*subsection {* calculations for univariate polynomials *}*)
7.166 -
7.167 - fun lc (uvp: unipoly) =
7.168 - if nth uvp (length uvp- 1) <>0
7.169 - then nth uvp (length uvp- 1)
7.170 - else lc (nth_drop (length uvp- 1) uvp);
7.171 -
7.172 - fun deg (uvp: unipoly) =
7.173 - if nth uvp (length uvp- 1) <>0
7.174 - then length uvp - 1
7.175 - else deg (nth_drop (length uvp- 1) uvp)
7.176 -
7.177 - fun lc_of_unipoly_not_0 [] = [](* and delete lc=0*)
7.178 - | lc_of_unipoly_not_0 (uvp: unipoly) =
7.179 - if nth uvp (length uvp - 1) =0
7.180 - then lc_of_unipoly_not_0 (nth_drop (length uvp- 1) uvp)
7.181 - else uvp;
7.182 -
7.183 - fun normirt_polynom (poly1: unipoly) (m: int) =
7.184 - let
7.185 - val poly1 = lc_of_unipoly_not_0 poly1
7.186 - val lc_a=lc poly1;
7.187 - fun normirt poly1 b m lc_a i =
7.188 - if i=0
7.189 - then [mod_div (nth poly1 i) lc_a m]@b
7.190 - else normirt poly1 ( [mod_div(nth poly1 i) lc_a m]@b) m lc_a (i- 1) ;
7.191 - in
7.192 - if length(poly1)=0
7.193 - then poly1
7.194 - else normirt poly1 [] m lc_a (length(poly1)- 1)
7.195 - end
7.196 -
7.197 - infix %*
7.198 - fun (a: unipoly) %* b = map2 Integer.mult a (replicate (length (a)) b)
7.199 -
7.200 - infix %/
7.201 - fun (poly: unipoly) %/ b = (* =quotient*)
7.202 - let fun division poly b = poly div' b;
7.203 - in map2 division poly (replicate (length (poly)) b) end
7.204 -
7.205 - infix %-
7.206 - fun (poly: unipoly) %- b =
7.207 - let fun minus poly b = poly - b;
7.208 - in map2 minus poly (replicate (length (poly)) b) end
7.209 -
7.210 - infix %+%
7.211 - fun (a: unipoly) %+% (b: unipoly) = map2 Integer.add a b
7.212 -
7.213 - infix %-%
7.214 - fun (a: unipoly) %-% (b: unipoly) =
7.215 - let fun minus a b = a-b;
7.216 - in map2 minus a b end
7.217 -
7.218 - (* if poly2|poly1 *)
7.219 - infix %|%
7.220 - fun [b: int] %|% [a: int] =
7.221 - if abs b<= abs a andalso a mod b = 0
7.222 - then true
7.223 - else false
7.224 - | (poly2: unipoly) %|% (poly1: unipoly) =
7.225 - let
7.226 - val b = (replicate (length poly1 - length(lc_of_unipoly_not_0 poly2)) 0) @ lc_of_unipoly_not_0 poly2 ;
7.227 - val c = lc poly1 div' lc b;
7.228 - val rest = lc_of_unipoly_not_0 (poly1 %-% (b %* c));
7.229 - in
7.230 - if rest = []
7.231 - then true
7.232 - else
7.233 - if c<>0 andalso length rest >= length poly2
7.234 - then poly2 %|% rest
7.235 - else false
7.236 - end
7.237 -
7.238 - fun poly_centr (poly: unipoly) (m: int) =
7.239 - let
7.240 - fun midle m =
7.241 - let val mid = m div' 2
7.242 - in
7.243 - if m mod mid = 0
7.244 - then mid
7.245 - else mid+1
7.246 - end
7.247 - fun centr a m mid =
7.248 - if a > mid
7.249 - then a - m
7.250 - else a
7.251 - fun polyCentr poly poly' m mid counter =
7.252 - if length(poly) > counter
7.253 - then polyCentr poly (poly' @ [centr (nth poly counter) m mid]) m mid (counter+1)
7.254 - else (poly': unipoly)
7.255 - in polyCentr poly [] m (midle m) 0
7.256 - end
7.257 -
7.258 - fun sum_lmb (a: unipoly) exp =
7.259 - let fun sum' [a] _ = a
7.260 - | sum' (a: int list) exp =
7.261 - sum' ([((nth a 0)) + ( Integer.pow exp (nth a 1))] @ nth_drop 0 (nth_drop 0 a)) exp
7.262 - in sum' ([Integer.pow exp (nth a 0)] @ nth_drop 0 a) exp
7.263 - end;
7.264 -Integer.min ;
7.265 - fun landau_mignotte_bound a b =
7.266 - (Integer.pow (Integer.min (deg a) (deg b)) 2) * (abs (Integer.gcd (lc a) (lc b)))*
7.267 - (Int.min (abs((aprox_root (sum_lmb a 2)) div ~(lc a)), abs(( aprox_root (sum_lmb b 2)) div ~(lc b))));
7.268 -
7.269 -(*subsection {* modulo calculations for polynomials *}*)
7.270 -
7.271 - fun CRA_2_poly (ma, mb) (r1, r2) =
7.272 - let
7.273 - fun CRA_2' (ma, mb) (r1, r2) = CRA_2 (r1, ma,r2, mb)
7.274 - in map (CRA_2' (ma, mb)) (r1 ~~ r2) end
7.275 -
7.276 - infix poly_mod
7.277 - fun uvp poly_mod m =
7.278 - let fun poly_mod' uvp m n =
7.279 - if n < length (uvp)
7.280 - then poly_mod' ((nth_drop 0 uvp)@[(nth uvp 0) mod m]) m ( n + 1 )
7.281 - else uvp (*end of poly_mod'*)
7.282 - in
7.283 - poly_mod' uvp m 0 end
7.284 -
7.285 - fun mod_poly_gcd (upoly1: unipoly) (upoly2: unipoly) (m: int) =
7.286 - let
7.287 - val moda = upoly1 poly_mod m;
7.288 - val modb = (replicate (length upoly1 - length(lc_of_unipoly_not_0 (upoly2 poly_mod m))) 0)
7.289 - @ (lc_of_unipoly_not_0 (upoly2 poly_mod m)) ;
7.290 - val c = mod_div (lc moda) (lc modb) m
7.291 - val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m)
7.292 - in if rest = []
7.293 - then [upoly2, [0]]
7.294 - else
7.295 - if length rest < length upoly2
7.296 - then mod_poly_gcd upoly2 rest m
7.297 - else mod_poly_gcd rest upoly2 m
7.298 - end
7.299 -
7.300 - fun mod_gcd_p (poly1: unipoly) (poly2: unipoly) (p: int) =
7.301 - let val gcd_p = mod_poly_gcd poly1 poly2 p
7.302 - in normirt_polynom (nth gcd_p 0) p
7.303 - end
7.304 -
7.305 - fun gcd_n (up: unipoly) =
7.306 - if length up = 2
7.307 - then abs (Integer.gcd (nth up 0)(nth up 1))
7.308 - else gcd_n ([abs (Integer.gcd (nth up 0)(nth up 1))]@(nth_drop 0 (nth_drop 0 up)))
7.309 -
7.310 -fun gcd_n up = abs (Integer.gcds up);
7.311 -
7.312 - fun pp [_: int] = [1]
7.313 - | pp (poly1: unipoly) =
7.314 - if (poly1 poly_mod (gcd_n poly1)) = (replicate (length (poly1)) 0)
7.315 - then poly1 %/ (gcd_n poly1)
7.316 - else poly1;
7.317 -
7.318 -"=========== code for [1] p.93: univariate ==============";
7.319 -"=========== code for [1] p.93: univariate ==============";
7.320 -"=========== code for [1] p.93: univariate ==============";
7.321 -(*subsection {* GCD_MOD Algorgithmus, code for [1] p.93: univariate *}*)
7.322 -
7.323 - fun GCD_MOD (a: unipoly) (b: unipoly) =
7.324 - if a = [0] orelse b = [0] then [1] else
7.325 - if a = b then a else
7.326 - let
7.327 -(*1*)val d = abs (Integer.gcd (lc a) (lc b));
7.328 - val M = 2*d*landau_mignotte_bound a b;
7.329 - fun GOTO2 a b d M p = (*==============================*)
7.330 - let
7.331 -(*2*) val p = find_next_prime_not_divide d p
7.332 - val c_p = normirt_polynom ( mod_gcd_p a b p) p
7.333 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
7.334 - fun GOTO3 a b d M p g_p = (*~~~~~~~~~~~~~~~~~~~~~~*)
7.335 -(*3*) if (deg g_p) = 0
7.336 - then [1]
7.337 - else
7.338 - let
7.339 - val P = p
7.340 - val g = g_p
7.341 - fun WHILE a b d M P g p = (*------------------*)
7.342 -(*4*) if P > M
7.343 - then g
7.344 - else
7.345 - let
7.346 - val p = find_next_prime_not_divide d p
7.347 - val c_p = normirt_polynom ( mod_gcd_p a b p) p
7.348 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
7.349 - in
7.350 - if deg g_p < deg g
7.351 - then GOTO3 a b d M p g_p
7.352 - else
7.353 - if deg g_p = deg g
7.354 - then
7.355 - let
7.356 - val g = CRA_2_poly (P,p)(g,g_p)
7.357 - val P = P*p
7.358 - val g = poly_centr(g poly_mod P)P
7.359 - in WHILE a b d M P g p end
7.360 - else WHILE a b d M P g p
7.361 - end (*----------------------------------*)
7.362 -
7.363 -
7.364 - val g = WHILE a b d M P g p (* << 1.Mal -----*)
7.365 - in g end (*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
7.366 -
7.367 - val g = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)
7.368 -(*5*) val g = pp g
7.369 - in
7.370 - if (g %|% a) andalso (g %|% b)
7.371 - then g
7.372 - else GOTO2 a b d M p
7.373 - end (*==============================================*)
7.374 -
7.375 -
7.376 - val g = GOTO2 a b d M (*p*)1(* << 1. Mal =============*)
7.377 - in g end;
7.378 -
7.379 -"----------- auxiliary functions multivariate -----------";
7.380 -"----------- auxiliary functions multivariate -----------";
7.381 -"----------- auxiliary functions multivariate -----------";
7.382 -(*subsection {* GCD_MODm Algorgithmus, auxiliary functions multivariate *}*)
7.383 -
7.384 -type monom = (int * int list);
7.385 -type multipoly = monom list;
7.386 -
7.387 -(*subsection {* calculations for multivariate polynomials *}*)
7.388 -
7.389 - fun listgreater a b =
7.390 - let fun gr a b = a>=b;
7.391 - fun all [] = true
7.392 - | all list =
7.393 - if nth list 0 then all (nth_drop 0 list) else false
7.394 - in all (map2 gr a b)
7.395 - end
7.396 -
7.397 - fun get_coef (mvp: multipoly) (place: int) =
7.398 - let fun coef ((coef,_): monom) = coef;
7.399 - in coef (nth mvp place) end
7.400 -
7.401 - fun get_varlist (mvp: multipoly) (place: int) =
7.402 - let fun varlist ((_,var): monom) = var;
7.403 - in varlist (nth mvp place) end
7.404 -
7.405 - fun get_coeflist (poly: multipoly) =
7.406 - let
7.407 - fun get_coeflist' (poly: multipoly) list =
7.408 - if poly = []
7.409 - then
7.410 - list
7.411 - else
7.412 - get_coeflist' (nth_drop 0 poly) (list @ [get_coef poly 0])
7.413 - in
7.414 - get_coeflist' poly []
7.415 - end
7.416 -
7.417 - fun add_var_to_multipoly (mvp: multipoly) (order: int) =
7.418 - let fun add ([]: multipoly) (_: int) (new: multipoly) = new
7.419 - | add (mvp: multipoly) (order: int) (new: multipoly) =
7.420 - let val (first, last) = chop order (get_varlist mvp 0)
7.421 - in add (nth_drop 0 mvp) (order) (new @ [((get_coef mvp 0),(first @ [0] @ last))]) end
7.422 - in add mvp order [] end
7.423 -
7.424 - fun cero_multipoly 1 = [(0,[0])]
7.425 - | cero_multipoly (diferent_var: int) =
7.426 - add_var_to_multipoly (cero_multipoly (diferent_var- 1)) 0;
7.427 -
7.428 - fun short_mv (mvp: multipoly) =
7.429 - let fun short (mvp: multipoly) (new: multipoly) =
7.430 - if length mvp =1
7.431 - then if (get_coef mvp 0) = 0
7.432 - then if new = [] then cero_multipoly (length (get_varlist mvp 0)) else new
7.433 - else new @ mvp
7.434 - else if get_varlist mvp 0 = get_varlist mvp 1
7.435 - then short ( [(get_coef mvp 0 + get_coef mvp 1,get_varlist mvp 0)]
7.436 - @ nth_drop 0 (nth_drop 0 mvp)) new
7.437 - else if (get_coef mvp 0) = 0 then short (nth_drop 0 mvp) new
7.438 - else short (nth_drop 0 mvp) (new @ [nth mvp 0])
7.439 - in short mvp [] end
7.440 -
7.441 - (* if a is greater than b *)
7.442 - fun greater_var [a: int] [b: int] =
7.443 - a > b
7.444 - | greater_var (a: int list) (b: int list) =
7.445 - if (nth a (length a - 1))= (nth b (length b - 1))
7.446 - then greater_var (nth_drop (length a- 1) a) (nth_drop (length b- 1) b)
7.447 - else (nth a (length a - 1)) > (nth b (length b - 1))
7.448 -
7.449 - fun order_multipoly (a: multipoly)=
7.450 - let
7.451 - val ordered = [nth a 0];
7.452 - fun order_mp [] [] = cero_multipoly (length (get_varlist a 0))
7.453 - | order_mp [] (ordered: multipoly) = short_mv ordered
7.454 - | order_mp a (ordered: multipoly) =
7.455 - if greater_var (get_varlist a 0) (get_varlist ordered (length ordered - 1))
7.456 - then order_mp (nth_drop 0 a)(ordered @ [nth a 0])
7.457 - else let
7.458 - val rest = [nth ordered (length ordered - 1)];
7.459 - fun order_mp' [] (new: multipoly) (rest: multipoly) = new @ rest
7.460 - | order_mp' (ordered': multipoly) (new: multipoly) (rest: multipoly) =
7.461 - if greater_var (get_varlist new 0) (get_varlist ordered' (length ordered' - 1))
7.462 - then ordered' @ new @ rest
7.463 - else order_mp' (nth_drop (length ordered' - 1) ordered') new
7.464 - ([nth ordered' (length ordered' - 1)]@ rest)
7.465 - in order_mp (nth_drop 0 a)
7.466 - (order_mp' (nth_drop (length ordered - 1) ordered) [nth a 0] rest)
7.467 - end
7.468 - in order_mp (nth_drop 0 a) ordered
7.469 - end
7.470 -
7.471 - fun lc_multipoly (mpoly: multipoly) =
7.472 - get_coef (order_multipoly mpoly) (length (order_multipoly mpoly) - 1);
7.473 -
7.474 - (* greatest variablegroup *)
7.475 - fun deg_multipoly (mvp: multipoly) =
7.476 - get_varlist (order_multipoly mvp) (length (order_multipoly mvp) - 1)
7.477 -
7.478 - fun max_deg_var [m: monom] (x: int) = nth (get_varlist [m] 0) x |
7.479 - max_deg_var (mvp: multipoly) (x: int) =
7.480 - let
7.481 - fun max_monom (m1: monom) (m2: monom) (x: int)=
7.482 - if nth (get_varlist [m1] 0) x < nth (get_varlist [m2] 0) x then m2 else m1;
7.483 - in
7.484 - max_deg_var ([max_monom (nth(mvp) 0)( nth(mvp) 1) x] @ nth_drop 0 (nth_drop 0 mvp)) x
7.485 - end
7.486 -
7.487 - fun greater_deg (monom1: monom) (monom2: monom)=
7.488 - if (max_deg_var [monom1] (length(get_varlist [monom1] 0) - 1)) >
7.489 - (max_deg_var [monom2] (length (get_varlist [monom2] 0) - 1))
7.490 - then 1
7.491 - else if (max_deg_var [monom1] (length(get_varlist [monom1] 0) - 1)) <
7.492 - (max_deg_var [monom2] (length (get_varlist [monom2] 0) - 1))
7.493 - then 2
7.494 - else if length (get_varlist [monom1] 0) >1
7.495 - then
7.496 - greater_deg (1,(nth_drop 0 (get_varlist [monom1] 0))) (1,(nth_drop 0 (get_varlist [monom2] 0)))
7.497 - else 0
7.498 -
7.499 - infix %%/
7.500 - fun (poly: multipoly) %%/ b = (* =quotient*)
7.501 - let fun division monom b = (((get_coef [monom] 0) div' b),get_varlist [monom] 0);
7.502 - in order_multipoly(map2 division poly (replicate (length(poly)) b)) end
7.503 -
7.504 - fun cont [a: int] = a
7.505 - | cont (poly1: unipoly) = gcd_n poly1
7.506 - fun cont_multipoly [(a,_): monom] = a
7.507 - | cont_multipoly (poly: multipoly) = gcd_n (get_coeflist poly)
7.508 -
7.509 - fun evaluate_mv (mvp: multipoly) (var: int) (value: int) =
7.510 - let fun eval ([]: multipoly) (_: int) (_: int) (new: multipoly) = order_multipoly new |
7.511 - eval (mvp: multipoly) (var: int) (value: int) (new: multipoly) =
7.512 - eval (nth_drop 0 mvp) var value
7.513 - (new @ [((Integer.pow (nth (get_varlist mvp 0) var)value) * get_coef mvp 0,
7.514 - nth_drop var (get_varlist mvp 0))]);
7.515 - in eval mvp var value [] end
7.516 -
7.517 - fun multi_to_uni (mvp: multipoly) =
7.518 - if length (get_varlist mvp 0) = 1
7.519 - then let fun mtu ([]: multipoly) (uvp: unipoly) = uvp |
7.520 - mtu (mvp: multipoly) (uvp: unipoly) =
7.521 - if length uvp = (nth(get_varlist mvp 0) 0)
7.522 - then mtu (nth_drop 0 mvp) (uvp @ [get_coef mvp 0])
7.523 - else mtu mvp (uvp @ [0])
7.524 - in mtu (order_multipoly mvp) [] end
7.525 - else error "Polynom has more than one variable!";
7.526 -
7.527 - fun uni_to_multi (uvp: unipoly) =
7.528 - let fun utm ([]: unipoly) (mvp: multipoly) (_: int)= short_mv mvp
7.529 - | utm (uvp: unipoly) (mvp: multipoly) (counter: int) =
7.530 - utm (nth_drop 0 uvp) (mvp @ [((nth uvp 0),[counter])]) (counter+1)
7.531 - in utm uvp [] 0 end
7.532 -
7.533 - fun find_new_r (mvp1: multipoly) (mvp2: multipoly) (old: int) =
7.534 - let val poly1 = evaluate_mv mvp1 (length (get_varlist mvp1 0) - 2) (old + 1)
7.535 - val poly2 = evaluate_mv mvp2 (length (get_varlist mvp2 0) - 2) (old + 1);
7.536 - in
7.537 - if max_deg_var poly1 (length (get_varlist poly1 0) - 1)= max_deg_var mvp1 (length (get_varlist mvp1 0) - 1) orelse
7.538 - max_deg_var poly2 (length (get_varlist poly2 0) - 1)= max_deg_var mvp2 (length (get_varlist mvp2 0) - 1)
7.539 - then old + 1
7.540 - else find_new_r (mvp1) (mvp2) (old + 1)
7.541 - end
7.542 -
7.543 - fun mult_with_new_var ([]: multipoly) (_: unipoly) (_: int) = []
7.544 - | mult_with_new_var (mvp: multipoly) (uvp: unipoly) (order: int) =
7.545 - let fun mult ([]: multipoly) (_: unipoly) (_: int) (new: multipoly) (_: int) = short_mv new
7.546 - | mult (mvp: multipoly) (uvp: unipoly) (order: int) (new: multipoly) (_: int) =
7.547 - let fun mult' (_: multipoly) ([]: unipoly) (_: int) (new': multipoly) (_) = order_multipoly new'
7.548 - | mult' (mvp': multipoly) (uvp': unipoly) (order: int) (new': multipoly) (c2': int) =
7.549 - let val (first, last) = chop order (get_varlist mvp' 0)
7.550 - in mult' mvp' (nth_drop 0 uvp') order
7.551 - (new' @ [(((get_coef mvp' 0) * (nth uvp' 0)),(first @ [c2'] @ last))]) (c2'+1) end
7.552 - in mult (nth_drop 0 mvp) uvp order (new @ (mult' mvp uvp order [] 0)) (0) end
7.553 - in mult mvp uvp order [] 0 end
7.554 -
7.555 - fun greater_deg_multipoly (var1: int list) (var2: int list) =
7.556 - if var1 = [] andalso var2 =[] then 0
7.557 - else if (nth var1 (length var1 - 1)) = (nth var2 (length var1 - 1) )
7.558 - then greater_deg_multipoly (nth_drop (length var1 - 1) var1) (nth_drop (length var1 - 1) var2)
7.559 - else if (nth var1 (length var1 - 1)) > (nth var2 (length var1 - 1))
7.560 - then 1 else 2 ;
7.561 -
7.562 - infix %%+%%
7.563 - fun ([]: multipoly) %%+%% (mvp2: multipoly) = mvp2
7.564 - | (mvp1: multipoly) %%+%% ([]: multipoly) = mvp1
7.565 - | (mvp1: multipoly) %%+%% (mvp2: multipoly) =
7.566 - let fun plus ([]: multipoly) (mvp2: multipoly) (new: multipoly) = order_multipoly (new @ mvp2)
7.567 - | plus (mvp1: multipoly) ([]: multipoly) (new: multipoly) = order_multipoly (new @ mvp1)
7.568 - | plus (mvp1: multipoly) (mvp2: multipoly) (new: multipoly) =
7.569 - if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 0
7.570 - then plus (nth_drop 0 mvp1) (nth_drop 0 mvp2)
7.571 - (new @ [(((get_coef mvp1 0) + (get_coef mvp2 0)), get_varlist mvp1 0)])
7.572 - else if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 1
7.573 - then plus mvp1 (nth_drop 0 mvp2) (new @ [nth mvp2 0])
7.574 - else plus (nth_drop 0 mvp1) mvp2 (new @ [nth mvp1 0])
7.575 - in plus mvp1 mvp2 [] end
7.576 -
7.577 - infix %%-%%
7.578 - fun (mvp1: multipoly) %%-%% (mvp2: multipoly) =
7.579 - let fun neg ([]: multipoly) (new: multipoly) = order_multipoly new
7.580 - | neg (mvp: multipoly) (new: multipoly) =
7.581 - neg (nth_drop 0 mvp) (new @ [(((get_coef mvp 0) * ~1), get_varlist mvp 0)])
7.582 - val neg_mvp2 = neg mvp2
7.583 - in mvp1 %%+%% (neg_mvp2 []) end
7.584 -
7.585 - infix %%*%%
7.586 - fun (mvp1: multipoly) %%*%% (mvp2: multipoly) =
7.587 - let fun mult ([]: multipoly) (_: multipoly) (_: multipoly) (new: multipoly) = order_multipoly new
7.588 - | mult (mvp1: multipoly) ([]: multipoly) (regular_mvp2: multipoly) (new: multipoly) = mult (nth_drop 0 mvp1) regular_mvp2 regular_mvp2 new
7.589 - | mult (mvp1: multipoly) (mvp2: multipoly) (regular_mvp2: multipoly) (new: multipoly) =
7.590 - mult mvp1 (nth_drop 0 mvp2) regular_mvp2 (new @ [(((get_coef mvp1 0) * (get_coef mvp2 0)), ((get_varlist mvp1 0) %+% (get_varlist mvp2 0)))])
7.591 - in if (length mvp1) > (length mvp2) then mult mvp1 mvp2 mvp2 [] else mult mvp2 mvp1 mvp1 [] end
7.592 -
7.593 - (* if poly1|poly2 *)
7.594 - infix %%|%%
7.595 - fun [(coef2, var2): monom] %%|%% [(coef1, var1): monom] =
7.596 - ( (listgreater var1 var2)
7.597 - andalso (coef1 mod coef2) = 0)
7.598 - | (_: multipoly) %%|%% [(0,_)]= true
7.599 - | (poly2: multipoly) %%|%% (poly1: multipoly) =
7.600 - if [nth poly2 (length poly2 - 1)] %%|%% [nth poly1 (length poly1 - 1)]
7.601 - then poly2 %%|%% (poly1 %%-%%
7.602 - ([((get_coef poly1 (length poly1 - 1)) div (get_coef poly2 (length poly2 - 1)),
7.603 - (get_varlist poly1 (length poly1 - 1)) %-%(get_varlist poly2 (length poly2 - 1)))] %%*%%
7.604 - poly2))
7.605 - else false
7.606 -
7.607 -(*subsection {* Newtoninterpolation *}*)
7.608 - (* first step *)
7.609 - fun newton_first (x: int list) (f: multipoly list) (order: int) =
7.610 - let val polynom =(add_var_to_multipoly (nth f 0) order) %%+%%
7.611 - ((mult_with_new_var (((nth f 1)%%-%% (nth f 0))%%/
7.612 - (nth x 1) - (nth x 0))) [(nth x 0) * ~1, 1] order);
7.613 - val new_value_poly = [(nth x 0) * ~1, 1];
7.614 - val steps = [((nth f 1) %%-%%(nth f 0))%%/ ((nth x 1) - (nth x 0))]
7.615 - in (polynom, new_value_poly, steps) end
7.616 -
7.617 - fun newton_mv (x: int list) (f: multipoly list) (steps: multipoly list) (t: unipoly) (p: multipoly) (order: int) =
7.618 - if length x = 2
7.619 - then let val (polynom, new_value_poly, steps) = newton_first x [(nth f 0), (nth f 1)] order
7.620 - in (polynom, new_value_poly, steps) end
7.621 - else let val new_value_poly = multi_to_uni((uni_to_multi t) %%*%% (uni_to_multi [(nth x (length x - 2) )* ~1, 1]));
7.622 - val new_steps = [((nth f (length f - 1)) %%-%% (nth f (length f - 2))) %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))];
7.623 - fun next_step ([]: multipoly list) (new_steps: multipoly list) (_: int list) = new_steps
7.624 - | next_step (steps: multipoly list) (new_steps: multipoly list) (x': int list) =
7.625 - next_step (nth_drop 0 steps)
7.626 - (new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0)) %%/
7.627 - ((nth x' (length x' - 1)) - (nth x' (length x' - 3))))])
7.628 - ( nth_drop (length x' - 2) x')
7.629 - val steps = next_step steps new_steps x;
7.630 - val polynom' = (p %%+%% (mult_with_new_var (nth steps (length steps - 1)) new_value_poly order));
7.631 - in (order_multipoly(polynom'), new_value_poly, steps) end;
7.632 -
7.633 -"=========== code for [1] p.95: multivariate ============";
7.634 -"=========== code for [1] p.95: multivariate ============";
7.635 -"=========== code for [1] p.95: multivariate ============";
7.636 -(*subsection {* GCD_MODm algorithm, code for [1] p.95: multivariate *}*)
7.637 -
7.638 -fun GCD_MODm a b n s r=
7.639 - if greater_var (deg_multipoly b) (deg_multipoly a) then GCD_MODm b a n s r else
7.640 -(*0*) if s = 0
7.641 - then uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %*
7.642 - (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))
7.643 - else
7.644 - let
7.645 -(*1*) val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
7.646 -(*2*) fun GOTO2 a b n s M r_list steps = (*==============================*)
7.647 - let
7.648 - val r = find_new_r a b r;
7.649 - val r_list = r_list @ [r];
7.650 - val g_r = GCD_MODm (order_multipoly (evaluate_mv a (s- 1) r)) ( order_multipoly (evaluate_mv b (s- 1) r)) (n- 1) (s- 1) 0
7.651 - val steps = steps @ [g_r];
7.652 -(*3*) fun GOTO3 a b n s M g_r r r_list steps = (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)
7.653 - let
7.654 - val m = 1;
7.655 - fun WHILE a b n s M m r r_list steps g g_n mult = (* ----------------------- *)
7.656 - if m > M
7.657 - then
7.658 - if g_n %%|%% a andalso g_n %%|%% b
7.659 - then g_n
7.660 - else GOTO2 a b n s M r_list steps
7.661 - else
7.662 - let
7.663 - val r = find_new_r a b r;
7.664 - val r_list = r_list @ [r];
7.665 - val g_r = GCD_MODm (order_multipoly (evaluate_mv a (s- 1) r))
7.666 - (order_multipoly (evaluate_mv b (s- 1) r)) (n- 1) (s- 1) 0
7.667 - in
7.668 - if greater_var (deg_multipoly g) (deg_multipoly g_r)
7.669 - then GOTO3 a b n s M g_r r r_list steps
7.670 - else if (deg_multipoly g)= (deg_multipoly g_r)
7.671 - then
7.672 - let
7.673 - val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1)
7.674 - in if (nth steps (length steps - 1)) = (cero_multipoly s)
7.675 - then WHILE a b n s M (M+1) r r_list steps g_r g_n new
7.676 - else WHILE a b n s M (m+1) r r_list steps g_r g_n new
7.677 - end
7.678 - else WHILE a b n s M (m+1) r r_list steps g g_n mult
7.679 - end (* WHILE *)
7.680 - in (* GOTO3*) (* ----------------------- *)
7.681 - WHILE a b n s M m r r_list steps g_r ( cero_multipoly (s+1)) [1]
7.682 - end (*GOTO3*)
7.683 - in (* GOTO2*) (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)
7.684 - GOTO3 a b n s M g_r r r_list steps
7.685 - end (*GOTO2*)
7.686 - in (*==============================*)
7.687 - GOTO2 a b n s M [] []
7.688 - end (*end 0*);
7.689 -
7.690 -"--- begin of univariate polynomials --------------------";
7.691 -"--- begin of univariate polynomials --------------------";
7.692 -"--- begin of univariate polynomials --------------------";
7.693 -
7.694 -"----------- fun mod_inv --------------------------------";
7.695 -"----------- fun mod_inv --------------------------------";
7.696 -"----------- fun mod_inv--- -----------------------------";
7.697 -"~~~~~ fun mod_inv, args:"; val (r, m) = (5,7);
7.698 -"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (r, r, m, 1);
7.699 -rinv < m; (*=true*)
7.700 -r mod m = 1; (*=false*)
7.701 -"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
7.702 -rinv < m;(*=true*)
7.703 -r mod m = 1; (*=false*)
7.704 -"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
7.705 -rinv < m;(*=true*)
7.706 -r mod m = 1; (*=true*)
7.707 -"~~~~~ to mod_inv return val:"; val (rinv) = (rinv);
7.708 -
7.709 -if mod_inv 5 7 = 3 then () else error "mod_inv 5 7 = 3 changed";
7.710 -if mod_inv 3 7 = 5 then () else error "mod_inv 3 7 = 5 changed";
7.711 -if mod_inv 4 339 = 85 then () else error "mod_inv 4 339 = 85 changed";
7.712 -
7.713 -"----------- fun CRA_2 ----------------------------------";
7.714 -"----------- fun CRA_2 ----------------------------------";
7.715 -"----------- fun CRA_2 ----------------------------------";
7.716 -
7.717 -if CRA_2(17,9,3,4) = 35 then () else error "CRA_2(17,9,3,4)=35 changed";
7.718 -if CRA_2(7,2,6,11) = 17 then () else error "CRA_2(7,2,6,11) = 17 changed";
7.719 -
7.720 -"----------- fun lc -------------------------------------";
7.721 -"----------- fun lc -------------------------------------";
7.722 -"----------- fun lc -------------------------------------";
7.723 -if lc [3,4,5,6] = 6 then () else error "lc (3,4,5,6) = 6 changed" ;
7.724 -if lc [3,4,5,6,0] = 6 then () else error "lc (3,4,5,6,0) = 6 changed" ;
7.725 -
7.726 -"----------- fun deg ------------------------------------";
7.727 -"----------- fun deg ------------------------------------";
7.728 -"----------- fun deg ------------------------------------";
7.729 -if deg [3,4,5,6] = 3 then () else error "lc (3,4,5,6) = 3 changed" ;
7.730 -if deg [3,4,5,6,0] = 3 then () else error "lc (3,4,5,6) = 6 changed";
7.731 -
7.732 -"----------- fun primeList ------------------------------";
7.733 -"----------- fun primeList ------------------------------";
7.734 -"----------- fun primeList ------------------------------";
7.735 -if primeList 1 = [2] then () else error " primeList 1 = [2] changed";
7.736 -if primeList 3 = [2,3] then () else error " primeList 3 = [2,3] changed";
7.737 -if primeList 6 = [2,3,5,7] then () else error " primeList 6 = [2,3,5,7] changed";
7.738 -if primeList 7 = [2,3,5,7] then () else error " primeList 7 = [2,3,5,7] changed";
7.739 -
7.740 -"----------- fun find_next_prime_not_divide -----------------";
7.741 -"----------- fun find_next_prime_not_divide -----------------";
7.742 -"----------- fun find_next_prime_not_divide -----------------";
7.743 -if find_next_prime_not_divide 15 2 = 7 then () else error "findNextPrimeNotdivide 15 2 = 7 changed";
7.744 -if find_next_prime_not_divide 90 1 = 7 then () else error "findNextPrimeNotdivide 90 1 = 7 changed";
7.745 -
7.746 -"----------- fun poly_mod -------------------------------";
7.747 -"----------- fun poly_mod -------------------------------";
7.748 -"----------- fun poly_mod -------------------------------";
7.749 -if ([5,4,7,8,1] poly_mod 5) = [0, 4, 2, 3, 1] then ()
7.750 -else error "[5,4,7,8,1] poly_mod 5 = [0, 4, 2, 3, 1] changed" ;
7.751 -if ([5,4,~7,8,~1] poly_mod 5) = [0, 4, 3, 3, 4] then ()
7.752 -else error "[5,4,~7,8,~1] poly_mod 5 = [0, 4, 3, 3, 4] changed" ;
7.753 -
7.754 -"----------- fun %* ------------------------------";
7.755 -"----------- fun %* ------------------------------";
7.756 -"----------- fun %* ------------------------------";
7.757 -if ([5,4,7,8,1] %* 5) = [25, 20, 35, 40, 5] then ()
7.758 -else error "[5,4,7,8,1] %* 5 = [25, 20, 35, 40, 5] changed";
7.759 -if ([5,4,~7,8,~1] %* 5) = [25, 20, ~35, 40, ~5] then ()
7.760 -else error "[5,4,~7,8,~1] %* 5 = [25, 20, ~35, 40, ~5] changed";
7.761 -
7.762 -"----------- fun %/ -------------------------------";
7.763 -"----------- fun %/ -------------------------------";
7.764 -"----------- fun %/ -------------------------------";
7.765 -if ([4,3,2,5,6] %/ 3) =[1, 1, 0, 1, 2] then ()
7.766 -else error "%/ [4,3,2,5,6] 3 =[1, 1, 0, 1, 2] changed";
7.767 -if ([4,3,2,0] %/ 3) =[1, 1, 0, 0] then ()
7.768 -else error "%/ [4,3,2,0] 3 =[1, 1, 0, 0] changed";
7.769 -
7.770 -"----------- fun %-% ------------------------";
7.771 -"----------- fun %-% ------------------------";
7.772 -"----------- fun %-% ------------------------";
7.773 -if ([8,~7,0,1] %-% [~2,2,3,0])=[10,~9,~3,1] then ()
7.774 -else error "[8,~7,0,1] %-% [~2,2,3,0]=[10,~9,~3,1] changed";
7.775 -if ([8,7,6,5,4] %-% [2,2,3,1,1])=[6,5,3,4,3] then ()
7.776 -else error "[8,7,6,5,4] %-% [2,2,3,1,1]=[6,5,3,4,3] changed";
7.777 -
7.778 -"----------- fun %+% ------------------------------";
7.779 -"----------- fun %+% ------------------------------";
7.780 -"----------- fun %+% ------------------------------";
7.781 -if ([8,~7,0,1] %+% [~2,2,3,0])=[6,~5,3,1] then ()
7.782 -else error "[8,~7,0,1] %+% [~2,2,3,0]=[6,~5,3,1] changed";
7.783 -if ([8,7,6,5,4] %+% [2,2,3,1,1])=[10,9,9,6,5] then ()
7.784 -else error "[8,7,6,5,4] %+% [2,2,3,1,1]=[10,9,9,6,5] changed";
7.785 -
7.786 -"----------- fun CRA_2_poly -----------------------------";
7.787 -"----------- fun CRA_2_poly -----------------------------";
7.788 -"----------- fun CRA_2_poly -----------------------------";
7.789 -if (CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()
7.790 -else error "CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";
7.791 -
7.792 -"----------- fun mod_div --------------------------------";
7.793 -"----------- fun mod_div --------------------------------";
7.794 -"----------- fun mod_div --------------------------------";
7.795 -if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";
7.796 -if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";
7.797 -if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";
7.798 -
7.799 -"----------- fun lc_of_unipoly_not_0 --------------------";
7.800 -"----------- fun lc_of_unipoly_not_0 --------------------";
7.801 -"----------- fun lc_of_unipoly_not_0 --------------------";
7.802 -if lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] then ()
7.803 -else error "lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] changed";
7.804 -if lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0,1, 2, 3, 4, 5] then ()
7.805 -else error "lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0, 1, 2, 3, 4, 5] changed";
7.806 -
7.807 -"----------- fun %|% -------------------------------";
7.808 -"----------- fun %|% -------------------------------";
7.809 -"----------- fun %|% -------------------------------";
7.810 -"~~~~~ fun %|% , args:"; val (poly1, poly2) = ([9,12,4],[3,2]);
7.811 -val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
7.812 -val c = lc poly1 div lc b;
7.813 -val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
7.814 -rest = [];(*=false*)
7.815 -c<>0 andalso length rest >= length poly2; (*=true*)
7.816 -"~~~~~ fun %|% , args:"; val (poly1, poly2) = (rest, poly2);
7.817 -val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
7.818 -val c = lc poly1 div lc b;
7.819 -val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));
7.820 -rest = [];(*=true*)
7.821 -"~~~~~ to return val:"; val (divide) = (true);
7.822 -
7.823 -if ([4] %|% [6]) = false then () else error "[4] %|% [6] = false changed";
7.824 -if ( [8] %|%[16,0]) = true then () else error "[8] %|%[16,0] = true changed";
7.825 -if ([3,2] %|%[0,0,9,12,4] ) = true then ()
7.826 -else error "[3,2] %|%[0,0,9,12,4] = true changed";
7.827 -if ([8,0] %|% [16]) = true then () else error "[8,0] %|% [16] = true changed";
7.828 -
7.829 -"----------- fun mod_poly_gcd ------------------------";
7.830 -"----------- fun mod_poly_gcd ------------------------";
7.831 -"----------- fun mod_poly_gcd ------------------------";
7.832 -
7.833 -if ( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7)=[[2, 6, 0, 2, 6], [0]]
7.834 -then () else error "( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [ [5, 1, 0, 5, 1], [0]] changed";
7.835 -if ( mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7)=[[2, 6, 0, 2, 6], [0]] then ()
7.836 -else error "mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7=[[1, 3, 0, 1, 3], [0]] changed";
7.837 -if mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] then ()
7.838 -else error " mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] changed";
7.839 -
7.840 -"~~~~~ fun mod_poly_gcd , args:";
7.841 -val (a,b,m) = ([ ~13, 2,12],[ ~14, 2],5);
7.842 -val moda = a poly_mod m;
7.843 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
7.844 -val c = mod_div (lc moda) (lc modb) m;
7.845 -val rest = lc_of_unipoly_not_0 (moda %-% (modb %* c) poly_mod m) poly_mod m;
7.846 -rest = []; (*=false*)
7.847 -length rest < length b; (*=false*)
7.848 -"~~~~~ fun mod_poly_gcd , args:";
7.849 -val (a,b,m,d) = (rest, b , m , [c]);
7.850 -val moda = a poly_mod m
7.851 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
7.852 -val c = mod_div (lc moda) (lc modb) m
7.853 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
7.854 -rest = [];(*=flase*)
7.855 -length rest < length b; (*=true*)
7.856 -"~~~~~ fun mod_poly_gcd , args:";
7.857 -val (a,b,m,d) = (b, rest, m, [c] @ d);
7.858 -val moda = a poly_mod m
7.859 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
7.860 -val c = mod_div (lc moda) (lc modb) m
7.861 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
7.862 -rest = [];(*=flase*)
7.863 -length rest < length b; (*=false*)
7.864 -"~~~~~ fun mod_poly_gcd , args:";
7.865 -val (a,b,m,d) = (b, rest, m, [c] @ d);
7.866 -val moda = a poly_mod m
7.867 -val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;
7.868 -val c = mod_div (lc moda) (lc modb) m
7.869 -val rest = lc_of_unipoly_not_0 ((moda %-% (modb %* c)) poly_mod m);
7.870 -rest = [];(*=true*)
7.871 -"~~~~~ to return val:"; val ([b, rest, lsit])=[b, [0], [c] @ d];
7.872 -
7.873 -
7.874 -"----------- fun normirt_polynom ------------------------";
7.875 -"----------- fun normirt_polynom ------------------------";
7.876 -"----------- fun normirt_polynom ------------------------";
7.877 -if (normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] then ()
7.878 -else error "(normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] changed"
7.879 -
7.880 -"----------- fun poly_centr -----------------------------";
7.881 -"----------- fun poly_centr -----------------------------";
7.882 -"----------- fun poly_centr -----------------------------";
7.883 -if (poly_centr [7,3,5,8,1,3] 10) = [~3, 3, 5, ~2, 1, 3] then ()
7.884 -else error "poly_centr [7,3,5,8,1,3] 10 = [~3, 3, 5, ~2, 1, 3] cahnged";
7.885 -
7.886 -"----------- fun mod_gcd_p -------------------------------";
7.887 -"----------- fun mod_gcd_p -------------------------------";
7.888 -"----------- fun mod_gcd_p -------------------------------";
7.889 -if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5) = [1, 1, 1, 1] then ()
7.890 -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";
7.891 -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 ()
7.892 -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";
7.893 -
7.894 -
7.895 -"----------- fun gcd_n ----------------------------";
7.896 -"----------- fun gcd_n ----------------------------";
7.897 -"----------- fun gcd_n ----------------------------";
7.898 -if gcd_n [3,9,12,21] = 3 then () else error "gcd_n [3,9,12,21] = 3 changed";
7.899 -if gcd_n [12,16,32,44] = 4 then () else error "gcd_n [12,16,32,44] = 4 changed";
7.900 -if gcd_n [4,5,12] = 1 then () else error "gcd_n [4,5,12] = 1 changed";
7.901 -
7.902 -"----------- fun pp -------------------------------------";
7.903 -"----------- fun pp -------------------------------------";
7.904 -"----------- fun pp -------------------------------------";
7.905 -if pp [12,16,32,44] = [3, 4, 8, 11] then () else error "pp [12,16,32,44] = [3, 4, 8, 11] changed";
7.906 -if pp [4,5,12] = [4, 5, 12] then () else error "pp [4,5,12] = [4, 5, 12] changed";
7.907 -
7.908 -"----------- fun GCD_MOD --------------------------------";
7.909 -"----------- fun GCD_MOD --------------------------------";
7.910 -"----------- fun GCD_MOD --------------------------------";
7.911 -
7.912 -"~~~~~ fun GCD_MOD a b , args:"; val (a, b) = ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2]);
7.913 -val d = abs (Integer.gcd (lc a) (lc b));
7.914 -val M =2*d* landau_mignotte_bound a b;
7.915 -(*val g = GOTO2 a b d (*p*)1 M*)
7.916 - "~~~~~ fun GOTO2 a b d M p , args:"; val (a, b, d, p, M) = (a,b,d,3,M);
7.917 - val p = find_next_prime_not_divide d p
7.918 - val c_p = normirt_polynom ( mod_gcd_p a b p) p
7.919 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.920 - (*val (p, g) = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)*)
7.921 - "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
7.922 - (deg g_p) = 0; (*=false*)
7.923 - val P = p
7.924 - val g = g_p;
7.925 - (*(p, g) = WHILE a b d M P g (* << 1.Mal -----*)*)
7.926 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,p,g);
7.927 - P > M ;(*false*)
7.928 - val p = find_next_prime_not_divide d p
7.929 - val c_p = normirt_polynom ( mod_gcd_p a b p) p
7.930 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.931 - deg g_p < deg g; (*=fasle*)
7.932 - deg g_p = deg g; (*false*);
7.933 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.934 - P > M ;(*false*)
7.935 - val p = find_next_prime_not_divide d p
7.936 - val c_p = normirt_polynom ( mod_gcd_p a b p) p
7.937 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.938 - deg g_p < deg g; (*=false*)
7.939 - deg g_p = deg g; (*=true*)
7.940 - val g = CRA_2_poly (P,p)(g,g_p)
7.941 - val P = P*p;
7.942 - val g = poly_centr(g poly_mod P)P;
7.943 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.944 - P > M ;(*false*)
7.945 - val p = find_next_prime_not_divide d p
7.946 - val c_p = normirt_polynom (mod_gcd_p a b p) p
7.947 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.948 - deg g_p < deg g; (*=true*)
7.949 - "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
7.950 - (deg g_p) = 0; (*false*)
7.951 - val P = p
7.952 - val g = g_p;
7.953 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.954 - P > M ;(*false*)
7.955 - val p = find_next_prime_not_divide d p
7.956 - val c_p = normirt_polynom (mod_gcd_p a b p) p
7.957 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.958 - deg g_p < deg g; (*=fasle*)
7.959 - deg g_p = deg g;(*true*)
7.960 - val g = CRA_2_poly (P,p)(g,g_p)
7.961 - val P = P*p;
7.962 - val g = poly_centr(g poly_mod P)P;
7.963 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.964 - P > M ;(*false*)
7.965 - val p = find_next_prime_not_divide d p
7.966 - val c_p = normirt_polynom (mod_gcd_p a b p) p
7.967 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.968 - deg g_p < deg g; (*=false*)
7.969 - deg g_p = deg g;(*true*)
7.970 - val g = CRA_2_poly (P,p)(g,g_p)
7.971 - val P = P*p;
7.972 - val g = poly_centr(g poly_mod P)P;
7.973 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.974 - P > M ;(*false*)
7.975 - val p = find_next_prime_not_divide d p
7.976 - val c_p = normirt_polynom (mod_gcd_p a b p) p
7.977 - val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
7.978 - deg g_p < deg g; (*=false*)
7.979 - deg g_p = deg g;(*true*)
7.980 - val g = CRA_2_poly (P,p)(g,g_p)
7.981 - val P = P*p;
7.982 - val g = poly_centr(g poly_mod P)P;
7.983 - "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
7.984 - P > M ;(*true*)
7.985 - "~~~~~ to return WHILE val:"; val (g,p) = (g,p);
7.986 - "~~~~~ to return GOTO3 val:"; val (g,p) = (g,p);
7.987 - val g = pp g;
7.988 - (g %|% a) andalso (g %|% b);(*=true*)
7.989 -"~~~~~ to return GOTO2 val:"; val (g) = (g);
7.990 -"~~~~~ to return GCD_MOD val:"; val (g) = (g);
7.991 -
7.992 -"----------- fun GCD_MOD---------------------------------";
7.993 -"----------- fun GCD_MOD --------------------------------";
7.994 -"----------- fun GCD_MOD --------------------------------";
7.995 -if GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] then ()
7.996 -else error "GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] changed";
7.997 -
7.998 -"------------------------------------------------------------";
7.999 -"--------------------- Multivariate Case --------------------";
7.1000 -"------------------------------------------------------------";
7.1001 -
7.1002 -"----------- fun get_coef --------------------------- ---";
7.1003 -"----------- fun get_coef -------------------------------";
7.1004 -"----------- fun get_coef -------------------------------";
7.1005 -if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 then () else
7.1006 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 changed";
7.1007 -if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 then () else
7.1008 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 changed";
7.1009 -if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 then () else
7.1010 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 changed";
7.1011 -if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 then () else
7.1012 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 changed";
7.1013 -if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 then () else
7.1014 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 changed";
7.1015 -
7.1016 -
7.1017 -"----------- fun get_varlist ----------------------------";
7.1018 -"----------- fun get_varlist ----------------------------";
7.1019 -"----------- fun get_varlist ----------------------------";
7.1020 -if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] then () else
7.1021 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] changed";
7.1022 -if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] then () else
7.1023 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] changed";
7.1024 -if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] then () else
7.1025 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] changed";
7.1026 -if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] then () else
7.1027 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] changed";
7.1028 -if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] then () else
7.1029 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] changed";
7.1030 -
7.1031 -"----------- fun get_coeflist ---------------------------";
7.1032 -"----------- fun get_coeflist ---------------------------";
7.1033 -"----------- fun get_coeflist ---------------------------";
7.1034 -if get_coeflist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] = [1,2,3,4,5] then () else
7.1035 -error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] = [1,2,3,4,5] changed";
7.1036 -
7.1037 -"----------- fun add_var_to_multipoly -------------------";
7.1038 -"----------- fun add_var_to_multipoly -------------------";
7.1039 -"----------- fun add_var_to_multipoly -------------------";
7.1040 -if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] then () else
7.1041 -error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] changed";
7.1042 -if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3])] then () else
7.1043 -error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0, 2]), (2, [2, 0, 3]) changed";
7.1044 -if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] then () else
7.1045 -error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] changed";
7.1046 -
7.1047 -"----------- fun cero_multipoly -------------------------";
7.1048 -"----------- fun cero_multipoly -------------------------";
7.1049 -"----------- fun cero_multipoly -------------------------";
7.1050 -if cero_multipoly 1 = [(0, [0])] then () else
7.1051 -error "cero_multipoly 1 = [(0, [0])] changed";
7.1052 -if cero_multipoly 5 = [(0, [0, 0, 0, 0, 0])] then () else
7.1053 -error "cero_multipoly 1 = [(0, [0, 0, 0, 0, 0])] changed";
7.1054 -
7.1055 -"----------- fun short_mv -------------------------------";
7.1056 -"----------- fun short_mv -------------------------------";
7.1057 -"----------- fun short_mv -------------------------------";
7.1058 -"~~~~~ fun short_mv, args:"; val (mvp) = ([(~12, [1]), (2, [1]), (4, [0])]);
7.1059 -"~~~~~ fun short , args:"; val (mvp, new) = (mvp,([]:multipoly));
7.1060 -length mvp =1(*false*);
7.1061 - get_varlist mvp 0 = get_varlist mvp 1; (* true*)
7.1062 -"~~~~~ fun short , args:"; val (mvp, new) =
7.1063 -([(get_coef mvp 0 + get_coef mvp 1,get_varlist mvp 0)] @ nth_drop 0 (nth_drop 0 mvp), new );
7.1064 -length mvp =1(*false*);
7.1065 -get_varlist mvp 0 = get_varlist mvp 1;(*false*)
7.1066 -"~~~~~ fun short , args:"; val (mvp, new) = ((nth_drop 0 mvp),new @ [nth mvp 0]);
7.1067 -length mvp =1(*true*);
7.1068 -"~~~~~ to return val:"; val (shortform) = (new @ mvp);
7.1069 -
7.1070 -if short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])]
7.1071 -then () else error "short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])] changed"
7.1072 -
7.1073 -"----------- fun greater_var ----------------------------";
7.1074 -"----------- fun greater_var ----------------------------";
7.1075 -"----------- fun greater_var ----------------------------";
7.1076 -if greater_var [3] [4]
7.1077 -then error " greater_var [3] [4] changed = false" else ();
7.1078 -if greater_var [1,2] [0,3]
7.1079 -then error " greater_var [1,2] [0,3] changed = false" else ();
7.1080 -if greater_var [1,2] [3,0]
7.1081 -then () else error " greater_var [1,2] [3,0] = true changed";
7.1082 -if greater_var [1,2] [1,2]
7.1083 -then error " greater_var [1,2] [1,2] changed = false" else ();
7.1084 -
7.1085 -"----------- fun order_multipoly ------------------------";
7.1086 -"----------- fun order_multipoly ------------------------";
7.1087 -"----------- fun order_multipoly ------------------------";
7.1088 -if order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])]
7.1089 -then () else error "order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])] changed"
7.1090 -
7.1091 -"----------- fun lc_multipoly ---------------------------";
7.1092 -"----------- fun lc_multipoly ---------------------------";
7.1093 -"----------- fun lc_multipoly ---------------------------";
7.1094 -if lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1
7.1095 -then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1 changed";
7.1096 -if lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2
7.1097 -then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2 changed";
7.1098 -
7.1099 -"----------- fun deg_multipoly -------------------------";
7.1100 -"----------- fun deg_multipoly -------------------------";
7.1101 -"----------- fun deg_multipoly -------------------------";
7.1102 -if deg_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [1,2]
7.1103 -then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1 changed";
7.1104 -if deg_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = [1,2]
7.1105 -then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2 changed";
7.1106 -
7.1107 -
7.1108 -
7.1109 -"----------- fun max_deg_var ----------------------------";
7.1110 -"----------- fun max_deg_var ----------------------------";
7.1111 -"----------- fun max_deg_var ----------------------------";
7.1112 -if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 then ()
7.1113 -else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 changed";
7.1114 -if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 then ()
7.1115 -else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 changed";
7.1116 -
7.1117 -"----------- infix %%/ ----------------------------------";
7.1118 -"----------- infix %%/ ----------------------------------";
7.1119 -"----------- infix %%/ ----------------------------------";
7.1120 -if ([(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2) = [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] then ()
7.1121 -else error "[(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2 = [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] changed";
7.1122 -
7.1123 -"----------- fun cont -----------------------------------";
7.1124 -"----------- fun cont -----------------------------------";
7.1125 -"----------- fun cont -----------------------------------";
7.1126 -if cont [4,8,12,~2,0,~4] = 2 then () else error " cont [4,8,12,~2,0,~4] = 2 changed";
7.1127 -if cont [3,6,9,19] = 1 then () else error " cont [3,6,9,19] = 1 changed";
7.1128 -
7.1129 -"----------- fun cont_multipoly -------------------------";
7.1130 -"----------- fun cont_multipoly -------------------------";
7.1131 -"----------- fun cont_multipoly -------------------------";
7.1132 -if cont_multipoly [(3,[1,2,3,1])] = 3 then () else error " cont_multipoly [(3,[1,2,3,1])] = 3 changed";
7.1133 -if cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 then ()
7.1134 -else error "cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 changed";
7.1135 -if cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 then ()
7.1136 -else error "cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 changed";
7.1137 -
7.1138 -"----------- fun evaluate_mv ----------------------------";
7.1139 -"----------- fun evaluate_mv ----------------------------";
7.1140 -"----------- fun evaluate_mv ----------------------------";
7.1141 -if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])]
7.1142 -then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])] changed";
7.1143 -if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 = [(2, [0]), (~2, [2])]
7.1144 -then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 =[(2, [0]), (~2, [2])] changed";
7.1145 -if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])]
7.1146 -then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])] changed";
7.1147 -if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])]
7.1148 -then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])] changed";
7.1149 -
7.1150 -"----------- fun mutli_to_uni ---------------------------";
7.1151 -"----------- fun mutli_to_uni ---------------------------";
7.1152 -"----------- fun mutli_to_uni ---------------------------";
7.1153 -if multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1]
7.1154 -then () else error " multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1] changed";
7.1155 -
7.1156 -"~~~~~ fun multi_to_uni , args:"; val (mvp) = ([(~3, [0]), (1, [0]), (3, [1]), (~6, [1]),(2,[3])]);
7.1157 -val mvp = order_multipoly mvp;
7.1158 -"~~~~~ fun short, args:"; val (mvp, uvp) = ((short_mv mvp), ([]:unipoly));
7.1159 -length uvp = (nth(get_varlist mvp 0) 0);(*true*)
7.1160 -"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
7.1161 -length uvp = (nth(get_varlist mvp 0) 0);(*true*)
7.1162 -"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
7.1163 -length uvp = (nth(get_varlist mvp 0) 0);(*false*)
7.1164 -"~~~~~ fun short, args:"; val (mvp, uvp) = (mvp, (uvp @ [0]));
7.1165 -length uvp = (nth(get_varlist mvp 0) 0);(*true*)
7.1166 -"~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));
7.1167 -"~~~~~ to return val:"; val (unipoly) = (uvp);
7.1168 -
7.1169 -
7.1170 -
7.1171 -"----------- fun find_new_r ----------------------------";
7.1172 -"----------- fun find_new_r ----------------------------";
7.1173 -"----------- fun find_new_r ----------------------------";
7.1174 -if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 then ()
7.1175 -else error " find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 changed";
7.1176 -if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 then ()
7.1177 -else error "find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 changed";
7.1178 -if find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 then ()
7.1179 -else error "find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 changed";
7.1180 -
7.1181 -"----------- fun newton_mv -----------------------------";
7.1182 -"----------- fun newton_mv -----------------------------";
7.1183 -"----------- fun newton_mv -----------------------------";
7.1184 -if uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] then ()
7.1185 -else error "uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] changed";
7.1186 -if uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] then ()
7.1187 -else error "uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] changed";
7.1188 -
7.1189 -"----------- fun mult_with_new_var ---------------------";
7.1190 -"----------- fun mult_with_new_var ---------------------";
7.1191 -"----------- fun mult_with_new_var ---------------------";
7.1192 -if mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] then ()
7.1193 -else error "mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] changed";
7.1194 -
7.1195 -"----------- fun greater_deg_multipoly -----------------";
7.1196 -"----------- fun greater_deg_multipoly -----------------";
7.1197 -"----------- fun greater_deg_multipoly -----------------";
7.1198 -greater_deg_multipoly [1,2,3] [1,2,4] =2;
7.1199 -
7.1200 -"----------- infix %%+%% -------------------------------";
7.1201 -"----------- infix %%+%% -------------------------------";
7.1202 -"----------- infix %%+%% -------------------------------";
7.1203 -if ([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])]
7.1204 -then () else error "([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] changed";
7.1205 -
7.1206 -"----------- infix %%-%% -------------------------------";
7.1207 -"----------- infix %%-%% -------------------------------";
7.1208 -"----------- infix %%-%% -------------------------------";
7.1209 -if ([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])]
7.1210 -then () else error "([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])] changed";
7.1211 -
7.1212 -"----------- infix %%*%% -------------------------------";
7.1213 -"----------- infix %%*%% -------------------------------";
7.1214 -"----------- infix %%*%% -------------------------------";
7.1215 -if ([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])]
7.1216 -then () else error "([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])] changed";
7.1217 -if ([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] then ()
7.1218 -else error "([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] changed";
7.1219 -if ([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0
7.1220 -then () else error
7.1221 -"([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 changed";
7.1222 -
7.1223 -"----------- fun newton_first --------------------------";
7.1224 -"----------- fun newton_first --------------------------";
7.1225 -"----------- fun newton_first --------------------------";
7.1226 -if newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 =
7.1227 -([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])
7.1228 -then () else error
7.1229 -"newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 = "
7.1230 -"([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";
7.1231 -
7.1232 -"----------- fun newton_mv -----------------------------";
7.1233 -"----------- fun newton_mv -----------------------------";
7.1234 -"----------- fun newton_mv -----------------------------";
7.1235 -if newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 =
7.1236 -([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])
7.1237 -then () else error
7.1238 -"newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 = "
7.1239 -"([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";
7.1240 -if newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 =
7.1241 -([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]])
7.1242 -then () else error
7.1243 -"newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 ="
7.1244 -" ([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]]) changed";
7.1245 -
7.1246 -
7.1247 -"~~~~~ fun newton_mv, args:"; val (x,f,steps,t,p,order) = ([1, 2, 3, 4],[[(9, [0]), (5, [1])], [(16, [0]), (7, [1])]], [[(5, [0]), (2, [1])], [(1, [0])]], [2, ~3, 1], [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])], 0 );
7.1248 -length x = 2; (* false *)
7.1249 -val new_value_poly = multi_to_uni((uni_to_multi t) %%*%% (uni_to_multi [(nth x (length x - 2) )* ~1, 1]));
7.1250 -val new_steps = [((nth f (length f - 1)) %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))) %%-%% ((nth f (length f - 2)))];
7.1251 -
7.1252 -"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (steps, new_steps, x);
7.1253 -steps = []; (*false*)
7.1254 -
7.1255 -"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0))) %%/
7.1256 - ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' - 2) x');steps = []; (*false*)
7.1257 -
7.1258 -"~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps - 1)) %%-%%(nth steps 0))) %%/
7.1259 - ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' - 2) x');
7.1260 -steps = []; (*true*)
7.1261 -val steps = new_steps;
7.1262 -val polynom' = p %%+%% (mult_with_new_var (nth steps (length steps - 1)) new_value_poly order);
7.1263 -
7.1264 -"----------- fun listgreater ---------------------------";
7.1265 -"----------- fun listgreater ---------------------------";
7.1266 -"----------- fun listgreater ---------------------------";
7.1267 -if listgreater [1,2,3,4,5] [1,2,3,4,5] = true then ()
7.1268 -else error " listgreater [1,2,3,4,5] [1,2,3,4,5] = true changed";
7.1269 -if listgreater [1,2,3,4] [1,2,3,5] = false then ()
7.1270 -else error "listgreater [1,2,3,4] [1,2,3,5] = false changed" ;
7.1271 -if listgreater [1,4,5,4] [0,3,4,5] = false then ()
7.1272 -else error "listgreater [1,2,3,4] [0,3,4,5] = false changed ";
7.1273 -
7.1274 -"----------- fun greater_deg_multipoly -----------------";
7.1275 -"----------- fun greater_deg_multipoly -----------------";
7.1276 -"----------- fun greater_deg_multipoly -----------------";
7.1277 -if greater_deg_multipoly [1,2,3,4,5] [1,2,3,4,5] = 0 then ()
7.1278 -else error " greater_deg_multipoly [1,2,3,4,5] [1,2,3,4,5] = 0 changed";
7.1279 -if greater_deg_multipoly [1,2,3,4] [5,2,8,1] = 1 then ()
7.1280 -else error "greater_deg_multipoly [1,2,3,4] [1,2,3,5] = 1 changed" ;
7.1281 -if greater_deg_multipoly [1,4,5,4] [0,3,4,5] = 2 then ()
7.1282 -else error "greater_deg_multipoly [1,2,3,4] [0,3,4,5] = 2 changed ";
7.1283 -
7.1284 -"----------- finfix %%|%% ------------------------------";
7.1285 -"----------- finfix %%|%% ------------------------------";
7.1286 -"----------- finfix %%|%% ------------------------------";
7.1287 -if [(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] then ()
7.1288 -else error "[(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] = true changed";
7.1289 -if [(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] then ()
7.1290 -else error "[(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = true changed";
7.1291 -if [(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])]
7.1292 -then error "[(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = false changed" else ();
7.1293 -
7.1294 -"----------- fun GCD_MODm ------------------------------";
7.1295 -"----------- fun GCD_MODm ------------------------------";
7.1296 -"----------- fun GCD_MODm ------------------------------";
7.1297 -
7.1298 -if GCD_MODm
7.1299 - [(~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])]
7.1300 - [(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 1 0
7.1301 - = [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] then () else error
7.1302 -"GCD_MODm [(~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])]"
7.1303 -"[(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 1 0 "
7.1304 -"= [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] changed";
7.1305 -(* -xy +xy^2z+yz - 1*)(* xy +1*) (*=*) (*xy - 1*)
7.1306 -if GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0
7.1307 -= [(1, [0, 0, 0]), (1, [1, 1, 0])] then () else error
7.1308 -"GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0 "
7.1309 -"= [(1, [0, 0, 0]), (1, [1, 1, 0])] changed";
7.1310 -
7.1311 -" ~~~ 1 ~~~";
7.1312 -"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ([(1,[1,2,1])], [(1,[1,2,1])], 3, 2, 0);
7.1313 -s = 0; (*false*)
7.1314 -val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
7.1315 - " ~~~ 1 ~~~";
7.1316 -"~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);
7.1317 -val r = find_new_r a b r;
7.1318 -val r_list = r_list @ [r];
7.1319 -val (a',b',n',s',r',r_list',steps',M') = ( a,b,n,s,r,r_list,steps,M);
7.1320 -(*g_1=*)
7.1321 - " ~~~ 1.1 ~~~";
7.1322 -"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1323 -s = 0; (*false*)
7.1324 -
7.1325 -val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
7.1326 - " ~~~ 1.1 ~~~";
7.1327 -"~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);
7.1328 -val r = find_new_r a b r;
7.1329 -val r_list = r_list @ [r];
7.1330 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1331 -(*g_2=*)
7.1332 - " ~~~ 1.1.1 ~~~";
7.1333 -"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1334 -s = 0; (*true*)
7.1335 -val g_2 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))
7.1336 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1337 -val g_r = g_2;
7.1338 -val steps = steps @ [g_r];
7.1339 - " ~~~ 1.1 ~~~";
7.1340 -"~~~~~ fun GOTO3 , args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
7.1341 -(a, b, n, s, M, g_r, r, r_list, steps);
7.1342 -val m = 1;
7.1343 - " ~~~ 1.1.W1 ~~~";
7.1344 -"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
7.1345 -(a, b, n, s, M, m,r, r_list, steps, g_r, ( cero_multipoly (s+1)), [1]);
7.1346 -; (*false*)
7.1347 -val r = find_new_r a b r;
7.1348 -val r_list = r_list @ [r];
7.1349 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1350 -(*g_3=*)
7.1351 - " ~~~ 1.1.W1.1 ~~~";
7.1352 -"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1353 -s = 0; (*true*)
7.1354 -val g_3 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1355 - " ~~~ 1.1.W1 ~~~";
7.1356 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1357 -val g_r = g_3;
7.1358 -greater_var (deg_multipoly g) (deg_multipoly g_r) (*false*);
7.1359 -(deg_multipoly g)= (deg_multipoly g_r) (*true*);
7.1360 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1361 -if g_n = [(1,[1,1])] then () else error "g_n is false" ;
7.1362 -(nth steps (length steps - 1)) = (cero_multipoly s) (* false*);
7.1363 - " ~~~ 1.1.W2 ~~~";
7.1364 -"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
7.1365 -(a, b, n, s, M, (m+1), r, r_list, steps, g_r, g_n, new);
7.1366 -m > M; (*false*)
7.1367 -val r = find_new_r a b r;
7.1368 -val r_list = r_list @ [r];
7.1369 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1370 -(*g_4=*)
7.1371 -" ~~~ 1.1.W2.1 ~~~";
7.1372 -"~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s- 1) r)), ( order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1373 -s = 0; (*true*)
7.1374 -val g_r =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1375 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1376 -" ~~~ 1.1.W2 ~~~";
7.1377 -greater_var (deg_multipoly g) (deg_multipoly g_r) (*false*);
7.1378 -(deg_multipoly g)= (deg_multipoly g_r) (*true*);
7.1379 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1380 -if g_n = [(1,[1,1])] then () else error "g_n is false" ;
7.1381 -(nth steps (length steps - 1)) = (cero_multipoly s) (* true*);
7.1382 - " ~~~ 1.1.W3 ~~~";
7.1383 -"~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =
7.1384 -(a, b, n, s, M, (M+1), r, r_list, steps, g_r, g_n, new);
7.1385 -m > M; (*true*)
7.1386 -g_n %%|%% a andalso g_n %%|%% b(*true*);
7.1387 - " ~~~ 1 ~~~";
7.1388 -val ( a,b,n,s,r,r_list,steps,M) = (a',b',n',s',r',r_list',steps',M');
7.1389 -val g_1=g_n;
7.1390 -val g_r = g_1;
7.1391 -val steps = steps @ [g_r];
7.1392 - " ~~~ 1 ~~~";
7.1393 -"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
7.1394 -(a, b, n, s, M, g_r, r, r_list, steps);
7.1395 -val m = 1;
7.1396 - " ~~~ 1.W1 ~~~";
7.1397 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1398 -(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
7.1399 -val g' = g;
7.1400 -m > M; (*false*)
7.1401 -val r = find_new_r a b r;
7.1402 -val r_list = r_list @ [r];
7.1403 -val (a',b',n',s',r',r_list',steps',m') = ( a,b,n,s,r,r_list,steps,m);
7.1404 -(* g_2 = GCD_MODm *)
7.1405 - " ~~~ 1.W1.1 ~~~";
7.1406 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1407 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1408 -s = 0; (*false*)
7.1409 -val M = 1 + Int.min(max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
7.1410 - " ~~~ 1.W1.1 ~~~";
7.1411 -"~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =
7.1412 -(a, b, n, s, M, [], []);
7.1413 -val r = find_new_r a b r;
7.1414 -val r_list = r_list @ [r];
7.1415 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1416 - " ~~~ 1.W1.1.1 ~~~";
7.1417 - (*g_r=*)
7.1418 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1419 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1420 -s = 0; (*true*)
7.1421 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1422 - " ~~~ 1.W1.1 ~~~";
7.1423 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1424 -val steps = steps @ [g_r];
7.1425 - " ~~~ 1.W1.1 ~~~";
7.1426 -"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
7.1427 -(a, b, n, s, M, g_r, r, r_list, steps);
7.1428 -val m = 1;
7.1429 -val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);
7.1430 - " ~~~ 1.W1.1.W1 ~~~";
7.1431 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1432 -(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
7.1433 -
7.1434 -m > M; (*false*)
7.1435 -val r = find_new_r a b r;
7.1436 -val r_list = r_list @ [r];
7.1437 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1438 -(* g_r = GCD_MODm *)
7.1439 - " ~~~ 1.W1.1.W1.1 ~~~";
7.1440 - (*g_r=*)
7.1441 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1442 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1443 -s = 0; (*true*)
7.1444 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1445 - " ~~~ 1.W1.1.W1 ~~~";
7.1446 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1447 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1448 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1449 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1450 -(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
7.1451 - " ~~~ 1.W1.1.W2 ~~~";
7.1452 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1453 -(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
7.1454 -m > M; (*false*)
7.1455 -val r = find_new_r a b r;
7.1456 -val r_list = r_list @ [r];
7.1457 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1458 -(* g_r = GCD_MODm *)
7.1459 - " ~~~ 1.W1.1.W2.1 ~~~";
7.1460 - (*g_r=*)
7.1461 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1462 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1463 -s = 0; (*true*)
7.1464 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1465 - " ~~~ 1.W1.1.W2 ~~~";
7.1466 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1467 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1468 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1469 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1470 -(nth steps (length steps - 1)) = (cero_multipoly s); (*true*)
7.1471 - " ~~~ 1.W1.1.W3 ~~~";
7.1472 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1473 -(a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);
7.1474 -m>M; (*true*)
7.1475 -g_n %%|%% a andalso g_n %%|%% b; (*true*)
7.1476 - " ~~~ 1.W1.1 ~~~";
7.1477 -val ( a,b,n,s,r,r_list,g,M) = (a',b',n',s',r',r_list',g',M');
7.1478 -val g_r = g_n;
7.1479 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1480 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1481 -val (g_n,steps,m,new) = (g_n'',steps',m'',new'');
7.1482 - " ~~~ 1.W2 ~~~";
7.1483 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps new g_n (s- 1);
7.1484 -(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
7.1485 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1486 -(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
7.1487 -m > M; (*false*)
7.1488 -val r = find_new_r a b r;
7.1489 -val r_list = r_list @ [r];
7.1490 -val (a',b',n',s',r',r_list',steps',m',M',g') = ( a,b,n,s,r,r_list,steps,m,M,g);
7.1491 -(* g_2 = GCD_MODm *)
7.1492 - " ~~~ 1.W2.1 ~~~";
7.1493 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1494 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1495 -s = 0; (*false*)
7.1496 -val M = 1 + Int.min(max_deg_var a (length(get_varlist a 0)- 2), max_deg_var b (length(get_varlist b 0)- 2));
7.1497 - " ~~~ 1.W2.1 ~~~";
7.1498 -"~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =
7.1499 -(a, b, n, s, M, [], []);
7.1500 -val r = find_new_r a b r;
7.1501 -val r_list = r_list @ [r];
7.1502 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1503 - " ~~~ 1.W2.1.1 ~~~";
7.1504 - (*g_r=*)
7.1505 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1506 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1507 -s = 0; (*true*)
7.1508 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1509 - " ~~~ 1.W2.1 ~~~";
7.1510 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1511 -val steps = steps @ [g_r];
7.1512 - " ~~~ 1.W2.1 ~~~";
7.1513 -"~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =
7.1514 -(a, b, n, s, M, g_r, r, r_list, steps);
7.1515 -val m = 1;
7.1516 -val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);
7.1517 - " ~~~ 1.W2.1.W1 ~~~";
7.1518 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1519 -(a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);
7.1520 -m > M; (*false*)
7.1521 -val r = find_new_r a b r;
7.1522 -val r_list = r_list @ [r];
7.1523 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1524 -(* g_r = GCD_MODm *)
7.1525 - " ~~~ 1.W2.1.W1.1 ~~~";
7.1526 - (*g_r=*)
7.1527 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1528 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1529 -s = 0; (*true*)
7.1530 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1531 - " ~~~ 1.W2.1.W1 ~~~";
7.1532 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1533 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1534 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1535 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1536 -(nth steps (length steps - 1)) = (cero_multipoly s); (*false*)
7.1537 - " ~~~ 1.W2.1.W2 ~~~";
7.1538 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1539 -(a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);
7.1540 -m > M; (*false*)
7.1541 -val r = find_new_r a b r;
7.1542 -val r_list = r_list @ [r];
7.1543 -val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);
7.1544 -(* g_r = GCD_MODm *)
7.1545 - " ~~~ 1.W2.1.W2.1 ~~~";
7.1546 - (*g_r=*)
7.1547 -"~~~~~ fun GCD_MODm, args:"; val (a,b,n,s,r) =
7.1548 -((order_multipoly (evaluate_mv a (s- 1) r)),(order_multipoly (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);
7.1549 -s = 0; (*true*)
7.1550 -val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));
7.1551 - " ~~~ 1.W2.1.W2 ~~~";
7.1552 -val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');
7.1553 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1554 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1555 -val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1556 -(nth steps (length steps - 1)) = (cero_multipoly s); (*true*)
7.1557 - " ~~~ 1.W2.1.W3 ~~~";
7.1558 -"~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =
7.1559 -(a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);
7.1560 -m>M; (*true*)
7.1561 -g_n %%|%% a andalso g_n %%|%% b; (*true*)
7.1562 - " ~~~ 1.W2.1 ~~~";
7.1563 -val ( a,b,n,s,r,r_list,g) = (a',b',n',s',r',r_list',g');
7.1564 -val g_r = g_n;
7.1565 -greater_var (deg_multipoly g) (deg_multipoly g_r); (*false*)
7.1566 -(deg_multipoly g)= (deg_multipoly g_r); (*true*)
7.1567 -val (g_n,steps,m,mult) = (g_n'',steps',m'',new'');
7.1568 - " ~~~ 1.W2 ~~~";
7.1569 - val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s- 1);
7.1570 -"~~~~~ to return val:"; val (g) = (g_n);
7.1571 -
7.1572 -
7.1573 -
7.1574 -if g = [(1, [1, 2, 1])] then () else error "GCD_MODm [(1, [1, 2, 1])] [(1, [1, 2, 1])] has changes.";
7.1575 -
7.1576 -" ========================== END ========================== ";
7.1577 - fun nth _ [] = error "nth _ []" (*Isabelle2002, still saved the hours of update*)
7.1578 - | nth 1 (x::_) = x
7.1579 - | nth n (_::xs) = nth (n- 1) xs;
7.1580 -(*fun nth xs i = List.nth (xs, i); recent Isabelle: TODO update all isac code *)
7.1581 -
8.1 --- a/test/Tools/isac/Knowledge/rational-1.sml Fri Aug 06 12:09:06 2021 +0200
8.2 +++ b/test/Tools/isac/Knowledge/rational-1.sml Fri Aug 06 18:27:05 2021 +0200
8.3 @@ -119,8 +119,6 @@
8.4 "-------- fun cancel_p special cases -----------------------------------------------------------";
8.5 "-------- fun cancel_p special cases -----------------------------------------------------------";
8.6 val thy = @{theory Isac_Knowledge};
8.7 -(* cp from GCD_Poly_ML.thy; is overwritten by previous tests *)
8.8 - fun order (p: poly) = (add_monoms o (sort lex_ord')) p;
8.9
8.10 (*------- standard case: *)
8.11 val t = TermC.str2term "2 / 3 + 1 / 6 ::real";
8.12 @@ -204,7 +202,6 @@
8.13 (*+*)val [(2, [])] = b'; (* 2 * _2_ = 12 *)
8.14 (*+*)val [(12, [])] = c; (* 12 / 24 \<and> 12 / 12 ..gcd *)
8.15
8.16 -(*/------------------ TOODOO broken with "repair cancellation with zero polynomial" -----------\* )
8.17 val nomin = term_of_poly baseT expT vs
8.18 (((the (poly_of_term vs n1)) %%*%% b') %%+%% ((the (poly_of_term vs n2)) %%*%% a'));
8.19
8.20 @@ -228,7 +225,6 @@
8.21
8.22 val t' = HOLogic.mk_binop \<^const_name>\<open>divide\<close> (nomin, denom);
8.23 (*+*)val "0 / 24" = UnparseC.term t'
8.24 -( *\------------------ TOODOO broken with "repair cancellation with zero polynomial" -----------/*)
8.25
8.26 (*---------- fun cancel_p with Const AA *)
8.27 val thy = @{theory Partial_Fractions};
9.1 --- a/test/Tools/isac/Test_Isac.thy Fri Aug 06 12:09:06 2021 +0200
9.2 +++ b/test/Tools/isac/Test_Isac.thy Fri Aug 06 18:27:05 2021 +0200
9.3 @@ -308,7 +308,6 @@
9.4 ML_file "Knowledge/simplify.sml"
9.5 ML_file "Knowledge/poly.sml"
9.6 ML_file "Knowledge/gcd_poly_ml.sml"
9.7 - ML_file "Knowledge/gcd_poly_winkler.sml" (*must be after gcd_poly_ml.sml: redefines functions*)
9.8 ML_file "Knowledge/rational.sml" (*Test_Isac_Short*)
9.9 ML_file "Knowledge/equation.sml"
9.10 ML_file "Knowledge/root.sml"
10.1 --- a/test/Tools/isac/Test_Isac_Short.thy Fri Aug 06 12:09:06 2021 +0200
10.2 +++ b/test/Tools/isac/Test_Isac_Short.thy Fri Aug 06 18:27:05 2021 +0200
10.3 @@ -281,7 +281,6 @@
10.4 ML_file "Knowledge/poly-1.sml"
10.5 (*ML_file "Knowledge/poly-2.sml" Test_Isac_Short*)
10.6 ML_file "Knowledge/gcd_poly_ml.sml"
10.7 - ML_file "Knowledge/gcd_poly_winkler.sml" (*must be after gcd_poly_ml.sml: redefines functions*)
10.8 ML_file "Knowledge/rational-1.sml"
10.9 (*ML_file "Knowledge/rational-2.sml" Test_Isac_Short*)
10.10 ML_file "Knowledge/equation.sml"