1.1 --- a/src/Tools/isac/Knowledge/Rational.thy Mon Sep 16 11:28:43 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/Rational.thy Mon Sep 16 12:20:00 2013 +0200
1.3 @@ -1,77 +1,125 @@
1.4 -(* rationals, i.e. fractions of multivariate polynomials over the real field
1.5 +(* rationals, fractions of multivariate polynomials over the real field
1.6 author: isac team
1.7 - Copyright (c) isac team 2002
1.8 + Copyright (c) isac team 2002, 2013
1.9 Use is subject to license terms.
1.10
1.11 depends on Poly (and not on Atools), because
1.12 fractions with _normalized_ polynomials are canceled, added, etc.
1.13 -
1.14 - ATTENTION WN130616: WITH Unsynchronized.ref Rational.thy TAKES ~1MIN FOR EVALUATION
1.15 *)
1.16
1.17 theory Rational
1.18 imports Poly "~~/src/Tools/isac/Knowledge/GCD_Poly_ML"
1.19 begin
1.20
1.21 +section {* Constants for evaluation by "Calc" *}
1.22 consts
1.23
1.24 is'_expanded :: "real => bool" ("_ is'_expanded") (*RL->Poly.thy*)
1.25 is'_ratpolyexp :: "real => bool" ("_ is'_ratpolyexp")
1.26 get_denominator :: "real => real"
1.27 get_numerator :: "real => real"
1.28 -
1.29
1.30 -axioms(*axiomatization where*) (*.not contained in Isabelle2002,
1.31 - stated as axioms, TODO?: prove as theorems*)
1.32 +ML {*
1.33 +(*.the expression contains + - * ^ / only ?.*)
1.34 +fun is_ratpolyexp (Free _) = true
1.35 + | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ Free _ $ Free _) = true
1.36 + | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ Free _ $ Free _) = true
1.37 + | is_ratpolyexp (Const ("Groups.times_class.times",_) $ Free _ $ Free _) = true
1.38 + | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true
1.39 + | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ Free _ $ Free _) = true
1.40 + | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ t1 $ t2) =
1.41 + ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.42 + | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) =
1.43 + ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.44 + | is_ratpolyexp (Const ("Groups.times_class.times",_) $ t1 $ t2) =
1.45 + ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.46 + | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =
1.47 + ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.48 + | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ t1 $ t2) =
1.49 + ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.50 + | is_ratpolyexp _ = false;
1.51
1.52 - mult_cross: "[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)" (*and*)
1.53 - mult_cross1: " b ~= 0 ==> (a / b = c ) = (a = b * c)" (*and*)
1.54 - mult_cross2: " d ~= 0 ==> (a = c / d) = (a * d = c)" (*and*)
1.55 +(*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)
1.56 +fun eval_is_ratpolyexp (thmid:string) _
1.57 + (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =
1.58 + if is_ratpolyexp arg
1.59 + then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",
1.60 + Trueprop $ (mk_equality (t, @{term True})))
1.61 + else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",
1.62 + Trueprop $ (mk_equality (t, @{term False})))
1.63 + | eval_is_ratpolyexp _ _ _ _ = NONE;
1.64 +
1.65 +(*("get_denominator", ("Rational.get_denominator", eval_get_denominator ""))*)
1.66 +fun eval_get_denominator (thmid:string) _
1.67 + (t as Const ("Rational.get_denominator", _) $
1.68 + (Const ("Fields.inverse_class.divide", _) $ num $
1.69 + denom)) thy =
1.70 + SOME (mk_thmid thmid "" (term_to_string''' thy denom) "",
1.71 + Trueprop $ (mk_equality (t, denom)))
1.72 + | eval_get_denominator _ _ _ _ = NONE;
1.73 +
1.74 +(*("get_numerator", ("Rational.get_numerator", eval_get_numerator ""))*)
1.75 +fun eval_get_numerator (thmid:string) _
1.76 + (t as Const ("Rational.get_numerator", _) $
1.77 + (Const ("Fields.inverse_class.divide", _) $num
1.78 + $denom )) thy =
1.79 + SOME (mk_thmid thmid "" (term_to_string''' thy num) "",
1.80 + Trueprop $ (mk_equality (t, num)))
1.81 + | eval_get_numerator _ _ _ _ = NONE;
1.82 +*}
1.83 +
1.84 +section {* Theorems for rewriting *}
1.85 +
1.86 +axiomatization (* naming due to Isabelle2002, but not contained in Isabelle2002;
1.87 + many thms are due to RL and can be removed with updating the equation solver;
1.88 + TODO: replace by equivalent thms in recent Isabelle201x *)
1.89 +where
1.90 + mult_cross: "[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)" and
1.91 + mult_cross1: " b ~= 0 ==> (a / b = c ) = (a = b * c)" and
1.92 + mult_cross2: " d ~= 0 ==> (a = c / d) = (a * d = c)" and
1.93
1.94 - add_minus: "a + b - b = a"(*RL->Poly.thy*) (*and*)
1.95 - add_minus1: "a - b + b = a"(*RL->Poly.thy*) (*and*)
1.96 + add_minus: "a + b - b = a"(*RL->Poly.thy*) and
1.97 + add_minus1: "a - b + b = a"(*RL->Poly.thy*) and
1.98
1.99 - rat_mult: "a / b * (c / d) = a * c / (b * d)"(*?Isa02*) (*and*)
1.100 - rat_mult2: "a / b * c = a * c / b "(*?Isa02*) (*and*)
1.101 + rat_mult: "a / b * (c / d) = a * c / (b * d)"(*?Isa02*) and
1.102 + rat_mult2: "a / b * c = a * c / b "(*?Isa02*) and
1.103
1.104 - rat_mult_poly_l: "c is_polyexp ==> c * (a / b) = c * a / b" (*and*)
1.105 - rat_mult_poly_r: "c is_polyexp ==> (a / b) * c = a * c / b" (*and*)
1.106 + rat_mult_poly_l: "c is_polyexp ==> c * (a / b) = c * a / b" and
1.107 + rat_mult_poly_r: "c is_polyexp ==> (a / b) * c = a * c / b" and
1.108
1.109 (*real_times_divide1_eq .. Isa02*)
1.110 - real_times_divide_1_eq: "-1 * (c / d) =-1 * c / d " (*and*)
1.111 - real_times_divide_num: "a is_const ==>
1.112 - a * (c / d) = a * c / d " (*and*)
1.113 + real_times_divide_1_eq: "-1 * (c / d) = -1 * c / d " and
1.114 + real_times_divide_num: "a is_const ==> a * (c / d) = a * c / d " and
1.115
1.116 - real_mult_div_cancel2: "k ~= 0 ==> m * k / (n * k) = m / n" (*and*)
1.117 + real_mult_div_cancel2: "k ~= 0 ==> m * k / (n * k) = m / n" and
1.118 (*real_mult_div_cancel1: "k ~= 0 ==> k * m / (k * n) = m / n"..Isa02*)
1.119
1.120 - real_divide_divide1: "y ~= 0 ==> (u / v) / (y / z) = (u / v) * (z / y)" (*and*)
1.121 - real_divide_divide1_mg: "y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)" (*and*)
1.122 + real_divide_divide1: "y ~= 0 ==> (u / v) / (y / z) = (u / v) * (z / y)" and
1.123 + real_divide_divide1_mg: "y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)" and
1.124 (*real_divide_divide2_eq: "x / y / z = x / (y * z)"..Isa02*)
1.125
1.126 - rat_power: "(a / b)^^^n = (a^^^n) / (b^^^n)" (*and*)
1.127 -
1.128 + rat_power: "(a / b)^^^n = (a^^^n) / (b^^^n)" and
1.129
1.130 rat_add: "[| a is_const; b is_const; c is_const; d is_const |] ==>
1.131 - a / c + b / d = (a * d + b * c) / (c * d)" (*and*)
1.132 + a / c + b / d = (a * d + b * c) / (c * d)" and
1.133 rat_add_assoc: "[| a is_const; b is_const; c is_const; d is_const |] ==>
1.134 - a / c +(b / d + e) = (a * d + b * c)/(d * c) + e" (*and*)
1.135 + a / c +(b / d + e) = (a * d + b * c)/(d * c) + e" and
1.136 rat_add1: "[| a is_const; b is_const; c is_const |] ==>
1.137 - a / c + b / c = (a + b) / c" (*and*)
1.138 + a / c + b / c = (a + b) / c" and
1.139 rat_add1_assoc: "[| a is_const; b is_const; c is_const |] ==>
1.140 - a / c + (b / c + e) = (a + b) / c + e" (*and*)
1.141 + a / c + (b / c + e) = (a + b) / c + e" and
1.142 rat_add2: "[| a is_const; b is_const; c is_const |] ==>
1.143 - a / c + b = (a + b * c) / c" (*and*)
1.144 + a / c + b = (a + b * c) / c" and
1.145 rat_add2_assoc: "[| a is_const; b is_const; c is_const |] ==>
1.146 - a / c + (b + e) = (a + b * c) / c + e" (*and*)
1.147 + a / c + (b + e) = (a + b * c) / c + e" and
1.148 rat_add3: "[| a is_const; b is_const; c is_const |] ==>
1.149 - a + b / c = (a * c + b) / c" (*and*)
1.150 + a + b / c = (a * c + b) / c" and
1.151 rat_add3_assoc: "[| a is_const; b is_const; c is_const |] ==>
1.152 a + (b / c + e) = (a * c + b) / c + e"
1.153
1.154 section {* Cancellation and addition of fractions *}
1.155 -subsection {* Auxiliary functions and data *}
1.156 -subsubsection {* Conversion term <--> poly *}
1.157 +subsection {* Conversion term <--> poly *}
1.158 +subsubsection {* Convert a term to the internal representation of a multivariate polynomial *}
1.159 ML {*
1.160 fun monom_of_term vs (c, es) (Free (id, _)) =
1.161 if is_numeral id
1.162 @@ -118,9 +166,12 @@
1.163 in
1.164 case poly_of_term vs t of SOME _ => true | NONE => false
1.165 end
1.166 -val is_expanded = is_poly
1.167 +val is_expanded = is_poly (* TODO: check names *)
1.168 +val is_polynomial = is_poly (* TODO: check names *)
1.169 +*}
1.170
1.171 -(* convert internal representation of a multivariate polynomial to a term*)
1.172 +subsubsection {* Convert internal representation of a multivariate polynomial to a term *}
1.173 +ML {*
1.174 fun term_of_es _ _ _ [] = [] (*assumes same length for vs and es*)
1.175 | term_of_es baseT expT (_ :: vs) (0 :: es) =
1.176 [] @ term_of_es baseT expT vs es
1.177 @@ -147,149 +198,8 @@
1.178 in foldl (HOLogic.mk_binop "Groups.plus_class.plus") (hd monos, tl monos) end
1.179 *}
1.180
1.181 -text {*calculate in rationals: gcd, lcm, etc.
1.182 - (c) Stefan Karnel 2002
1.183 - Institute for Mathematics D and Institute for Software Technology,
1.184 - TU-Graz SS 2002 *}
1.185 -
1.186 -text {* Remark on notions in the documentation below:
1.187 - referring to the remark on 'polynomials' in Poly.sml we use
1.188 - [2] 'polynomial' normalform (Polynom)
1.189 - [3] 'expanded_term' normalform (Ausmultiplizierter Term),
1.190 - where normalform [2] is a special case of [3], i.e. [3] implies [2].
1.191 - Instead of
1.192 - 'fraction with numerator and nominator both in normalform [2]'
1.193 - 'fraction with numerator and nominator both in normalform [3]'
1.194 - we say:
1.195 - 'fraction in normalform [2]'
1.196 - 'fraction in normalform [3]'
1.197 - or
1.198 - 'fraction [2]'
1.199 - 'fraction [3]'.
1.200 - a 'simple fraction' is a term with '/' as outmost operator and
1.201 - numerator and nominator in normalform [2] or [3].
1.202 -*}
1.203 -
1.204 -ML {*
1.205 -val thy = @{theory};
1.206 -
1.207 -signature RATIONALI =
1.208 -sig
1.209 - type mv_monom
1.210 - type mv_poly
1.211 - val add_fraction_p_ : theory -> term -> (term * term list) option
1.212 - val calculate_Rational : rls
1.213 - val calc_rat_erls:rls
1.214 - val cancel_p : rls
1.215 - val cancel_p_ : theory -> term -> (term * term list) option
1.216 - val check_fraction : term -> (term * term) option
1.217 - val check_frac_sum : term -> ((term * term) * (term * term)) option
1.218 - val common_nominator_p : rls
1.219 - val common_nominator_p_ : theory -> term -> (term * term list) option
1.220 - val eval_is_expanded : string -> 'a -> term -> theory ->
1.221 - (string * term) option
1.222 - val expanded2polynomial : term -> term option
1.223 - val factout_p_ : theory -> term -> (term * term list) option
1.224 - val is_expanded : term -> bool
1.225 - val is_polynomial : term -> bool
1.226 -
1.227 - val mv_gcd : (int * int list) list -> mv_poly -> mv_poly
1.228 - val mv_lcm : mv_poly -> mv_poly -> mv_poly
1.229 -
1.230 - val norm_expanded_rat_ : theory -> term -> (term * term list) option
1.231 -(*WN0602.2.6.pull into struct !!!
1.232 - val norm_Rational : rls(*.normalizes an arbitrary rational term without
1.233 - roots into a simple and canceled fraction
1.234 - with normalform [2].*)
1.235 -*)
1.236 -(*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME
1.237 - rls (*.normalizes an rational term [2] without
1.238 - roots into a simple and canceled fraction
1.239 - with normalform [2].*)
1.240 -*)
1.241 - val norm_rational_ : theory -> term -> (term * term list) option
1.242 - val polynomial2expanded : term -> term option
1.243 - val rational_erls :
1.244 - rls (*.evaluates an arbitrary rational term with numerals.*)
1.245 -
1.246 -(*WN0210???SK: fehlen Funktionen, die exportiert werden sollen ? *)
1.247 -end (* sig *)
1.248 -
1.249 -(*.**************************************************************************
1.250 -survey on the functions
1.251 -~~~~~~~~~~~~~~~~~~~~~~~
1.252 - [2] 'polynomial' :rls | [3]'expanded_term':rls
1.253 ---------------------:------------------+-------------------:-----------------
1.254 - factout_p_ : | factout_ :
1.255 - cancel_p_ : | cancel_ :
1.256 - :cancel_p | :cancel
1.257 ---------------------:------------------+-------------------:-----------------
1.258 - common_nominator_p_: | common_nominator_ :
1.259 - :common_nominator_p| :common_nominator
1.260 - add_fraction_p_ : | add_fraction_ :
1.261 ---------------------:------------------+-------------------:-----------------
1.262 -???SK :norm_rational_p | :norm_rational
1.263 -
1.264 -This survey shows only the principal functions for reuse, and the identifiers
1.265 -of the rls exported. The list below shows some more useful functions.
1.266 -
1.267 -
1.268 -conversion from Isabelle-term to internal representation
1.269 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1.270 -
1.271 -... BITTE FORTSETZEN ...
1.272 -
1.273 -polynomial2expanded = ...
1.274 -expanded2polynomial = ...
1.275 -
1.276 -remark: polynomial2expanded o expanded2polynomial = I,
1.277 - where 'o' is function chaining, and 'I' is identity WN0210???SK
1.278 -
1.279 -functions for greatest common divisor and canceling
1.280 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1.281 -################################################################################
1.282 -## search Isabelle2009-2/src/HOL/Multivariate_Analysis
1.283 -## Amine Chaieb, Robert Himmelmann, John Harrison.
1.284 -################################################################################
1.285 -mv_gcd
1.286 -factout_
1.287 -factout_p_
1.288 -cancel_
1.289 -cancel_p_
1.290 -
1.291 -functions for least common multiple and addition of fractions
1.292 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1.293 -mv_lcm
1.294 -common_nominator_
1.295 -common_nominator_p_
1.296 -add_fraction_ (*.add 2 or more fractions.*)
1.297 -add_fraction_p_ (*.add 2 or more fractions.*)
1.298 -
1.299 -functions for normalform of rationals
1.300 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1.301 -WN0210???SK interne Funktionen f"ur norm_rational:
1.302 - schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?
1.303 -
1.304 -norm_rational_
1.305 -norm_expanded_rat_
1.306 -
1.307 -**************************************************************************.*)
1.308 -
1.309 -
1.310 -(*##*)
1.311 -structure RationalI : RATIONALI =
1.312 -struct
1.313 -(*##*)
1.314 -
1.315 -infix mem ins union; (*WN100819 updating to Isabelle2009-2*)
1.316 -fun x mem [] = false
1.317 - | x mem (y :: ys) = x = y orelse x mem ys;
1.318 -
1.319 -
1.320 -
1.321 -val is_expanded = is_poly
1.322 -val is_polynomial = is_poly
1.323 -
1.324 +subsection {* Apply gcd_poly for cancelling and adding fractions as terms *}
1.325 +ML {*
1.326 fun mk_noteq_0 baseT t =
1.327 Const ("HOL.Not", HOLogic.boolT --> HOLogic.boolT) $
1.328 (Const ("HOL.eq", [baseT, baseT] ---> HOLogic.boolT) $ t $ Free ("0", HOLogic.realT))
1.329 @@ -297,7 +207,10 @@
1.330 fun mk_asms baseT ts =
1.331 let val as' = filter_out is_num ts (* asm like "2 ~= 0" is needless *)
1.332 in map (mk_noteq_0 baseT) as' end
1.333 +*}
1.334
1.335 +subsubsection {* Factor out gcd for cancellation *}
1.336 +ML {*
1.337 fun check_fraction t =
1.338 let val Const ("Fields.inverse_class.divide", _) $ numerator $ denominator = t
1.339 in SOME (numerator, denominator) end
1.340 @@ -334,13 +247,15 @@
1.341 (HOLogic.mk_binop "Groups.times_class.times"
1.342 (term_of_poly baseT expT vs a', ct),
1.343 HOLogic.mk_binop "Groups.times_class.times" (b't, ct))
1.344 - val asm = mk_asms baseT [b't, ct]
1.345 - in SOME (t', asm) end
1.346 + in SOME (t', mk_asms baseT [b't, ct]) end
1.347 end
1.348 | _ => NONE : (term * term list) option
1.349 end
1.350 end
1.351 +*}
1.352
1.353 +subsubsection {* Cancel a fraction *}
1.354 +ML {*
1.355 (* cancel a term by the gcd ("" denote terms with internal algebraic structure)
1.356 cancel_p_ :: theory \<Rightarrow> term \<Rightarrow> (term \<times> term list) option
1.357 cancel_p_ thy "a / b" = SOME ("a' / b'", ["b' \<noteq> 0"])
1.358 @@ -382,7 +297,10 @@
1.359 | _ => NONE : (term * term list) option
1.360 end
1.361 end
1.362 +*}
1.363
1.364 +subsubsection {* Factor out to a common denominator for addition *}
1.365 +ML {*
1.366 (* addition of fractions allows (at most) one non-fraction (a monomial) *)
1.367 fun check_frac_sum
1.368 (Const ("Groups.plus_class.plus", _) $
1.369 @@ -440,10 +358,14 @@
1.370 end
1.371 end
1.372
1.373 +*}
1.374 +
1.375 +subsubsection {* Addition of at least one fraction within a sum *}
1.376 +ML {*
1.377 (* add fractions
1.378 assumes: is a term with outmost "+" and at least one outmost "/" in respective summands
1.379 NOTE: the case "(_ + _) + _" need not be considered due to iterated addition.*)
1.380 -fun add_fraction_p_ (thy: theory) t =
1.381 +fun add_fraction_p_ (_: theory) t =
1.382 case check_frac_sum t of
1.383 NONE => NONE
1.384 | SOME ((n1, d1), (n2, d2)) =>
1.385 @@ -463,2512 +385,71 @@
1.386 in SOME (t', mk_asms baseT [denom]) end
1.387 | _ => NONE : (term * term list) option
1.388 end
1.389 +*}
1.390
1.391 -fun (x ins xs) = if x mem xs then xs else x :: xs;
1.392 -fun xs union [] = xs
1.393 - | [] union ys = ys
1.394 - | (x :: xs) union ys = xs union (x ins ys);
1.395 +section {* Embed cancellation and addition into rewriting *}
1.396 +ML {* val thy = @{theory} *}
1.397 +subsection {* Rulesets and predicate for embedding *}
1.398 +ML {*
1.399 +(* evaluates conditions in calculate_Rational *)
1.400 +val calc_rat_erls =
1.401 + prep_rls
1.402 + (Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.403 + erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.404 + rules =
1.405 + [Calc ("HOL.eq", eval_equal "#equal_"),
1.406 + Calc ("Atools.is'_const", eval_const "#is_const_"),
1.407 + Thm ("not_true", num_str @{thm not_true}),
1.408 + Thm ("not_false", num_str @{thm not_false})],
1.409 + scr = EmptyScr});
1.410
1.411 -(*. gcd of integers .*)
1.412 -(* die gcd Funktion von Isabelle funktioniert nicht richtig !!! *)
1.413 -fun gcd_int a b = if b=0 then a
1.414 - else gcd_int b (a mod b);
1.415 +(* simplifies expressions with numerals;
1.416 + does NOT rearrange the term by AC-rewriting; thus terms with variables
1.417 + need to have constants to be commuted together respectively *)
1.418 +val calculate_Rational =
1.419 + prep_rls (merge_rls "calculate_Rational"
1.420 + (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.421 + erls = calc_rat_erls, srls = Erls,
1.422 + calc = [], errpatts = [],
1.423 + rules =
1.424 + [Calc ("Fields.inverse_class.divide", eval_cancel "#divide_e"),
1.425
1.426 -(*. univariate polynomials (uv) .*)
1.427 -(*. univariate polynomials are represented as a list
1.428 - of the coefficent in reverse maximum degree order .*)
1.429 -(*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*)
1.430 -type uv_poly = int list;
1.431 -
1.432 -(*. adds two uv polynomials .*)
1.433 -fun uv_mod_add_poly ([]:uv_poly,p2:uv_poly) = p2:uv_poly
1.434 - | uv_mod_add_poly (p1,[]) = p1
1.435 - | uv_mod_add_poly (x::p1,y::p2) = (x+y)::(uv_mod_add_poly(p1,p2));
1.436 -
1.437 -(*. multiplies a uv polynomial with a skalar s .*)
1.438 -fun uv_mod_smul_poly ([]:uv_poly,s:int) = []:uv_poly
1.439 - | uv_mod_smul_poly (x::p,s) = (x*s)::(uv_mod_smul_poly(p,s));
1.440 -
1.441 -(*. calculates the remainder of a polynomial divided by a skalar s .*)
1.442 -fun uv_mod_rem_poly ([]:uv_poly,s) = []:uv_poly
1.443 - | uv_mod_rem_poly (x::p,s) = (x mod s)::(uv_mod_smul_poly(p,s));
1.444 -
1.445 -(*. calculates the degree of a uv polynomial .*)
1.446 -fun uv_mod_deg ([]:uv_poly) = 0
1.447 - | uv_mod_deg p = length(p)-1;
1.448 -
1.449 -(*. calculates the remainder of x/p and represents it as
1.450 - value between -p/2 and p/2 .*)
1.451 -fun uv_mod_mod2(x,p)=
1.452 - let
1.453 - val y=(x mod p);
1.454 - in
1.455 - if (y)>(p div 2) then (y)-p else
1.456 - (
1.457 - if (y)<(~p div 2) then p+(y) else (y)
1.458 - )
1.459 - end;
1.460 -
1.461 -(*.calculates the remainder for each element of a integer list divided by p.*)
1.462 -fun uv_mod_list_modp [] p = []
1.463 - | uv_mod_list_modp (x::xs) p = (uv_mod_mod2(x,p))::(uv_mod_list_modp xs p);
1.464 -
1.465 -(*. appends an integer at the end of a integer list .*)
1.466 -fun uv_mod_null (p1:int list,0) = p1
1.467 - | uv_mod_null (p1:int list,n1:int) = uv_mod_null(p1,n1-1) @ [0];
1.468 -
1.469 -(*. uv polynomial division, result is (quotient, remainder) .*)
1.470 -(*. only for uv_mod_divides .*)
1.471 -(* FIXME: Division von x^9+x^5+1 durch x-1000 funktioniert nicht,
1.472 - integer zu klein *)
1.473 -fun uv_mod_pdiv (p1:uv_poly) ([]:uv_poly) =
1.474 - error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")
1.475 - | uv_mod_pdiv p1 [x] =
1.476 - let
1.477 - val xs= Unsynchronized.ref [];
1.478 - in
1.479 - if x<>0 then
1.480 - (
1.481 - xs:=(uv_mod_rem_poly(p1,x));
1.482 - while length(!xs)>0 andalso hd(!xs)=0 do xs:=tl(!xs)
1.483 - )
1.484 - else error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero");
1.485 - ([]:uv_poly,!xs:uv_poly)
1.486 - end
1.487 - | uv_mod_pdiv p1 p2 =
1.488 - let
1.489 - val n= uv_mod_deg(p2);
1.490 - val m= Unsynchronized.ref (uv_mod_deg(p1));
1.491 - val p1'= Unsynchronized.ref (rev(p1));
1.492 - val p2'=(rev(p2));
1.493 - val lc2=hd(p2');
1.494 - val q= Unsynchronized.ref [];
1.495 - val c= Unsynchronized.ref 0;
1.496 - val output= Unsynchronized.ref ([],[]);
1.497 - in
1.498 - (
1.499 - if (!m)=0 orelse p2=[0]
1.500 - then error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero")
1.501 - else
1.502 - (
1.503 - if (!m)<n then
1.504 - (
1.505 - output:=([0],p1)
1.506 - )
1.507 - else
1.508 - (
1.509 - while (!m)>=n do
1.510 - (
1.511 - c:=hd(!p1') div hd(p2');
1.512 - if !c<>0 then
1.513 - (
1.514 - p1':=uv_mod_add_poly(!p1',uv_mod_null(uv_mod_smul_poly(p2',~(!c)),!m-n));
1.515 - while length(!p1')>0 andalso hd(!p1')=0 do p1':= tl(!p1');
1.516 - m:=uv_mod_deg(!p1')
1.517 - )
1.518 - else m:=0
1.519 - );
1.520 - output:=(rev(!q),rev(!p1'))
1.521 - )
1.522 - );
1.523 - !output
1.524 - )
1.525 - end;
1.526 -
1.527 -(*. divides p1 by p2 in Zp .*)
1.528 -fun uv_mod_pdivp (p1:uv_poly) (p2:uv_poly) p =
1.529 - let
1.530 - val n=uv_mod_deg(p2);
1.531 - val m= Unsynchronized.ref (uv_mod_deg(uv_mod_list_modp p1 p));
1.532 - val p1'= Unsynchronized.ref (rev(p1));
1.533 - val p2'=(rev(uv_mod_list_modp p2 p));
1.534 - val lc2=hd(p2');
1.535 - val q= Unsynchronized.ref [];
1.536 - val c= Unsynchronized.ref 0;
1.537 - val output= Unsynchronized.ref ([],[]);
1.538 - in
1.539 - (
1.540 - if (!m)=0 orelse p2=[0] then error ("RATIONALS_UV_MOD_PDIVP_EXCEPTION: Division by zero")
1.541 - else
1.542 - (
1.543 - if (!m)<n then
1.544 - (
1.545 - output:=([0],p1)
1.546 - )
1.547 - else
1.548 - (
1.549 - while (!m)>=n do
1.550 - (
1.551 - c:=uv_mod_mod2(hd(!p1')*(power lc2 1), p);
1.552 - q:=(!c)::(!q);
1.553 - p1':=uv_mod_list_modp(tl(uv_mod_add_poly(uv_mod_smul_poly(!p1',lc2),
1.554 - uv_mod_smul_poly(uv_mod_smul_poly(p2',hd(!p1')),~1)))) p;
1.555 - m:=(!m)-1
1.556 - );
1.557 -
1.558 - while !p1'<>[] andalso hd(!p1')=0 do
1.559 - (
1.560 - p1':=tl(!p1')
1.561 - );
1.562 -
1.563 - output:=(rev(uv_mod_list_modp (!q) (p)),rev(!p1'))
1.564 - )
1.565 - );
1.566 - !output:uv_poly * uv_poly
1.567 - )
1.568 - end;
1.569 -
1.570 -(*. calculates the remainder of p1/p2 .*)
1.571 -fun uv_mod_prest (p1:uv_poly) ([]:uv_poly) = error("UV_MOD_PREST_EXCEPTION: Division by zero")
1.572 - | uv_mod_prest [] p2 = []:uv_poly
1.573 - | uv_mod_prest p1 p2 = (#2(uv_mod_pdiv p1 p2));
1.574 -
1.575 -(*. calculates the remainder of p1/p2 in Zp .*)
1.576 -fun uv_mod_prestp (p1:uv_poly) ([]:uv_poly) p= error("UV_MOD_PRESTP_EXCEPTION: Division by zero")
1.577 - | uv_mod_prestp [] p2 p= []:uv_poly
1.578 - | uv_mod_prestp p1 p2 p = #2(uv_mod_pdivp p1 p2 p);
1.579 -
1.580 -(*. calculates the content of a uv polynomial .*)
1.581 -fun uv_mod_cont ([]:uv_poly) = 0
1.582 - | uv_mod_cont (x::p)= gcd_int x (uv_mod_cont(p));
1.583 -
1.584 -(*. divides each coefficient of a uv polynomial by y .*)
1.585 -fun uv_mod_div_list (p:uv_poly,0) = error("UV_MOD_DIV_LIST_EXCEPTION: Division by zero")
1.586 - | uv_mod_div_list ([],y) = []:uv_poly
1.587 - | uv_mod_div_list (x::p,y) = (x div y)::uv_mod_div_list(p,y);
1.588 -
1.589 -(*. calculates the primitiv part of a uv polynomial .*)
1.590 -fun uv_mod_pp ([]:uv_poly) = []:uv_poly
1.591 - | uv_mod_pp p =
1.592 - let
1.593 - val c= Unsynchronized.ref 0;
1.594 - in
1.595 - (
1.596 - c:=uv_mod_cont(p);
1.597 -
1.598 - if !c=0 then error ("RATIONALS_UV_MOD_PP_EXCEPTION: content is 0")
1.599 - else uv_mod_div_list(p,!c)
1.600 - )
1.601 - end;
1.602 -
1.603 -(*. gets the leading coefficient of a uv polynomial .*)
1.604 -fun uv_mod_lc ([]:uv_poly) = 0
1.605 - | uv_mod_lc p = hd(rev(p));
1.606 -
1.607 -(*. calculates the euklidean polynomial remainder sequence in Zp .*)
1.608 -fun uv_mod_prs_euklid_p(p1:uv_poly,p2:uv_poly,p)=
1.609 - let
1.610 - val f = Unsynchronized.ref [];
1.611 - val f'= Unsynchronized.ref p2;
1.612 - val fi= Unsynchronized.ref [];
1.613 - in
1.614 - (
1.615 - f:=p2::p1::[];
1.616 - while uv_mod_deg(!f')>0 do
1.617 - (
1.618 - f':=uv_mod_prestp (hd(tl(!f))) (hd(!f)) p;
1.619 - if (!f')<>[] then
1.620 - (
1.621 - fi:=(!f');
1.622 - f:=(!fi)::(!f)
1.623 - )
1.624 - else ()
1.625 - );
1.626 - (!f)
1.627 -
1.628 - )
1.629 - end;
1.630 -
1.631 -(*. calculates the gcd of p1 and p2 in Zp .*)
1.632 -fun uv_mod_gcd_modp ([]:uv_poly) (p2:uv_poly) p = p2:uv_poly
1.633 - | uv_mod_gcd_modp p1 [] p= p1
1.634 - | uv_mod_gcd_modp p1 p2 p=
1.635 - let
1.636 - val p1'= Unsynchronized.ref [];
1.637 - val p2'= Unsynchronized.ref [];
1.638 - val pc= Unsynchronized.ref [];
1.639 - val g= Unsynchronized.ref [];
1.640 - val d= Unsynchronized.ref 0;
1.641 - val prs= Unsynchronized.ref [];
1.642 - in
1.643 - (
1.644 - if uv_mod_deg(p1)>=uv_mod_deg(p2) then
1.645 - (
1.646 - p1':=uv_mod_list_modp (uv_mod_pp(p1)) p;
1.647 - p2':=uv_mod_list_modp (uv_mod_pp(p2)) p
1.648 - )
1.649 - else
1.650 - (
1.651 - p1':=uv_mod_list_modp (uv_mod_pp(p2)) p;
1.652 - p2':=uv_mod_list_modp (uv_mod_pp(p1)) p
1.653 - );
1.654 - d:=uv_mod_mod2((gcd_int (uv_mod_cont(p1))) (uv_mod_cont(p2)), p) ;
1.655 - if !d>(p div 2) then d:=(!d)-p else ();
1.656 -
1.657 - prs:=uv_mod_prs_euklid_p(!p1',!p2',p);
1.658 -
1.659 - if hd(!prs)=[] then pc:=hd(tl(!prs))
1.660 - else pc:=hd(!prs);
1.661 -
1.662 - g:=uv_mod_smul_poly(uv_mod_pp(!pc),!d);
1.663 - !g
1.664 - )
1.665 - end;
1.666 -
1.667 -(*. calculates the minimum of two real values x and y .*)
1.668 -fun uv_mod_r_min(x,y):Real.real = if x>y then y else x;
1.669 -
1.670 -(*. calculates the minimum of two integer values x and y .*)
1.671 -fun uv_mod_min(x,y) = if x>y then y else x;
1.672 -
1.673 -(*. adds the squared values of a integer list .*)
1.674 -fun uv_mod_add_qu [] = 0.0
1.675 - | uv_mod_add_qu (x::p) = Real.fromInt(x)*Real.fromInt(x) + uv_mod_add_qu p;
1.676 -
1.677 -(*. calculates the euklidean norm .*)
1.678 -fun uv_mod_norm ([]:uv_poly) = 0.0
1.679 - | uv_mod_norm p = Math.sqrt(uv_mod_add_qu(p));
1.680 -
1.681 -(*. multipies two values a and b .*)
1.682 -fun uv_mod_multi a b = a * b;
1.683 -
1.684 -(*. decides if x is a prim, the list contains all primes which are lower then x .*)
1.685 -fun uv_mod_prim(x,[])= false
1.686 - | uv_mod_prim(x,[y])=if ((x mod y) <> 0) then true
1.687 - else false
1.688 - | uv_mod_prim(x,y::ys) = if uv_mod_prim(x,[y])
1.689 - then
1.690 - if uv_mod_prim(x,ys) then true
1.691 - else false
1.692 - else false;
1.693 -
1.694 -(*. gets the first prime, which is greater than p and does not divide g .*)
1.695 -fun uv_mod_nextprime(g,p)=
1.696 - let
1.697 - val list= Unsynchronized.ref [2];
1.698 - val exit= Unsynchronized.ref 0;
1.699 - val i = Unsynchronized.ref 2
1.700 - in
1.701 - while (!i<p) do (* calculates the primes lower then p *)
1.702 - (
1.703 - if uv_mod_prim(!i,!list) then
1.704 - (
1.705 - if (p mod !i <> 0)
1.706 - then
1.707 - (
1.708 - list:= (!i)::(!list);
1.709 - i:= (!i)+1
1.710 - )
1.711 - else i:=(!i)+1
1.712 - )
1.713 - else i:= (!i)+1
1.714 - );
1.715 - i:=(p+1);
1.716 - while (!exit=0) do (* calculate next prime which does not divide g *)
1.717 - (
1.718 - if uv_mod_prim(!i,!list) then
1.719 - (
1.720 - if (g mod !i <> 0)
1.721 - then
1.722 - (
1.723 - list:= (!i)::(!list);
1.724 - exit:= (!i)
1.725 - )
1.726 - else i:=(!i)+1
1.727 - )
1.728 - else i:= (!i)+1
1.729 - );
1.730 - !exit
1.731 - end;
1.732 -
1.733 -(*. decides if p1 is a factor of p2 in Zp .*)
1.734 -fun uv_mod_dividesp ([]:uv_poly) (p2:uv_poly) p= error("UV_MOD_DIVIDESP: Division by zero")
1.735 - | uv_mod_dividesp p1 p2 p= if uv_mod_prestp p2 p1 p = [] then true else false;
1.736 -
1.737 -(*. decides if p1 is a factor of p2 .*)
1.738 -fun uv_mod_divides ([]:uv_poly) (p2:uv_poly) = error("UV_MOD_DIVIDES: Division by zero")
1.739 - | uv_mod_divides p1 p2 = if uv_mod_prest p2 p1 = [] then true else false;
1.740 -
1.741 -(*. chinese remainder algorithm .*)
1.742 -fun uv_mod_cra2(r1,r2,m1,m2)=
1.743 - let
1.744 - val c= Unsynchronized.ref 0;
1.745 - val r1'= Unsynchronized.ref 0;
1.746 - val d= Unsynchronized.ref 0;
1.747 - val a= Unsynchronized.ref 0;
1.748 - in
1.749 - (
1.750 - while (uv_mod_mod2((!c)*m1,m2))<>1 do
1.751 - (
1.752 - c:=(!c)+1
1.753 - );
1.754 - r1':= uv_mod_mod2(r1,m1);
1.755 - d:=uv_mod_mod2(((r2-(!r1'))*(!c)),m2);
1.756 - !r1'+(!d)*m1
1.757 - )
1.758 - end;
1.759 -
1.760 -(*. applies the chinese remainder algorithmen to the coefficients of x1 and x2 .*)
1.761 -fun uv_mod_cra_2 ([],[],m1,m2) = []
1.762 - | uv_mod_cra_2 ([],x2,m1,m2) = error("UV_MOD_CRA_2_EXCEPTION: invalid call x1")
1.763 - | uv_mod_cra_2 (x1,[],m1,m2) = error("UV_MOD_CRA_2_EXCEPTION: invalid call x2")
1.764 - | uv_mod_cra_2 (x1::x1s,x2::x2s,m1,m2) = (uv_mod_cra2(x1,x2,m1,m2))::(uv_mod_cra_2(x1s,x2s,m1,m2));
1.765 -
1.766 -(*. calculates the gcd of two uv polynomials p1' and p2' with the modular algorithm .*)
1.767 -fun uv_mod_gcd (p1':uv_poly) (p2':uv_poly) =
1.768 - let
1.769 - val p1= Unsynchronized.ref (uv_mod_pp(p1'));
1.770 - val p2= Unsynchronized.ref (uv_mod_pp(p2'));
1.771 - val c=gcd_int (uv_mod_cont(p1')) (uv_mod_cont(p2'));
1.772 - val temp= Unsynchronized.ref [];
1.773 - val cp= Unsynchronized.ref [];
1.774 - val qp= Unsynchronized.ref [];
1.775 - val q= Unsynchronized.ref [];
1.776 - val pn= Unsynchronized.ref 0;
1.777 - val d= Unsynchronized.ref 0;
1.778 - val g1= Unsynchronized.ref 0;
1.779 - val p= Unsynchronized.ref 0;
1.780 - val m= Unsynchronized.ref 0;
1.781 - val exit= Unsynchronized.ref 0;
1.782 - val i= Unsynchronized.ref 1;
1.783 - in
1.784 - if length(!p1)>length(!p2) then ()
1.785 - else
1.786 - (
1.787 - temp:= !p1;
1.788 - p1:= !p2;
1.789 - p2:= !temp
1.790 - );
1.791 -
1.792 -
1.793 - d:=gcd_int (uv_mod_lc(!p1)) (uv_mod_lc(!p2));
1.794 - g1:=uv_mod_lc(!p1)*uv_mod_lc(!p2);
1.795 - p:=4;
1.796 -
1.797 - m:=Real.ceil(2.0 * Real.fromInt(!d) *
1.798 - Real.fromInt(power 2 (uv_mod_min(uv_mod_deg(!p2),uv_mod_deg(!p1)))) *
1.799 - Real.fromInt(!d) *
1.800 - uv_mod_r_min(uv_mod_norm(!p1) / Real.fromInt(abs(uv_mod_lc(!p1))),
1.801 - uv_mod_norm(!p2) / Real.fromInt(abs(uv_mod_lc(!p2)))));
1.802 -
1.803 - while (!exit=0) do
1.804 - (
1.805 - p:=uv_mod_nextprime(!d,!p);
1.806 - cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)) ;
1.807 - if abs(uv_mod_lc(!cp))<>1 then (* leading coefficient = 1 ? *)
1.808 - (
1.809 - i:=1;
1.810 - while (!i)<(!p) andalso (abs(uv_mod_mod2((uv_mod_lc(!cp)*(!i)),(!p)))<>1) do
1.811 - (
1.812 - i:=(!i)+1
1.813 - );
1.814 - cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)
1.815 - )
1.816 - else ();
1.817 -
1.818 - qp:= ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp));
1.819 -
1.820 - if uv_mod_deg(!qp)=0 then (q:=[1]; exit:=1) else ();
1.821 -
1.822 - pn:=(!p);
1.823 - q:=(!qp);
1.824 -
1.825 - while !pn<= !m andalso !m>(!p) andalso !exit=0 do
1.826 - (
1.827 - p:=uv_mod_nextprime(!d,!p);
1.828 - cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p));
1.829 - if uv_mod_lc(!cp)<>1 then (* leading coefficient = 1 ? *)
1.830 - (
1.831 - i:=1;
1.832 - while (!i)<(!p) andalso ((uv_mod_mod2((uv_mod_lc(!q)*(!i)),(!p)))<>1) do
1.833 - (
1.834 - i:=(!i)+1
1.835 - );
1.836 - cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)
1.837 - )
1.838 - else ();
1.839 -
1.840 - qp:=uv_mod_list_modp ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp) ) (!p);
1.841 - if uv_mod_deg(!qp)>uv_mod_deg(!q) then
1.842 - (
1.843 - (*print("degree to high!!!\n")*)
1.844 - )
1.845 - else
1.846 - (
1.847 - if uv_mod_deg(!qp)=uv_mod_deg(!q) then
1.848 - (
1.849 - q:=uv_mod_cra_2(!q,!qp,!pn,!p);
1.850 - pn:=(!pn) * !p;
1.851 - q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn)); (* found already gcd ? *)
1.852 - if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then (exit:=1) else ()
1.853 - )
1.854 - else
1.855 - (
1.856 - if uv_mod_deg(!qp)<uv_mod_deg(!q) then
1.857 - (
1.858 - pn:= !p;
1.859 - q:= !qp
1.860 - )
1.861 - else ()
1.862 - )
1.863 - )
1.864 - );
1.865 - q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn));
1.866 - if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then exit:=1 else ()
1.867 - );
1.868 - uv_mod_smul_poly(!q,c):uv_poly
1.869 - end;
1.870 -
1.871 -(*. multivariate polynomials .*)
1.872 -(*. multivariate polynomials are represented as a list of the pairs,
1.873 - first is the coefficent and the second is a list of the exponents .*)
1.874 -(*. 5 * x^5 * y^3 + 4 * x^3 * z^2 + 2 * x^2 * y * z^3 - z - 19
1.875 - => [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*)
1.876 -
1.877 -(*. global variables .*)
1.878 -(*. order indicators .*)
1.879 -val LEX_=0; (* lexicographical term order *)
1.880 -val GGO_=1; (* greatest degree order *)
1.881 -
1.882 -(*. datatypes for internal representation.*)
1.883 -type mv_monom = (int * (*.coefficient or the monom.*)
1.884 - int list); (*.list of exponents) .*)
1.885 -fun mv_monom2str (i, is) = "("^ int2str i^"," ^ ints2str' is ^ ")";
1.886 -
1.887 -type mv_poly = mv_monom list;
1.888 -fun mv_poly2str p = (strs2str' o (map mv_monom2str)) p;
1.889 -
1.890 -(*. help function for monom_greater and geq .*)
1.891 -fun mv_mg_hlp([]) = EQUAL
1.892 - | mv_mg_hlp(x::list)=if x<0 then LESS
1.893 - else if x>0 then GREATER
1.894 - else mv_mg_hlp(list);
1.895 -
1.896 -(*. adds a list of values .*)
1.897 -fun mv_addlist([]) = 0
1.898 - | mv_addlist(p1) = hd(p1)+mv_addlist(tl(p1));
1.899 -
1.900 -(*. tests if the monomial M1 is greater as the monomial M2 and returns a boolean value .*)
1.901 -(*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
1.902 -fun mv_monom_greater((M1x,M1l):mv_monom,(M2x,M2l):mv_monom,order)=
1.903 - if order=LEX_ then
1.904 - (
1.905 - if length(M1l)<>length(M2l) then error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
1.906 - else if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
1.907 - )
1.908 - else
1.909 - if order=GGO_ then
1.910 - (
1.911 - if length(M1l)<>length(M2l) then error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
1.912 - else
1.913 - if mv_addlist(M1l)=mv_addlist(M2l) then if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
1.914 - else if mv_addlist(M1l)>mv_addlist(M2l) then true else false
1.915 - )
1.916 - else error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");
1.917 -
1.918 -(*. tests if the monomial X is greater as the monomial Y and returns a order value (GREATER,EQUAL,LESS) .*)
1.919 -(*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
1.920 -fun mv_geq order ((x1,x):mv_monom,(x2,y):mv_monom) =
1.921 -let
1.922 - val temp= Unsynchronized.ref EQUAL;
1.923 -in
1.924 - if order=LEX_ then
1.925 - (
1.926 - if length(x)<>length(y) then
1.927 - error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
1.928 - else
1.929 - (
1.930 - temp:=mv_mg_hlp((map op- (x~~y)));
1.931 - if !temp=EQUAL then
1.932 - ( if x1=x2 then EQUAL
1.933 - else if x1>x2 then GREATER
1.934 - else LESS
1.935 - )
1.936 - else (!temp)
1.937 - )
1.938 - )
1.939 - else
1.940 - if order=GGO_ then
1.941 - (
1.942 - if length(x)<>length(y) then
1.943 - error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
1.944 - else
1.945 - if mv_addlist(x)=mv_addlist(y) then
1.946 - (mv_mg_hlp((map op- (x~~y))))
1.947 - else if mv_addlist(x)>mv_addlist(y) then GREATER else LESS
1.948 - )
1.949 - else error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")
1.950 -end;
1.951 -
1.952 -(*. cuts the first variable from a polynomial .*)
1.953 -fun mv_cut([]:mv_poly)=[]:mv_poly
1.954 - | mv_cut((x,[])::list) = error ("RATIONALS_MV_CUT_EXCEPTION: Invalid list ")
1.955 - | mv_cut((x,y::ys)::list)=(x,ys)::mv_cut(list);
1.956 -
1.957 -(*. leading power product .*)
1.958 -fun mv_lpp([]:mv_poly,order) = []
1.959 - | mv_lpp([(x,y)],order) = y
1.960 - | mv_lpp(p1,order) = #2(hd(rev(sort (mv_geq order) p1)));
1.961 -
1.962 -(*. leading monomial .*)
1.963 -fun mv_lm([]:mv_poly,order) = (0,[]):mv_monom
1.964 - | mv_lm([x],order) = x
1.965 - | mv_lm(p1,order) = hd(rev(sort (mv_geq order) p1));
1.966 -
1.967 -(*. leading coefficient in term order .*)
1.968 -fun mv_lc2([]:mv_poly,order) = 0
1.969 - | mv_lc2([(x,y)],order) = x
1.970 - | mv_lc2(p1,order) = #1(hd(rev(sort (mv_geq order) p1)));
1.971 -
1.972 -
1.973 -(*. reverse the coefficients in mv polynomial .*)
1.974 -fun mv_rev_to([]:mv_poly) = []:mv_poly
1.975 - | mv_rev_to((c,e)::xs) = (c,rev(e))::mv_rev_to(xs);
1.976 -
1.977 -(*. leading coefficient in reverse term order .*)
1.978 -fun mv_lc([]:mv_poly,order) = []:mv_poly
1.979 - | mv_lc([(x,y)],order) = mv_rev_to(mv_cut(mv_rev_to([(x,y)])))
1.980 - | mv_lc(p1,order) =
1.981 - let
1.982 - val p1o= Unsynchronized.ref (rev(sort (mv_geq order) (mv_rev_to(p1))));
1.983 - val lp=hd(#2(hd(!p1o)));
1.984 - val lc= Unsynchronized.ref [];
1.985 - in
1.986 - (
1.987 - while (length(!p1o)>0 andalso hd(#2(hd(!p1o)))=lp) do
1.988 - (
1.989 - lc:=hd(mv_cut([hd(!p1o)]))::(!lc);
1.990 - p1o:=tl(!p1o)
1.991 - );
1.992 - if !lc=[] then error ("RATIONALS_MV_LC_EXCEPTION: lc is empty") else ();
1.993 - mv_rev_to(!lc)
1.994 - )
1.995 - end;
1.996 -
1.997 -(*. compares two powerproducts .*)
1.998 -fun mv_monom_equal((_,xlist):mv_monom,(_,ylist):mv_monom) = (foldr and_) (((map op=) (xlist~~ylist)),true);
1.999 -
1.1000 -(*. help function for mv_add .*)
1.1001 -fun mv_madd([]:mv_poly,[]:mv_poly,order) = []:mv_poly
1.1002 - | mv_madd([(0,_)],p2,order) = p2
1.1003 - | mv_madd(p1,[(0,_)],order) = p1
1.1004 - | mv_madd([],p2,order) = p2
1.1005 - | mv_madd(p1,[],order) = p1
1.1006 - | mv_madd(p1,p2,order) =
1.1007 - (
1.1008 - if mv_monom_greater(hd(p1),hd(p2),order)
1.1009 - then hd(p1)::mv_madd(tl(p1),p2,order)
1.1010 - else if mv_monom_equal(hd(p1),hd(p2))
1.1011 - then if mv_lc2(p1,order)+mv_lc2(p2,order)<>0
1.1012 - then (mv_lc2(p1,order)+mv_lc2(p2,order),mv_lpp(p1,order))::mv_madd(tl(p1),tl(p2),order)
1.1013 - else mv_madd(tl(p1),tl(p2),order)
1.1014 - else hd(p2)::mv_madd(p1,tl(p2),order)
1.1015 - )
1.1016 -
1.1017 -(*. adds two multivariate polynomials .*)
1.1018 -fun mv_add([]:mv_poly,p2:mv_poly,order) = p2
1.1019 - | mv_add(p1,[],order) = p1
1.1020 - | mv_add(p1,p2,order) = mv_madd(rev(sort (mv_geq order) p1),rev(sort (mv_geq order) p2), order);
1.1021 -
1.1022 -(*. monom multiplication .*)
1.1023 -fun mv_mmul((x1,y1):mv_monom,(x2,y2):mv_monom)=(x1*x2,(map op+) (y1~~y2)):mv_monom;
1.1024 -
1.1025 -(*. deletes all monomials with coefficient 0 .*)
1.1026 -fun mv_shorten([]:mv_poly,order) = []:mv_poly
1.1027 - | mv_shorten(x::xs,order)=mv_madd([x],mv_shorten(xs,order),order);
1.1028 -
1.1029 -(*. zeros a list .*)
1.1030 -fun mv_null2([])=[]
1.1031 - | mv_null2(x::l)=0::mv_null2(l);
1.1032 -
1.1033 -(*. multiplies two multivariate polynomials .*)
1.1034 -fun mv_mul([]:mv_poly,[]:mv_poly,_) = []:mv_poly
1.1035 - | mv_mul([],y::p2,_) = [(0,mv_null2(#2(y)))]
1.1036 - | mv_mul(x::p1,[],_) = [(0,mv_null2(#2(x)))]
1.1037 - | mv_mul(x::p1,y::p2,order) = mv_shorten(rev(sort (mv_geq order) (mv_mmul(x,y) :: (mv_mul(p1,y::p2,order) @
1.1038 - mv_mul([x],p2,order)))),order);
1.1039 -
1.1040 -(*. gets the maximum value of a list .*)
1.1041 -fun mv_getmax([])=0
1.1042 - | mv_getmax(x::p1)= let
1.1043 - val m=mv_getmax(p1);
1.1044 - in
1.1045 - if m>x then m
1.1046 - else x
1.1047 - end;
1.1048 -(*. calculates the maximum degree of an multivariate polynomial .*)
1.1049 -fun mv_grad([]:mv_poly) = 0
1.1050 - | mv_grad(p1:mv_poly)= mv_getmax((map mv_addlist) ((map #2) p1));
1.1051 -
1.1052 -(*. converts the sign of a value .*)
1.1053 -fun mv_minus(x)=(~1) * x;
1.1054 -
1.1055 -(*. converts the sign of all coefficients of a polynomial .*)
1.1056 -fun mv_minus2([]:mv_poly)=[]:mv_poly
1.1057 - | mv_minus2(p1)=(mv_minus(#1(hd(p1))),#2(hd(p1)))::(mv_minus2(tl(p1)));
1.1058 -
1.1059 -(*. searches for a negativ value in a list .*)
1.1060 -fun mv_is_negativ([])=false
1.1061 - | mv_is_negativ(x::xs)=if x<0 then true else mv_is_negativ(xs);
1.1062 -
1.1063 -(*. division of monomials .*)
1.1064 -fun mv_mdiv((0,[]):mv_monom,_:mv_monom)=(0,[]):mv_monom
1.1065 - | mv_mdiv(_,(0,[]))= error ("RATIONALS_MV_MDIV_EXCEPTION Division by 0 ")
1.1066 - | mv_mdiv(p1:mv_monom,p2:mv_monom)=
1.1067 - let
1.1068 - val c= Unsynchronized.ref (#1(p2));
1.1069 - val pp= Unsynchronized.ref [];
1.1070 - in
1.1071 - (
1.1072 - if !c=0 then error("MV_MDIV_EXCEPTION Dividing by zero")
1.1073 - else c:=(#1(p1) div #1(p2));
1.1074 - if #1(p2)<>0 then
1.1075 - (
1.1076 - pp:=(#2(mv_mmul((1,#2(p1)),(1,(map mv_minus) (#2(p2))))));
1.1077 - if mv_is_negativ(!pp) then (0,!pp)
1.1078 - else (!c,!pp)
1.1079 - )
1.1080 - else error("MV_MDIV_EXCEPTION Dividing by empty Polynom")
1.1081 - )
1.1082 - end;
1.1083 -
1.1084 -(*. prints a polynom for (internal use only) .*)
1.1085 -fun mv_print_poly([]:mv_poly)=tracing("[]\n")
1.1086 - | mv_print_poly((x,y)::[])= tracing("("^Int.toString(x)^","^ints2str(y)^")\n")
1.1087 - | mv_print_poly((x,y)::p1) = (tracing("("^Int.toString(x)^","^ints2str(y)^"),");mv_print_poly(p1));
1.1088 -
1.1089 -
1.1090 -(*. division of two multivariate polynomials .*)
1.1091 -fun mv_division([]:mv_poly,g:mv_poly,order)=([]:mv_poly,[]:mv_poly)
1.1092 - | mv_division(f,[],order)= error ("RATIONALS_MV_DIVISION_EXCEPTION Division by zero")
1.1093 - | mv_division(f,g,order)=
1.1094 - let
1.1095 - val r= Unsynchronized.ref [];
1.1096 - val q= Unsynchronized.ref [];
1.1097 - val g'= Unsynchronized.ref ([] : mv_monom list);
1.1098 - val k= Unsynchronized.ref 0;
1.1099 - val m= Unsynchronized.ref (0,[0]);
1.1100 - val exit= Unsynchronized.ref 0;
1.1101 - in
1.1102 - r := rev(sort (mv_geq order) (mv_shorten(f,order)));
1.1103 - g':= rev(sort (mv_geq order) (mv_shorten(g,order)));
1.1104 - if #1(hd(!g'))=0 then error("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero") else ();
1.1105 - if (mv_monom_greater (hd(!g'),hd(!r),order)) then ([(0,mv_null2(#2(hd(f))))],(!r))
1.1106 - else
1.1107 - (
1.1108 - exit:=0;
1.1109 - while (if (!exit)=0 then not(mv_monom_greater (hd(!g'),hd(!r),order)) else false) do
1.1110 - (
1.1111 - if (#1(mv_lm(!g',order)))<>0 then m:=mv_mdiv(mv_lm(!r,order),mv_lm(!g',order))
1.1112 - else error ("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero");
1.1113 - if #1(!m)<>0 then
1.1114 - (
1.1115 - q:=(!m)::(!q);
1.1116 - r:=mv_add((!r),mv_minus2(mv_mul(!g',[!m],order)),order)
1.1117 - )
1.1118 - else exit:=1;
1.1119 - if (if length(!r)<>0 then length(!g')<>0 else false) then ()
1.1120 - else (exit:=1)
1.1121 - );
1.1122 - (rev(!q),!r)
1.1123 - )
1.1124 - end;
1.1125 -
1.1126 -(*. multiplies a polynomial with an integer .*)
1.1127 -fun mv_skalar_mul([]:mv_poly,c) = []:mv_poly
1.1128 - | mv_skalar_mul((x,y)::p1,c) = ((x * c),y)::mv_skalar_mul(p1,c);
1.1129 -
1.1130 -(*. inserts the a first variable into an polynomial with exponent v .*)
1.1131 -fun mv_correct([]:mv_poly,v:int)=[]:mv_poly
1.1132 - | mv_correct((x,y)::list,v:int)=(x,v::y)::mv_correct(list,v);
1.1133 -
1.1134 -(*. multivariate case .*)
1.1135 -
1.1136 -(*. decides if x is a factor of y .*)
1.1137 -fun mv_divides([]:mv_poly,[]:mv_poly)= error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
1.1138 - | mv_divides(x,[]) = error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
1.1139 - | mv_divides(x:mv_poly,y:mv_poly) = #2(mv_division(y,x,LEX_))=[];
1.1140 -
1.1141 -(*. gets the maximum of a and b .*)
1.1142 -fun mv_max(a,b) = if a>b then a else b;
1.1143 -
1.1144 -(*. gets the maximum exponent of a mv polynomial in the lexicographic term order .*)
1.1145 -fun mv_deg([]:mv_poly) = 0
1.1146 - | mv_deg(p1)=
1.1147 - let
1.1148 - val p1'=mv_shorten(p1,LEX_);
1.1149 - in
1.1150 - if length(p1')=0 then 0
1.1151 - else mv_max(hd(#2(hd(p1'))),mv_deg(tl(p1')))
1.1152 - end;
1.1153 -
1.1154 -(*. gets the maximum exponent of a mv polynomial in the reverse lexicographic term order .*)
1.1155 -fun mv_deg2([]:mv_poly) = 0
1.1156 - | mv_deg2(p1)=
1.1157 - let
1.1158 - val p1'=mv_shorten(p1,LEX_);
1.1159 - in
1.1160 - if length(p1')=0 then 0
1.1161 - else mv_max(hd(rev(#2(hd(p1')))),mv_deg2(tl(p1')))
1.1162 - end;
1.1163 -
1.1164 -(*. evaluates the mv polynomial at the value v of the main variable .*)
1.1165 -fun mv_subs([]:mv_poly,v) = []:mv_poly
1.1166 - | mv_subs((c,e)::p1:mv_poly,v) = mv_skalar_mul(mv_cut([(c,e)]),power v (hd(e))) @ mv_subs(p1,v);
1.1167 -
1.1168 -(*. calculates the content of a uv-polynomial in mv-representation .*)
1.1169 -fun uv_content2([]:mv_poly) = 0
1.1170 - | uv_content2((c,e)::p1) = (gcd_int c (uv_content2(p1)));
1.1171 -
1.1172 -(*. converts a uv-polynomial from mv-representation to uv-representation .*)
1.1173 -fun uv_to_list ([]:mv_poly)=[]:uv_poly
1.1174 - | uv_to_list ((c1,e1)::others) =
1.1175 - let
1.1176 - val count= Unsynchronized.ref 0;
1.1177 - val max=mv_grad((c1,e1)::others);
1.1178 - val help= Unsynchronized.ref ((c1,e1)::others);
1.1179 - val list= Unsynchronized.ref [];
1.1180 - in
1.1181 - if length(e1)>1 then error ("RATIONALS_TO_LIST_EXCEPTION: not univariate")
1.1182 - else if length(e1)=0 then [c1]
1.1183 - else
1.1184 - (
1.1185 - count:=0;
1.1186 - while (!count)<=max do
1.1187 - (
1.1188 - if length(!help)>0 andalso hd(#2(hd(!help)))=max-(!count) then
1.1189 - (
1.1190 - list:=(#1(hd(!help)))::(!list);
1.1191 - help:=tl(!help)
1.1192 - )
1.1193 - else
1.1194 - (
1.1195 - list:= 0::(!list)
1.1196 - );
1.1197 - count := (!count) + 1
1.1198 - );
1.1199 - (!list)
1.1200 - )
1.1201 - end;
1.1202 -
1.1203 -(*. converts a uv-polynomial from uv-representation to mv-representation .*)
1.1204 -fun uv_to_poly ([]:uv_poly) = []:mv_poly
1.1205 - | uv_to_poly p1 =
1.1206 - let
1.1207 - val count= Unsynchronized.ref 0;
1.1208 - val help= Unsynchronized.ref p1;
1.1209 - val list= Unsynchronized.ref [];
1.1210 - in
1.1211 - while length(!help)>0 do
1.1212 - (
1.1213 - if hd(!help)=0 then ()
1.1214 - else list:=(hd(!help),[!count])::(!list);
1.1215 - count:=(!count)+1;
1.1216 - help:=tl(!help)
1.1217 - );
1.1218 - (!list)
1.1219 - end;
1.1220 -
1.1221 -(*. univariate gcd calculation from polynomials in multivariate representation .*)
1.1222 -fun uv_gcd ([]:mv_poly) p2 = p2
1.1223 - | uv_gcd p1 ([]:mv_poly) = p1
1.1224 - | uv_gcd p1 [(c,[e])] =
1.1225 - let
1.1226 - val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
1.1227 - val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
1.1228 - in
1.1229 - [(gcd_int (uv_content2(p1)) c,[min])]
1.1230 - end
1.1231 - | uv_gcd [(c,[e])] p2 =
1.1232 - let
1.1233 - val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p2,LEX_))));
1.1234 - val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
1.1235 - in
1.1236 - [(gcd_int (uv_content2(p2)) c,[min])]
1.1237 - end
1.1238 - | uv_gcd p11 p22 = uv_to_poly(uv_mod_gcd (uv_to_list(mv_shorten(p11,LEX_))) (uv_to_list(mv_shorten(p22,LEX_))));
1.1239 -
1.1240 -(*. help function for the newton interpolation .*)
1.1241 -fun mv_newton_help ([]:mv_poly list,k:int) = []:mv_poly list
1.1242 - | mv_newton_help (pl:mv_poly list,k) =
1.1243 - let
1.1244 - val x= Unsynchronized.ref (rev(pl));
1.1245 - val t= Unsynchronized.ref [];
1.1246 - val y= Unsynchronized.ref [];
1.1247 - val n= Unsynchronized.ref 1;
1.1248 - val n1= Unsynchronized.ref [];
1.1249 - in
1.1250 - (
1.1251 - while length(!x)>1 do
1.1252 - (
1.1253 - if length(hd(!x))>0 then n1:=mv_null2(#2(hd(hd(!x))))
1.1254 - else if length(hd(tl(!x)))>0 then n1:=mv_null2(#2(hd(hd(tl(!x)))))
1.1255 - else n1:=[];
1.1256 - t:= #1(mv_division(mv_add(hd(!x),mv_skalar_mul(hd(tl(!x)),~1),LEX_),[(k,!n1)],LEX_));
1.1257 - y:=(!t)::(!y);
1.1258 - x:=tl(!x)
1.1259 - );
1.1260 - (!y)
1.1261 - )
1.1262 - end;
1.1263 -
1.1264 -(*. help function for the newton interpolation .*)
1.1265 -fun mv_newton_add ([]:mv_poly list) t = []:mv_poly
1.1266 - | mv_newton_add [x:mv_poly] t = x
1.1267 - | mv_newton_add (pl:mv_poly list) t =
1.1268 - let
1.1269 - val expos= Unsynchronized.ref [];
1.1270 - val pll= Unsynchronized.ref pl;
1.1271 - in
1.1272 - (
1.1273 -
1.1274 - while length(!pll)>0 andalso hd(!pll)=[] do
1.1275 - (
1.1276 - pll:=tl(!pll)
1.1277 - );
1.1278 - if length(!pll)>0 then expos:= #2(hd(hd(!pll))) else expos:=[];
1.1279 - mv_add(hd(pl),
1.1280 - mv_mul(
1.1281 - mv_add(mv_correct(mv_cut([(1,mv_null2(!expos))]),1),[(~t,mv_null2(!expos))],LEX_),
1.1282 - mv_newton_add (tl(pl)) (t+1),
1.1283 - LEX_
1.1284 - ),
1.1285 - LEX_)
1.1286 - )
1.1287 - end;
1.1288 -
1.1289 -(*. calculates the newton interpolation with polynomial coefficients .*)
1.1290 -(*. step-depth is 1 and if the result is not an integerpolynomial .*)
1.1291 -(*. this function returns [] .*)
1.1292 -fun mv_newton ([]:(mv_poly) list) = []:mv_poly
1.1293 - | mv_newton ([mp]:(mv_poly) list) = mp:mv_poly
1.1294 - | mv_newton pl =
1.1295 - let
1.1296 - val c= Unsynchronized.ref pl;
1.1297 - val c1= Unsynchronized.ref [];
1.1298 - val n=length(pl);
1.1299 - val k= Unsynchronized.ref 1;
1.1300 - val i= Unsynchronized.ref n;
1.1301 - val ppl= Unsynchronized.ref [];
1.1302 - in
1.1303 - c1:=hd(pl)::[];
1.1304 - c:=mv_newton_help(!c,!k);
1.1305 - c1:=(hd(!c))::(!c1);
1.1306 - while(length(!c)>1 andalso !k<n) do
1.1307 - (
1.1308 - k:=(!k)+1;
1.1309 - while length(!c)>0 andalso hd(!c)=[] do c:=tl(!c);
1.1310 - if !c=[] then () else c:=mv_newton_help(!c,!k);
1.1311 - ppl:= !c;
1.1312 - if !c=[] then () else c1:=(hd(!c))::(!c1)
1.1313 - );
1.1314 - while hd(!c1)=[] do c1:=tl(!c1);
1.1315 - c1:=rev(!c1);
1.1316 - ppl:= !c1;
1.1317 - mv_newton_add (!c1) 1
1.1318 - end;
1.1319 -
1.1320 -(*. sets the exponents of the first variable to zero .*)
1.1321 -fun mv_null3([]:mv_poly) = []:mv_poly
1.1322 - | mv_null3((x,y)::xs) = (x,0::tl(y))::mv_null3(xs);
1.1323 -
1.1324 -(*. calculates the minimum exponents of a multivariate polynomial .*)
1.1325 -fun mv_min_pp([]:mv_poly)=[]
1.1326 - | mv_min_pp((c,e)::xs)=
1.1327 - let
1.1328 - val y= Unsynchronized.ref xs;
1.1329 - val x= Unsynchronized.ref [];
1.1330 - in
1.1331 - (
1.1332 - x:=e;
1.1333 - while length(!y)>0 do
1.1334 - (
1.1335 - x:=(map uv_mod_min) ((!x) ~~ (#2(hd(!y))));
1.1336 - y:=tl(!y)
1.1337 - );
1.1338 - !x
1.1339 - )
1.1340 - end;
1.1341 -
1.1342 -(*. checks if all elements of the list have value zero .*)
1.1343 -fun list_is_null [] = true
1.1344 - | list_is_null (x::xs) = (x=0 andalso list_is_null(xs));
1.1345 -
1.1346 -(* check if main variable is zero*)
1.1347 -fun main_zero (ms : mv_poly) = (list_is_null o (map (hd o #2))) ms;
1.1348 -
1.1349 -(*. calculates the content of an polynomial .*)
1.1350 -fun mv_content([]:mv_poly) = []:mv_poly
1.1351 - | mv_content(p1) =
1.1352 - let
1.1353 - val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
1.1354 - val test= Unsynchronized.ref (hd(#2(hd(!list))));
1.1355 - val result= Unsynchronized.ref [];
1.1356 - val min=(hd(#2(hd(rev(!list)))));
1.1357 - in
1.1358 - (
1.1359 - if length(!list)>1 then
1.1360 - (
1.1361 - while (if length(!list)>0 then (hd(#2(hd(!list)))=(!test)) else false) do
1.1362 - (
1.1363 - result:=(#1(hd(!list)),tl(#2(hd(!list))))::(!result);
1.1364 -
1.1365 - if length(!list)<1 then list:=[]
1.1366 - else list:=tl(!list)
1.1367 -
1.1368 - );
1.1369 - if length(!list)>0 then
1.1370 - (
1.1371 - list:=mv_gcd (!result) (mv_cut(mv_content(!list)))
1.1372 - )
1.1373 - else list:=(!result);
1.1374 - list:=mv_correct(!list,0);
1.1375 - (!list)
1.1376 - )
1.1377 - else
1.1378 - (
1.1379 - mv_null3(!list)
1.1380 - )
1.1381 - )
1.1382 - end
1.1383 -
1.1384 -(*. calculates the primitiv part of a polynomial .*)
1.1385 -and mv_pp([]:mv_poly) = []:mv_poly
1.1386 - | mv_pp(p1) = let
1.1387 - val cont= Unsynchronized.ref [];
1.1388 - val pp= Unsynchronized.ref [];
1.1389 - in
1.1390 - cont:=mv_content(p1);
1.1391 - pp:=(#1(mv_division(p1,!cont,LEX_)));
1.1392 - if !pp=[]
1.1393 - then error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")
1.1394 - else (!pp)
1.1395 - end
1.1396 -
1.1397 -(*. calculates the gcd of two multivariate polynomials with a modular approach .*)
1.1398 -and mv_gcd ([]:mv_poly) ([]:mv_poly) :mv_poly= []:mv_poly
1.1399 - | mv_gcd ([]:mv_poly) (p2) :mv_poly= p2:mv_poly
1.1400 - | mv_gcd (p1:mv_poly) ([]) :mv_poly= p1:mv_poly
1.1401 - | mv_gcd ([(x,xs)]:mv_poly) ([(y,ys)]):mv_poly =
1.1402 - let
1.1403 - val xpoly:mv_poly = [(x,xs)];
1.1404 - val ypoly:mv_poly = [(y,ys)];
1.1405 - in
1.1406 - (
1.1407 - if xs=ys then [((gcd_int x y),xs)]
1.1408 - else [((gcd_int x y),(map uv_mod_min)(xs~~ys))]:mv_poly
1.1409 - )
1.1410 - end
1.1411 - | mv_gcd (p1:mv_poly) ([(y,ys)]) :mv_poly=
1.1412 - (
1.1413 - [(gcd_int (uv_content2(p1)) (y),(map uv_mod_min)(mv_min_pp(p1)~~ys))]:mv_poly
1.1414 - )
1.1415 - | mv_gcd ([(y,ys)]:mv_poly) (p2):mv_poly =
1.1416 - (
1.1417 - [(gcd_int (uv_content2(p2)) (y),(map uv_mod_min)(mv_min_pp(p2)~~ys))]:mv_poly
1.1418 - )
1.1419 - | mv_gcd (p1':mv_poly) (p2':mv_poly):mv_poly=
1.1420 - let
1.1421 - val vc=length(#2(hd(p1')));
1.1422 - val cont =
1.1423 - (
1.1424 - if main_zero(mv_content(p1')) andalso
1.1425 - (main_zero(mv_content(p2'))) then
1.1426 - mv_correct((mv_gcd (mv_cut(mv_content(p1'))) (mv_cut(mv_content(p2')))),0)
1.1427 - else
1.1428 - mv_gcd (mv_content(p1')) (mv_content(p2'))
1.1429 - );
1.1430 - val p1= #1(mv_division(p1',mv_content(p1'),LEX_));
1.1431 - val p2= #1(mv_division(p2',mv_content(p2'),LEX_));
1.1432 - val gcd= Unsynchronized.ref [];
1.1433 - val candidate= Unsynchronized.ref [];
1.1434 - val interpolation_list= Unsynchronized.ref [];
1.1435 - val delta= Unsynchronized.ref [];
1.1436 - val p1r = Unsynchronized.ref [];
1.1437 - val p2r = Unsynchronized.ref [];
1.1438 - val p1r' = Unsynchronized.ref [];
1.1439 - val p2r' = Unsynchronized.ref [];
1.1440 - val factor= Unsynchronized.ref [];
1.1441 - val r= Unsynchronized.ref 0;
1.1442 - val gcd_r= Unsynchronized.ref [];
1.1443 - val d= Unsynchronized.ref 0;
1.1444 - val exit= Unsynchronized.ref 0;
1.1445 - val current_degree= Unsynchronized.ref 99999; (*. FIXME: unlimited ! .*)
1.1446 - in
1.1447 - (
1.1448 - if vc<2 then (* areUnivariate(p1',p2') *)
1.1449 - (
1.1450 - gcd:=uv_gcd (mv_shorten(p1',LEX_)) (mv_shorten(p2',LEX_))
1.1451 - )
1.1452 - else
1.1453 - (
1.1454 - while !exit=0 do
1.1455 - (
1.1456 - r:=(!r)+1;
1.1457 - p1r := mv_lc(p1,LEX_);
1.1458 - p2r := mv_lc(p2,LEX_);
1.1459 - if main_zero(!p1r) andalso
1.1460 - main_zero(!p2r)
1.1461 - then
1.1462 - (
1.1463 - delta := mv_correct((mv_gcd (mv_cut (!p1r)) (mv_cut (!p2r))),0)
1.1464 - )
1.1465 - else
1.1466 - (
1.1467 - delta := mv_gcd (!p1r) (!p2r)
1.1468 - );
1.1469 - (*if mv_shorten(mv_subs(!p1r,!r),LEX_)=[] andalso
1.1470 - mv_shorten(mv_subs(!p2r,!r),LEX_)=[] *)
1.1471 - if mv_lc2(mv_shorten(mv_subs(!p1r,!r),LEX_),LEX_)=0 andalso
1.1472 - mv_lc2(mv_shorten(mv_subs(!p2r,!r),LEX_),LEX_)=0
1.1473 - then
1.1474 - (
1.1475 - )
1.1476 - else
1.1477 - (
1.1478 - gcd_r:=mv_shorten(mv_gcd (mv_shorten(mv_subs(p1,!r),LEX_))
1.1479 - (mv_shorten(mv_subs(p2,!r),LEX_)) ,LEX_);
1.1480 - gcd_r:= #1(mv_division(mv_mul(mv_correct(mv_subs(!delta,!r),0),!gcd_r,LEX_),
1.1481 - mv_correct(mv_lc(!gcd_r,LEX_),0),LEX_));
1.1482 - d:=mv_deg2(!gcd_r); (* deg(gcd_r,z) *)
1.1483 - if (!d < !current_degree) then
1.1484 - (
1.1485 - current_degree:= !d;
1.1486 - interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
1.1487 - )
1.1488 - else
1.1489 - (
1.1490 - if (!d = !current_degree) then
1.1491 - (
1.1492 - interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
1.1493 - )
1.1494 - else ()
1.1495 - )
1.1496 - );
1.1497 - if length(!interpolation_list)> uv_mod_min(mv_deg(p1),mv_deg(p2)) then
1.1498 - (
1.1499 - candidate := mv_newton(rev(!interpolation_list));
1.1500 - if !candidate=[] then ()
1.1501 - else
1.1502 - (
1.1503 - candidate:=mv_pp(!candidate);
1.1504 - if mv_divides(!candidate,p1) andalso mv_divides(!candidate,p2) then
1.1505 - (
1.1506 - gcd:= mv_mul(!candidate,cont,LEX_);
1.1507 - exit:=1
1.1508 - )
1.1509 - else ()
1.1510 - );
1.1511 - interpolation_list:=[mv_correct(!gcd_r,0)]
1.1512 - )
1.1513 - else ()
1.1514 - )
1.1515 - );
1.1516 - (!gcd):mv_poly
1.1517 - )
1.1518 - end;
1.1519 -
1.1520 -
1.1521 -(*. calculates the least common divisor of two polynomials .*)
1.1522 -fun mv_lcm (p1:mv_poly) (p2:mv_poly) :mv_poly =
1.1523 - (
1.1524 - #1(mv_division(mv_mul(p1,p2,LEX_),mv_gcd p1 p2,LEX_))
1.1525 - );
1.1526 -
1.1527 -(* gets the variables (strings) of a term *)
1.1528 -
1.1529 -fun get_vars(term1) = (map free2str) (vars term1); (*["a","b","c"]; *)
1.1530 -
1.1531 -(*. counts the negative coefficents in a polynomial .*)
1.1532 -fun count_neg ([]:mv_poly) = 0
1.1533 - | count_neg ((c,e)::xs) = if c<0 then 1+count_neg xs
1.1534 - else count_neg xs;
1.1535 -
1.1536 -(*. help function for is_polynomial
1.1537 - checks the order of the operators .*)
1.1538 -fun test_polynomial (Const ("uminus",_) $ Free (str,_)) _ = true (*WN.13.3.03*)
1.1539 - | test_polynomial (t as Free(str,_)) v = true
1.1540 - | test_polynomial (t as Const ("Groups.times_class.times",_) $ t1 $ t2) v = if v="^" then false
1.1541 - else (test_polynomial t1 "*") andalso (test_polynomial t2 "*")
1.1542 - | test_polynomial (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1.1543 - else (test_polynomial t1 " ") andalso (test_polynomial t2 " ")
1.1544 - | test_polynomial (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_polynomial t1 "^") andalso (test_polynomial t2 "^")
1.1545 - | test_polynomial _ v = false;
1.1546 -
1.1547 -(*. tests if a term is a polynomial .*)
1.1548 -fun is_polynomial t = test_polynomial t " ";
1.1549 -
1.1550 -(*. help function for is_expanded
1.1551 - checks the order of the operators .*)
1.1552 -fun test_exp (t as Free(str,_)) v = true
1.1553 - | test_exp (t as Const ("Groups.times_class.times",_) $ t1 $ t2) v = if v="^" then false
1.1554 - else (test_exp t1 "*") andalso (test_exp t2 "*")
1.1555 - | test_exp (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1.1556 - else (test_exp t1 " ") andalso (test_exp t2 " ")
1.1557 - | test_exp (t as Const ("Groups.minus_class.minus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1.1558 - else (test_exp t1 " ") andalso (test_exp t2 " ")
1.1559 - | test_exp (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_exp t1 "^") andalso (test_exp t2 "^")
1.1560 - | test_exp _ v = false;
1.1561 -
1.1562 -
1.1563 -(*. help function for check_coeff:
1.1564 - converts the term to a list of coefficients .*)
1.1565 -fun term2coef' (t as Free(str,_(*typ*))) v :mv_poly option =
1.1566 - let
1.1567 - val x= Unsynchronized.ref NONE;
1.1568 - val len= Unsynchronized.ref 0;
1.1569 - val vl= Unsynchronized.ref [];
1.1570 - val vh= Unsynchronized.ref [];
1.1571 - val i= Unsynchronized.ref 0;
1.1572 - in
1.1573 - if is_numeral str then
1.1574 - (
1.1575 - SOME [(((the o int_of_str) str),mv_null2(v))] handle _ => NONE
1.1576 - )
1.1577 - else (* variable *)
1.1578 - (
1.1579 - len:=length(v);
1.1580 - vh:=v;
1.1581 - while ((!len)>(!i)) do
1.1582 - (
1.1583 - if str=hd((!vh)) then
1.1584 - (
1.1585 - vl:=1::(!vl)
1.1586 - )
1.1587 - else
1.1588 - (
1.1589 - vl:=0::(!vl)
1.1590 - );
1.1591 - vh:=tl(!vh);
1.1592 - i:=(!i)+1
1.1593 - );
1.1594 - SOME [(1,rev(!vl))] handle _ => NONE
1.1595 - )
1.1596 - end
1.1597 - | term2coef' (Const ("Groups.times_class.times",_) $ t1 $ t2) v :mv_poly option=
1.1598 - let
1.1599 - val t1pp= Unsynchronized.ref [];
1.1600 - val t2pp= Unsynchronized.ref [];
1.1601 - val t1c= Unsynchronized.ref 0;
1.1602 - val t2c= Unsynchronized.ref 0;
1.1603 - in
1.1604 - (
1.1605 - t1pp:=(#2(hd(the(term2coef' t1 v))));
1.1606 - t2pp:=(#2(hd(the(term2coef' t2 v))));
1.1607 - t1c:=(#1(hd(the(term2coef' t1 v))));
1.1608 - t2c:=(#1(hd(the(term2coef' t2 v))));
1.1609 -
1.1610 - SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] handle _ => NONE
1.1611 -
1.1612 - )
1.1613 - end
1.1614 - | term2coef' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ (t2 as Free (str2,_))) v :mv_poly option=
1.1615 - let
1.1616 - val x= Unsynchronized.ref NONE;
1.1617 - val len= Unsynchronized.ref 0;
1.1618 - val vl= Unsynchronized.ref [];
1.1619 - val vh= Unsynchronized.ref [];
1.1620 - val vtemp= Unsynchronized.ref [];
1.1621 - val i= Unsynchronized.ref 0;
1.1622 - in
1.1623 - (
1.1624 - if (not o is_numeral) str1 andalso is_numeral str2 then
1.1625 - (
1.1626 - len:=length(v);
1.1627 - vh:=v;
1.1628 -
1.1629 - while ((!len)>(!i)) do
1.1630 - (
1.1631 - if str1=hd((!vh)) then
1.1632 - (
1.1633 - vl:=((the o int_of_str) str2)::(!vl)
1.1634 - )
1.1635 - else
1.1636 - (
1.1637 - vl:=0::(!vl)
1.1638 - );
1.1639 - vh:=tl(!vh);
1.1640 - i:=(!i)+1
1.1641 - );
1.1642 - SOME [(1,rev(!vl))] handle _ => NONE
1.1643 - )
1.1644 - else error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term")
1.1645 - )
1.1646 - end
1.1647 - | term2coef' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option=
1.1648 - (
1.1649 - SOME ((the(term2coef' t1 v)) @ (the(term2coef' t2 v))) handle _ => NONE
1.1650 - )
1.1651 - | term2coef' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option=
1.1652 - (
1.1653 - SOME ((the(term2coef' t1 v)) @ mv_skalar_mul((the(term2coef' t2 v)),1)) handle _ => NONE
1.1654 - )
1.1655 - | term2coef' (term) v = error ("RATIONALS_TERM2COEF_EXCEPTION 2: Invalid term");
1.1656 -
1.1657 -(*. checks if all coefficients of a polynomial are positiv (except the first) .*)
1.1658 -fun check_coeff t = (* erste Koeffizient kann <0 sein !!! *)
1.1659 - if count_neg(tl(the(term2coef' t (get_vars(t)))))=0 then true
1.1660 - else false;
1.1661 -
1.1662 -(*WN.7.3.03 Hilfsfunktion f"ur term2poly'*)
1.1663 -fun mk_monom v' p vs =
1.1664 - let fun conv p (v: string) = if v'= v then p else 0
1.1665 - in map (conv p) vs end;
1.1666 -(* mk_monom "y" 5 ["a","b","x","y","z"];
1.1667 -val it = [0,0,0,5,0] : int list*)
1.1668 -
1.1669 -(*. this function converts the term representation into the internal representation mv_poly .*)
1.1670 -fun term2poly' (Const ("uminus",_) $ Free (str,_)) v = (*WN.7.3.03*)
1.1671 - if is_numeral str
1.1672 - then SOME [((the o int_of_str) ("-"^str), mk_monom "#" 0 v)]
1.1673 - else SOME [(~1, mk_monom str 1 v)]
1.1674 -
1.1675 - | term2poly' (Free(str,_)) v :mv_poly option =
1.1676 - let
1.1677 - val x= Unsynchronized.ref NONE;
1.1678 - val len= Unsynchronized.ref 0;
1.1679 - val vl= Unsynchronized.ref [];
1.1680 - val vh= Unsynchronized.ref [];
1.1681 - val i= Unsynchronized.ref 0;
1.1682 - in
1.1683 - if is_numeral str then
1.1684 - (
1.1685 - SOME [(((the o int_of_str) str),mv_null2 v)] handle _ => NONE
1.1686 - )
1.1687 - else (* variable *)
1.1688 - (
1.1689 - len:=length v;
1.1690 - vh:= v;
1.1691 - while ((!len)>(!i)) do
1.1692 - (
1.1693 - if str=hd((!vh)) then
1.1694 - (
1.1695 - vl:=1::(!vl)
1.1696 - )
1.1697 - else
1.1698 - (
1.1699 - vl:=0::(!vl)
1.1700 - );
1.1701 - vh:=tl(!vh);
1.1702 - i:=(!i)+1
1.1703 - );
1.1704 - SOME [(1,rev(!vl))] handle _ => NONE
1.1705 - )
1.1706 - end
1.1707 - | term2poly' (Const ("Groups.times_class.times",_) $ t1 $ t2) v :mv_poly option=
1.1708 - let
1.1709 - val t1pp= Unsynchronized.ref [];
1.1710 - val t2pp= Unsynchronized.ref [];
1.1711 - val t1c= Unsynchronized.ref 0;
1.1712 - val t2c= Unsynchronized.ref 0;
1.1713 - in
1.1714 - (
1.1715 - t1pp:=(#2(hd(the(term2poly' t1 v))));
1.1716 - t2pp:=(#2(hd(the(term2poly' t2 v))));
1.1717 - t1c:=(#1(hd(the(term2poly' t1 v))));
1.1718 - t2c:=(#1(hd(the(term2poly' t2 v))));
1.1719 -
1.1720 - SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )]
1.1721 - handle _ => NONE
1.1722 -
1.1723 - )
1.1724 - end
1.1725 - | term2poly' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $
1.1726 - (t2 as Free (str2,_))) v :mv_poly option=
1.1727 - let
1.1728 - val x= Unsynchronized.ref NONE;
1.1729 - val len= Unsynchronized.ref 0;
1.1730 - val vl= Unsynchronized.ref [];
1.1731 - val vh= Unsynchronized.ref [];
1.1732 - val vtemp= Unsynchronized.ref [];
1.1733 - val i= Unsynchronized.ref 0;
1.1734 - in
1.1735 - (
1.1736 - if (not o is_numeral) str1 andalso is_numeral str2 then
1.1737 - (
1.1738 - len:=length(v);
1.1739 - vh:=v;
1.1740 -
1.1741 - while ((!len)>(!i)) do
1.1742 - (
1.1743 - if str1=hd((!vh)) then
1.1744 - (
1.1745 - vl:=((the o int_of_str) str2)::(!vl)
1.1746 - )
1.1747 - else
1.1748 - (
1.1749 - vl:=0::(!vl)
1.1750 - );
1.1751 - vh:=tl(!vh);
1.1752 - i:=(!i)+1
1.1753 - );
1.1754 - SOME [(1,rev(!vl))] handle _ => NONE
1.1755 - )
1.1756 - else error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term")
1.1757 - )
1.1758 - end
1.1759 - | term2poly' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option =
1.1760 - (
1.1761 - SOME ((the(term2poly' t1 v)) @ (the(term2poly' t2 v))) handle _ => NONE
1.1762 - )
1.1763 - | term2poly' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option =
1.1764 - (
1.1765 - SOME ((the(term2poly' t1 v)) @ mv_skalar_mul((the(term2poly' t2 v)),~1)) handle _ => NONE
1.1766 - )
1.1767 - | term2poly' (term) v = error ("RATIONALS_TERM2POLY_EXCEPTION 2: Invalid term");
1.1768 -
1.1769 -(*. translates an Isabelle term into internal representation.
1.1770 - term2poly
1.1771 - fn : term -> (*normalform [2] *)
1.1772 - string list -> (*for ...!!! BITTE DIE ERKLÄRUNG,
1.1773 - DIE DU MIR LETZTES MAL GEGEBEN HAST*)
1.1774 - mv_monom list (*internal representation *)
1.1775 - option (*the translation may fail with NONE*)
1.1776 -.*)
1.1777 -fun term2poly (t:term) v =
1.1778 - if is_polynomial t then term2poly' t v
1.1779 - else error ("term2poly: invalid = "^(term2str t));
1.1780 -
1.1781 -(*. same as term2poly with automatic detection of the variables .*)
1.1782 -fun term2polyx t = term2poly t (((map free2str) o vars) t);
1.1783 -
1.1784 -(*. checks if the term is in expanded polynomial form and converts it into the internal representation .*)
1.1785 -fun expanded2poly (t:term) v =
1.1786 - (*if is_expanded t then*) term2poly' t v
1.1787 - (*else error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*);
1.1788 -
1.1789 -(*. same as expanded2poly with automatic detection of the variables .*)
1.1790 -fun expanded2polyx t = expanded2poly t (((map free2str) o vars) t);
1.1791 -
1.1792 -(*. converts a powerproduct into term representation .*)
1.1793 -fun powerproduct2term(xs,v) =
1.1794 - let
1.1795 - val xss= Unsynchronized.ref xs;
1.1796 - val vv= Unsynchronized.ref v;
1.1797 - in
1.1798 - (
1.1799 - while hd(!xss)=0 do
1.1800 - (
1.1801 - xss:=tl(!xss);
1.1802 - vv:=tl(!vv)
1.1803 - );
1.1804 -
1.1805 - if list_is_null(tl(!xss)) then
1.1806 - (
1.1807 - if hd(!xss)=1 then Free(hd(!vv), HOLogic.realT)
1.1808 - else
1.1809 - (
1.1810 - Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1811 - Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT)
1.1812 - )
1.1813 - )
1.1814 - else
1.1815 - (
1.1816 - if hd(!xss)=1 then
1.1817 - (
1.1818 - Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1819 - Free(hd(!vv), HOLogic.realT) $
1.1820 - powerproduct2term(tl(!xss),tl(!vv))
1.1821 - )
1.1822 - else
1.1823 - (
1.1824 - Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1825 - (
1.1826 - Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1827 - Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT)
1.1828 - ) $
1.1829 - powerproduct2term(tl(!xss),tl(!vv))
1.1830 - )
1.1831 - )
1.1832 - )
1.1833 - end;
1.1834 -
1.1835 -(*. converts a monom into term representation .*)
1.1836 -(*fun monom2term ((c,e):mv_monom, v:string list) =
1.1837 - if c=0 then Free(str_of_int 0,HOLogic.realT)
1.1838 - else
1.1839 - (
1.1840 - if list_is_null(e) then
1.1841 - (
1.1842 - Free(str_of_int c,HOLogic.realT)
1.1843 - )
1.1844 - else
1.1845 - (
1.1846 - if c=1 then
1.1847 - (
1.1848 - powerproduct2term(e,v)
1.1849 - )
1.1850 - else
1.1851 - (
1.1852 - Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1853 - Free(str_of_int c,HOLogic.realT) $
1.1854 - powerproduct2term(e,v)
1.1855 - )
1.1856 - )
1.1857 - );*)
1.1858 -
1.1859 -
1.1860 -(*fun monom2term ((i, is):mv_monom, v) =
1.1861 - if list_is_null is
1.1862 - then
1.1863 - if i >= 0
1.1864 - then Free (str_of_int i, HOLogic.realT)
1.1865 - else Const ("uminus", HOLogic.realT --> HOLogic.realT) $
1.1866 - Free ((str_of_int o abs) i, HOLogic.realT)
1.1867 - else
1.1868 - if i > 0
1.1869 - then Const ("Groups.times_class.times", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
1.1870 - (Free (str_of_int i, HOLogic.realT)) $
1.1871 - powerproduct2term(is, v)
1.1872 - else Const ("Groups.times_class.times", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
1.1873 - (Const ("uminus", HOLogic.realT --> HOLogic.realT) $
1.1874 - Free ((str_of_int o abs) i, HOLogic.realT)) $
1.1875 - powerproduct2term(is, vs);---------------------------*)
1.1876 -fun monom2term ((i, is) : mv_monom, vs) =
1.1877 - if list_is_null is
1.1878 - then Free (str_of_int i, HOLogic.realT)
1.1879 - else if i = 1
1.1880 - then powerproduct2term (is, vs)
1.1881 - else Const ("Groups.times_class.times", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
1.1882 - (Free (str_of_int i, HOLogic.realT)) $
1.1883 - powerproduct2term (is, vs);
1.1884 -
1.1885 -(*. converts the internal polynomial representation into an Isabelle term.*)
1.1886 -fun poly2term' ([] : mv_poly, vs) = Free(str_of_int 0, HOLogic.realT)
1.1887 - | poly2term' ([(c, e) : mv_monom], vs) = monom2term ((c, e), vs)
1.1888 - | poly2term' ((c, e) :: ces, vs) =
1.1889 - Const("Groups.plus_class.plus", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
1.1890 - poly2term (ces, vs) $ monom2term ((c, e), vs)
1.1891 -and poly2term (xs, vs) = poly2term' (rev (sort (mv_geq LEX_) (xs)), vs);
1.1892 -
1.1893 -
1.1894 -(*. converts a monom into term representation .*)
1.1895 -(*. ignores the sign of the coefficients => use only for exp-poly functions .*)
1.1896 -fun monom2term2((c,e):mv_monom, v:string list) =
1.1897 - if c=0 then Free(str_of_int 0,HOLogic.realT)
1.1898 - else
1.1899 - (
1.1900 - if list_is_null(e) then
1.1901 - (
1.1902 - Free(str_of_int (abs(c)),HOLogic.realT)
1.1903 - )
1.1904 - else
1.1905 - (
1.1906 - if abs(c)=1 then
1.1907 - (
1.1908 - powerproduct2term(e,v)
1.1909 - )
1.1910 - else
1.1911 - (
1.1912 - Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1913 - Free(str_of_int (abs(c)),HOLogic.realT) $
1.1914 - powerproduct2term(e,v)
1.1915 - )
1.1916 - )
1.1917 - );
1.1918 -
1.1919 -(*. converts the expanded polynomial representation into the term representation .*)
1.1920 -fun exp2term' ([]:mv_poly,vars) = Free(str_of_int 0,HOLogic.realT)
1.1921 - | exp2term' ([(c,e)],vars) = monom2term((c,e),vars)
1.1922 - | exp2term' ((c1,e1)::others,vars) =
1.1923 - if c1<0 then
1.1924 - Const("Groups.minus_class.minus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1925 - exp2term'(others,vars) $
1.1926 - (
1.1927 - monom2term2((c1,e1),vars)
1.1928 - )
1.1929 - else
1.1930 - Const("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1931 - exp2term'(others,vars) $
1.1932 - (
1.1933 - monom2term2((c1,e1),vars)
1.1934 - );
1.1935 -
1.1936 -(*. sorts the powerproduct by lexicographic termorder and converts them into
1.1937 - a term in polynomial representation .*)
1.1938 -fun poly2expanded (xs,vars) = exp2term'(rev(sort (mv_geq LEX_) (xs)),vars);
1.1939 -
1.1940 -(*. converts a polynomial into expanded form .*)
1.1941 -fun polynomial2expanded t =
1.1942 - (let
1.1943 - val vars=(((map free2str) o vars) t);
1.1944 - in
1.1945 - SOME (poly2expanded (the (term2poly t vars), vars))
1.1946 - end) handle _ => NONE;
1.1947 -
1.1948 -(*. converts a polynomial into polynomial form .*)
1.1949 -fun expanded2polynomial t =
1.1950 - (let
1.1951 - val vars=(((map free2str) o vars) t);
1.1952 - in
1.1953 - SOME (poly2term (the (expanded2poly t vars), vars))
1.1954 - end) handle _ => NONE;
1.1955 -
1.1956 -
1.1957 -(*. calculates the greatest common divisor of numerator and denominator and seperates it from each .*)
1.1958 -fun step_cancel (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =
1.1959 - let
1.1960 - val p1' = Unsynchronized.ref [];
1.1961 - val p2' = Unsynchronized.ref [];
1.1962 - val p3 = Unsynchronized.ref []
1.1963 - val vars = rev(get_vars(p1) union get_vars(p2));
1.1964 - in
1.1965 - (
1.1966 - p1':= sort (mv_geq LEX_) (the (term2poly p1 vars ));
1.1967 - p2':= sort (mv_geq LEX_) (the (term2poly p2 vars ));
1.1968 - p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1.1969 - if (!p3)=[(1,mv_null2(vars))] then
1.1970 - (
1.1971 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2
1.1972 - )
1.1973 - else
1.1974 - (
1.1975 -
1.1976 - p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1.1977 - p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1.1978 -
1.1979 - if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1.1980 - (
1.1981 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.1982 - $
1.1983 - (
1.1984 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1985 - poly2term(!p1',vars) $
1.1986 - poly2term(!p3,vars)
1.1987 - )
1.1988 - $
1.1989 - (
1.1990 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.1991 - poly2term(!p2',vars) $
1.1992 - poly2term(!p3,vars)
1.1993 - )
1.1994 - )
1.1995 - else
1.1996 - (
1.1997 - p1':=mv_skalar_mul(!p1',~1);
1.1998 - p2':=mv_skalar_mul(!p2',~1);
1.1999 - p3:=mv_skalar_mul(!p3,~1);
1.2000 - (
1.2001 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2002 - $
1.2003 - (
1.2004 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2005 - poly2term(!p1',vars) $
1.2006 - poly2term(!p3,vars)
1.2007 - )
1.2008 - $
1.2009 - (
1.2010 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2011 - poly2term(!p2',vars) $
1.2012 - poly2term(!p3,vars)
1.2013 - )
1.2014 - )
1.2015 - )
1.2016 - )
1.2017 - )
1.2018 - end
1.2019 -| step_cancel _ = error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction");
1.2020 -
1.2021 -(*. calculates the greatest common divisor of numerator and denominator and divides each through it .*)
1.2022 -fun direct_cancel (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =
1.2023 - let
1.2024 - val p1' = Unsynchronized.ref [];
1.2025 - val p2' = Unsynchronized.ref [];
1.2026 - val p3 = Unsynchronized.ref []
1.2027 - val vars = rev(get_vars(p1) union get_vars(p2));
1.2028 - in
1.2029 - (
1.2030 - p1':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p1 vars )),LEX_));
1.2031 - p2':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p2 vars )),LEX_));
1.2032 - p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1.2033 -
1.2034 - if (!p3)=[(1,mv_null2(vars))] then
1.2035 - (
1.2036 - (Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
1.2037 - )
1.2038 - else
1.2039 - (
1.2040 - p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1.2041 - p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1.2042 - if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1.2043 - (
1.2044 - (
1.2045 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2046 - $
1.2047 - (
1.2048 - poly2term((!p1'),vars)
1.2049 - )
1.2050 - $
1.2051 - (
1.2052 - poly2term((!p2'),vars)
1.2053 - )
1.2054 - )
1.2055 - ,
1.2056 - if mv_grad(!p3)>0 then
1.2057 - [
1.2058 - (
1.2059 - Const ("HOL.Not",[bool]--->bool) $
1.2060 - (
1.2061 - Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $
1.2062 - poly2term((!p3),vars) $
1.2063 - Free("0",HOLogic.realT)
1.2064 - )
1.2065 - )
1.2066 - ]
1.2067 - else
1.2068 - []
1.2069 - )
1.2070 - else
1.2071 - (
1.2072 - p1':=mv_skalar_mul(!p1',~1);
1.2073 - p2':=mv_skalar_mul(!p2',~1);
1.2074 - if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);
1.2075 - (
1.2076 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2077 - $
1.2078 - (
1.2079 - poly2term((!p1'),vars)
1.2080 - )
1.2081 - $
1.2082 - (
1.2083 - poly2term((!p2'),vars)
1.2084 - )
1.2085 - ,
1.2086 - if mv_grad(!p3)>0 then
1.2087 - [
1.2088 - (
1.2089 - Const ("HOL.Not",[bool]--->bool) $
1.2090 - (
1.2091 - Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $
1.2092 - poly2term((!p3),vars) $
1.2093 - Free("0",HOLogic.realT)
1.2094 - )
1.2095 - )
1.2096 - ]
1.2097 - else
1.2098 - []
1.2099 - )
1.2100 - )
1.2101 - )
1.2102 - )
1.2103 - end
1.2104 - | direct_cancel _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");
1.2105 -
1.2106 -(*. same es direct_cancel, this time for expanded forms (input+output).*)
1.2107 -fun direct_cancel_expanded (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =
1.2108 - let
1.2109 - val p1' = Unsynchronized.ref [];
1.2110 - val p2' = Unsynchronized.ref [];
1.2111 - val p3 = Unsynchronized.ref []
1.2112 - val vars = rev(get_vars(p1) union get_vars(p2));
1.2113 - in
1.2114 - (
1.2115 - p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_));
1.2116 - p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_));
1.2117 - p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1.2118 -
1.2119 - if (!p3)=[(1,mv_null2(vars))] then
1.2120 - (
1.2121 - (Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
1.2122 - )
1.2123 - else
1.2124 - (
1.2125 - p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1.2126 - p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1.2127 - if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1.2128 - (
1.2129 - (
1.2130 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2131 - $
1.2132 - (
1.2133 - poly2expanded((!p1'),vars)
1.2134 - )
1.2135 - $
1.2136 - (
1.2137 - poly2expanded((!p2'),vars)
1.2138 - )
1.2139 - )
1.2140 - ,
1.2141 - if mv_grad(!p3)>0 then
1.2142 - [
1.2143 - (
1.2144 - Const ("HOL.Not",[bool]--->bool) $
1.2145 - (
1.2146 - Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $
1.2147 - poly2expanded((!p3),vars) $
1.2148 - Free("0",HOLogic.realT)
1.2149 - )
1.2150 - )
1.2151 - ]
1.2152 - else
1.2153 - []
1.2154 - )
1.2155 - else
1.2156 - (
1.2157 - p1':=mv_skalar_mul(!p1',~1);
1.2158 - p2':=mv_skalar_mul(!p2',~1);
1.2159 - if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);
1.2160 - (
1.2161 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2162 - $
1.2163 - (
1.2164 - poly2expanded((!p1'),vars)
1.2165 - )
1.2166 - $
1.2167 - (
1.2168 - poly2expanded((!p2'),vars)
1.2169 - )
1.2170 - ,
1.2171 - if mv_grad(!p3)>0 then
1.2172 - [
1.2173 - (
1.2174 - Const ("HOL.Not",[bool]--->bool) $
1.2175 - (
1.2176 - Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $
1.2177 - poly2expanded((!p3),vars) $
1.2178 - Free("0",HOLogic.realT)
1.2179 - )
1.2180 - )
1.2181 - ]
1.2182 - else
1.2183 - []
1.2184 - )
1.2185 - )
1.2186 - )
1.2187 - )
1.2188 - end
1.2189 - | direct_cancel_expanded _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");
1.2190 -
1.2191 -
1.2192 -(*. adds two fractions .*)
1.2193 -fun add_fract ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =
1.2194 - let
1.2195 - val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
1.2196 - val t11'= Unsynchronized.ref (the(term2poly t11 vars));
1.2197 -(* stopped Test_Isac.thy ...
1.2198 -val _= tracing"### add_fract: done t11"
1.2199 -*)
1.2200 - val t12'= Unsynchronized.ref (the(term2poly t12 vars));
1.2201 -(* stopped Test_Isac.thy ...
1.2202 -val _= tracing"### add_fract: done t12"
1.2203 -*)
1.2204 - val t21'= Unsynchronized.ref (the(term2poly t21 vars));
1.2205 -(* stopped Test_Isac.thy ...
1.2206 -val _= tracing"### add_fract: done t21"
1.2207 -*)
1.2208 - val t22'= Unsynchronized.ref (the(term2poly t22 vars));
1.2209 -(* stopped Test_Isac.thy ...
1.2210 -val _= tracing"### add_fract: done t22"
1.2211 -*)
1.2212 - val den= Unsynchronized.ref [];
1.2213 - val nom= Unsynchronized.ref [];
1.2214 - val m1= Unsynchronized.ref [];
1.2215 - val m2= Unsynchronized.ref [];
1.2216 - in
1.2217 -
1.2218 - (
1.2219 - den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
1.2220 -tracing"### add_fract: done sort mv_lcm";
1.2221 - m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
1.2222 -tracing"### add_fract: done sort mv_division t12";
1.2223 - m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
1.2224 -tracing"### add_fract: done sort mv_division t22";
1.2225 - nom :=sort (mv_geq LEX_)
1.2226 - (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),
1.2227 - mv_mul(!t21',!m2,LEX_),
1.2228 - LEX_),
1.2229 - LEX_));
1.2230 -tracing"### add_fract: done sort mv_add";
1.2231 - (
1.2232 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2233 - $
1.2234 - (
1.2235 - poly2term((!nom),vars)
1.2236 - )
1.2237 - $
1.2238 - (
1.2239 - poly2term((!den),vars)
1.2240 - )
1.2241 - )
1.2242 - )
1.2243 - end
1.2244 - | add_fract (_,_) = error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call");
1.2245 -
1.2246 -(*. adds two expanded fractions .*)
1.2247 -fun add_fract_exp ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =
1.2248 - let
1.2249 - val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
1.2250 - val t11'= Unsynchronized.ref (the(expanded2poly t11 vars));
1.2251 - val t12'= Unsynchronized.ref (the(expanded2poly t12 vars));
1.2252 - val t21'= Unsynchronized.ref (the(expanded2poly t21 vars));
1.2253 - val t22'= Unsynchronized.ref (the(expanded2poly t22 vars));
1.2254 - val den= Unsynchronized.ref [];
1.2255 - val nom= Unsynchronized.ref [];
1.2256 - val m1= Unsynchronized.ref [];
1.2257 - val m2= Unsynchronized.ref [];
1.2258 - in
1.2259 -
1.2260 - (
1.2261 - den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
1.2262 - m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
1.2263 - m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
1.2264 - nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_));
1.2265 - (
1.2266 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2267 - $
1.2268 - (
1.2269 - poly2expanded((!nom),vars)
1.2270 - )
1.2271 - $
1.2272 - (
1.2273 - poly2expanded((!den),vars)
1.2274 - )
1.2275 - )
1.2276 - )
1.2277 - end
1.2278 - | add_fract_exp (_,_) = error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call");
1.2279 -
1.2280 -(*. adds a list of terms .*)
1.2281 -fun add_list_of_fractions []= (Free("0",HOLogic.realT),[])
1.2282 - | add_list_of_fractions [x]= direct_cancel x
1.2283 - | add_list_of_fractions (x::y::xs) =
1.2284 - let
1.2285 - val (t1a,rest1)=direct_cancel(x);
1.2286 -val _= tracing"### add_list_of_fractions xs: has done direct_cancel(x)";
1.2287 - val (t2a,rest2)=direct_cancel(y);
1.2288 -val _= tracing"### add_list_of_fractions xs: has done direct_cancel(y)";
1.2289 - val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs));
1.2290 -val _= tracing"### add_list_of_fractions xs: has done add_list_of_fraction xs";
1.2291 - val (t4a,rest4)=direct_cancel(t3a);
1.2292 -val _= tracing"### add_list_of_fractions xs: has done direct_cancel(t3a)";
1.2293 - val rest=rest1 union rest2 union rest3 union rest4;
1.2294 - in
1.2295 - (tracing"### add_list_of_fractions in";
1.2296 - (
1.2297 - (t4a,rest)
1.2298 - )
1.2299 - )
1.2300 - end;
1.2301 -
1.2302 -(*. adds a list of expanded terms .*)
1.2303 -fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])
1.2304 - | add_list_of_fractions_exp [x]= direct_cancel_expanded x
1.2305 - | add_list_of_fractions_exp (x::y::xs) =
1.2306 - let
1.2307 - val (t1a,rest1)=direct_cancel_expanded(x);
1.2308 - val (t2a,rest2)=direct_cancel_expanded(y);
1.2309 - val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));
1.2310 - val (t4a,rest4)=direct_cancel_expanded(t3a);
1.2311 - val rest=rest1 union rest2 union rest3 union rest4;
1.2312 - in
1.2313 - (
1.2314 - (t4a,rest)
1.2315 - )
1.2316 - end;
1.2317 -
1.2318 -(*. calculates the lcm of a list of mv_poly .*)
1.2319 -fun calc_lcm ([x],var)= (x,var)
1.2320 - | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
1.2321 -
1.2322 -(*. converts a list of terms to a list of mv_poly .*)
1.2323 -fun t2d([],_)=[]
1.2324 - | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);
1.2325 -
1.2326 -(*. same as t2d, this time for expanded forms .*)
1.2327 -fun t2d_exp([],_)=[]
1.2328 - | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
1.2329 -
1.2330 -(*. converts a list of fract terms to a list of their denominators .*)
1.2331 -fun termlist2denominators [] = ([],[])
1.2332 - | termlist2denominators xs =
1.2333 - let
1.2334 - val xxs= Unsynchronized.ref xs;
1.2335 - val var= Unsynchronized.ref [];
1.2336 - in
1.2337 - var:=[];
1.2338 - while length(!xxs)>0 do
1.2339 - (
1.2340 - let
1.2341 - val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);
1.2342 - in
1.2343 - (
1.2344 - xxs:=tl(!xxs);
1.2345 - var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
1.2346 - )
1.2347 - end
1.2348 - );
1.2349 - (t2d(xs,!var),!var)
1.2350 - end;
1.2351 -
1.2352 -(*. calculates the lcm of a list of mv_poly .*)
1.2353 -fun calc_lcm ([x],var)= (x,var)
1.2354 - | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
1.2355 -
1.2356 -(*. converts a list of terms to a list of mv_poly .*)
1.2357 -fun t2d([],_)=[]
1.2358 - | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);
1.2359 -
1.2360 -(*. same as t2d, this time for expanded forms .*)
1.2361 -fun t2d_exp([],_)=[]
1.2362 - | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
1.2363 -
1.2364 -(*. converts a list of fract terms to a list of their denominators .*)
1.2365 -fun termlist2denominators [] = ([],[])
1.2366 - | termlist2denominators xs =
1.2367 - let
1.2368 - val xxs= Unsynchronized.ref xs;
1.2369 - val var= Unsynchronized.ref [];
1.2370 - in
1.2371 - var:=[];
1.2372 - while length(!xxs)>0 do
1.2373 - (
1.2374 - let
1.2375 - val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);
1.2376 - in
1.2377 - (
1.2378 - xxs:=tl(!xxs);
1.2379 - var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
1.2380 - )
1.2381 - end
1.2382 - );
1.2383 - (t2d(xs,!var),!var)
1.2384 - end;
1.2385 -
1.2386 -(*. same as termlist2denminators, this time for expanded forms .*)
1.2387 -fun termlist2denominators_exp [] = ([],[])
1.2388 - | termlist2denominators_exp xs =
1.2389 - let
1.2390 - val xxs= Unsynchronized.ref xs;
1.2391 - val var= Unsynchronized.ref [];
1.2392 - in
1.2393 - var:=[];
1.2394 - while length(!xxs)>0 do
1.2395 - (
1.2396 - let
1.2397 - val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);
1.2398 - in
1.2399 - (
1.2400 - xxs:=tl(!xxs);
1.2401 - var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
1.2402 - )
1.2403 - end
1.2404 - );
1.2405 - (t2d_exp(xs,!var),!var)
1.2406 - end;
1.2407 -
1.2408 -(*. reduces all fractions to the least common denominator .*)
1.2409 -fun com_den(x::xs,denom,den,var)=
1.2410 - let
1.2411 - val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;
1.2412 - val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
1.2413 - val p3= #1(mv_division(denom,p2,LEX_));
1.2414 - val p1var=get_vars(p1');
1.2415 - in
1.2416 - if length(xs)>0 then
1.2417 - if p3=[(1,mv_null2(var))] then
1.2418 - (
1.2419 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2420 - $
1.2421 - (
1.2422 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2423 - $
1.2424 - poly2term(the (term2poly p1' p1var),p1var)
1.2425 - $
1.2426 - den
1.2427 - )
1.2428 - $
1.2429 - #1(com_den(xs,denom,den,var))
1.2430 - ,
1.2431 - []
1.2432 - )
1.2433 - else
1.2434 - (
1.2435 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2436 - $
1.2437 - (
1.2438 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2439 - $
1.2440 - (
1.2441 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2442 - poly2term(the (term2poly p1' p1var),p1var) $
1.2443 - poly2term(p3,var)
1.2444 - )
1.2445 - $
1.2446 - (
1.2447 - den
1.2448 - )
1.2449 - )
1.2450 - $
1.2451 - #1(com_den(xs,denom,den,var))
1.2452 - ,
1.2453 - []
1.2454 - )
1.2455 - else
1.2456 - if p3=[(1,mv_null2(var))] then
1.2457 - (
1.2458 - (
1.2459 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2460 - $
1.2461 - poly2term(the (term2poly p1' p1var),p1var)
1.2462 - $
1.2463 - den
1.2464 - )
1.2465 - ,
1.2466 - []
1.2467 - )
1.2468 - else
1.2469 - (
1.2470 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2471 - $
1.2472 - (
1.2473 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2474 - poly2term(the (term2poly p1' p1var),p1var) $
1.2475 - poly2term(p3,var)
1.2476 - )
1.2477 - $
1.2478 - den
1.2479 - ,
1.2480 - []
1.2481 - )
1.2482 - end;
1.2483 -
1.2484 -(*. same as com_den, this time for expanded forms .*)
1.2485 -fun com_den_exp(x::xs,denom,den,var)=
1.2486 - let
1.2487 - val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;
1.2488 - val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
1.2489 - val p3= #1(mv_division(denom,p2,LEX_));
1.2490 - val p1var=get_vars(p1');
1.2491 - in
1.2492 - if length(xs)>0 then
1.2493 - if p3=[(1,mv_null2(var))] then
1.2494 - (
1.2495 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2496 - $
1.2497 - (
1.2498 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2499 - $
1.2500 - poly2expanded(the(expanded2poly p1' p1var),p1var)
1.2501 - $
1.2502 - den
1.2503 - )
1.2504 - $
1.2505 - #1(com_den_exp(xs,denom,den,var))
1.2506 - ,
1.2507 - []
1.2508 - )
1.2509 - else
1.2510 - (
1.2511 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2512 - $
1.2513 - (
1.2514 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2515 - $
1.2516 - (
1.2517 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2518 - poly2expanded(the(expanded2poly p1' p1var),p1var) $
1.2519 - poly2expanded(p3,var)
1.2520 - )
1.2521 - $
1.2522 - (
1.2523 - den
1.2524 - )
1.2525 - )
1.2526 - $
1.2527 - #1(com_den_exp(xs,denom,den,var))
1.2528 - ,
1.2529 - []
1.2530 - )
1.2531 - else
1.2532 - if p3=[(1,mv_null2(var))] then
1.2533 - (
1.2534 - (
1.2535 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2536 - $
1.2537 - poly2expanded(the(expanded2poly p1' p1var),p1var)
1.2538 - $
1.2539 - den
1.2540 - )
1.2541 - ,
1.2542 - []
1.2543 - )
1.2544 - else
1.2545 - (
1.2546 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1.2547 - $
1.2548 - (
1.2549 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2550 - poly2expanded(the(expanded2poly p1' p1var),p1var) $
1.2551 - poly2expanded(p3,var)
1.2552 - )
1.2553 - $
1.2554 - den
1.2555 - ,
1.2556 - []
1.2557 - )
1.2558 - end;
1.2559 -
1.2560 -(* wird aktuell nicht mehr gebraucht, bei rückänderung schon
1.2561 --------------------------------------------------------------
1.2562 -(* WN0210???SK brauch ma des überhaupt *)
1.2563 -fun com_den2(x::xs,denom,den,var)=
1.2564 - let
1.2565 - val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;
1.2566 - val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
1.2567 - val p3= #1(mv_division(denom,p2,LEX_));
1.2568 - val p1var=get_vars(p1');
1.2569 - in
1.2570 - if length(xs)>0 then
1.2571 - if p3=[(1,mv_null2(var))] then
1.2572 - (
1.2573 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2574 - poly2term(the(term2poly p1' p1var),p1var) $
1.2575 - com_den2(xs,denom,den,var)
1.2576 - )
1.2577 - else
1.2578 - (
1.2579 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2580 - (
1.2581 - let
1.2582 - val p3'=poly2term(p3,var);
1.2583 - val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
1.2584 - in
1.2585 - poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
1.2586 - end
1.2587 - ) $
1.2588 - com_den2(xs,denom,den,var)
1.2589 - )
1.2590 - else
1.2591 - if p3=[(1,mv_null2(var))] then
1.2592 - (
1.2593 - poly2term(the(term2poly p1' p1var),p1var)
1.2594 - )
1.2595 - else
1.2596 - (
1.2597 - let
1.2598 - val p3'=poly2term(p3,var);
1.2599 - val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
1.2600 - in
1.2601 - poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
1.2602 - end
1.2603 - )
1.2604 - end;
1.2605 -
1.2606 -(* WN0210???SK brauch ma des überhaupt *)
1.2607 -fun com_den_exp2(x::xs,denom,den,var)=
1.2608 - let
1.2609 - val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;
1.2610 - val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
1.2611 - val p3= #1(mv_division(denom,p2,LEX_));
1.2612 - val p1var=get_vars p1';
1.2613 - in
1.2614 - if length(xs)>0 then
1.2615 - if p3=[(1,mv_null2(var))] then
1.2616 - (
1.2617 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2618 - poly2expanded(the (expanded2poly p1' p1var),p1var) $
1.2619 - com_den_exp2(xs,denom,den,var)
1.2620 - )
1.2621 - else
1.2622 - (
1.2623 - Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2624 - (
1.2625 - let
1.2626 - val p3'=poly2expanded(p3,var);
1.2627 - val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
1.2628 - in
1.2629 - poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
1.2630 - end
1.2631 - ) $
1.2632 - com_den_exp2(xs,denom,den,var)
1.2633 - )
1.2634 - else
1.2635 - if p3=[(1,mv_null2(var))] then
1.2636 - (
1.2637 - poly2expanded(the (expanded2poly p1' p1var),p1var)
1.2638 - )
1.2639 - else
1.2640 - (
1.2641 - let
1.2642 - val p3'=poly2expanded(p3,var);
1.2643 - val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
1.2644 - in
1.2645 - poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
1.2646 - end
1.2647 - )
1.2648 - end;
1.2649 ----------------------------------------------------------*)
1.2650 -
1.2651 -
1.2652 -(*. searches for an element y of a list ys, which has an gcd not 1 with x .*)
1.2653 -fun exists_gcd (x,[]) = false
1.2654 - | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then exists_gcd (x,ys)
1.2655 - else true;
1.2656 -
1.2657 -(*. divides each element of the list xs with y .*)
1.2658 -fun list_div ([],y) = []
1.2659 - | list_div (x::xs,y) =
1.2660 - let
1.2661 - val (d,r)=mv_division(x,y,LEX_);
1.2662 - in
1.2663 - if r=[] then
1.2664 - d::list_div(xs,y)
1.2665 - else x::list_div(xs,y)
1.2666 - end;
1.2667 -
1.2668 -(*. checks if x is in the list ys .*)
1.2669 -fun in_list (x,[]) = false
1.2670 - | in_list (x,y::ys) = if x=y then true
1.2671 - else in_list(x,ys);
1.2672 -
1.2673 -(*. deletes all equal elements of the list xs .*)
1.2674 -fun kill_equal [] = []
1.2675 - | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)
1.2676 - else x::kill_equal(xs);
1.2677 -
1.2678 -(*. searches for new factors .*)
1.2679 -fun new_factors [] = []
1.2680 - | new_factors (list:mv_poly list):mv_poly list =
1.2681 - let
1.2682 - val l = kill_equal list;
1.2683 - val len = length(l);
1.2684 - in
1.2685 - if len>=2 then
1.2686 - (
1.2687 - let
1.2688 - val x::y::xs=l;
1.2689 - val gcd=mv_gcd x y;
1.2690 - in
1.2691 - if gcd=[(1,mv_null2(#2(hd(x))))] then
1.2692 - (
1.2693 - if exists_gcd(x,xs) then new_factors (y::xs @ [x])
1.2694 - else x::new_factors(y::xs)
1.2695 - )
1.2696 - else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))
1.2697 - end
1.2698 - )
1.2699 - else
1.2700 - if len=1 then [hd(l)]
1.2701 - else []
1.2702 - end;
1.2703 -
1.2704 -(*. gets the factors of a list .*)
1.2705 -fun get_factors x = new_factors x;
1.2706 -
1.2707 -(*. multiplies the elements of the list .*)
1.2708 -fun multi_list [] = []
1.2709 - | multi_list (x::xs) = if xs=[] then x
1.2710 - else mv_mul(x,multi_list xs,LEX_);
1.2711 -
1.2712 -(*. makes a term out of the elements of the list (polynomial representation) .*)
1.2713 -fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)
1.2714 - | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)
1.2715 - else
1.2716 - (
1.2717 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2718 - poly2term(sort (mv_geq LEX_) (x),vars) $
1.2719 - make_term(xs,vars)
1.2720 - );
1.2721 -
1.2722 -(*. factorizes the denominator (polynomial representation) .*)
1.2723 -fun factorize_den (l,den,vars) =
1.2724 - let
1.2725 - val factor_list=kill_equal( (get_factors l));
1.2726 - val mlist=multi_list(factor_list);
1.2727 - val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
1.2728 - in
1.2729 - if rest=[] then
1.2730 - (
1.2731 - if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)
1.2732 - else make_term(last::factor_list,vars)
1.2733 - )
1.2734 - else error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")
1.2735 - end;
1.2736 -
1.2737 -(*. makes a term out of the elements of the list (expanded polynomial representation) .*)
1.2738 -fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)
1.2739 - | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)
1.2740 - else
1.2741 - (
1.2742 - Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2743 - poly2expanded(sort (mv_geq LEX_) (x),vars) $
1.2744 - make_exp(xs,vars)
1.2745 - );
1.2746 -
1.2747 -(*. factorizes the denominator (expanded polynomial representation) .*)
1.2748 -fun factorize_den_exp (l,den,vars) =
1.2749 - let
1.2750 - val factor_list=kill_equal( (get_factors l));
1.2751 - val mlist=multi_list(factor_list);
1.2752 - val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
1.2753 - in
1.2754 - if rest=[] then
1.2755 - (
1.2756 - if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)
1.2757 - else make_exp(last::factor_list,vars)
1.2758 - )
1.2759 - else error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")
1.2760 - end;
1.2761 -
1.2762 -(*. calculates the common denominator of all elements of the list and multiplies .*)
1.2763 -(*. the nominators and denominators with the correct factor .*)
1.2764 -(*. (polynomial representation) .*)
1.2765 -fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)
1.2766 - | step_add_list_of_fractions [x]= error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")
1.2767 - | step_add_list_of_fractions (xs) =
1.2768 - let
1.2769 - val den_list=termlist2denominators (xs); (* list of denominators *)
1.2770 - val (denom,var)=calc_lcm(den_list); (* common denominator *)
1.2771 - val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
1.2772 - in
1.2773 - com_den(xs,denom,den,var)
1.2774 - end;
1.2775 -
1.2776 -(*. calculates the common denominator of all elements of the list and multiplies .*)
1.2777 -(*. the nominators and denominators with the correct factor .*)
1.2778 -(*. (expanded polynomial representation) .*)
1.2779 -fun step_add_list_of_fractions_exp [] = (Free("0",HOLogic.realT),[]:term list)
1.2780 - | step_add_list_of_fractions_exp [x] = error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")
1.2781 - | step_add_list_of_fractions_exp (xs)=
1.2782 - let
1.2783 - val den_list=termlist2denominators_exp (xs); (* list of denominators *)
1.2784 - val (denom,var)=calc_lcm(den_list); (* common denominator *)
1.2785 - val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
1.2786 - in
1.2787 - com_den_exp(xs,denom,den,var)
1.2788 - end;
1.2789 -
1.2790 -(* wird aktuell nicht mehr gebraucht, bei rückänderung schon
1.2791 --------------------------------------------------------------
1.2792 -(* WN0210???SK brauch ma des überhaupt *)
1.2793 -fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)
1.2794 - | step_add_list_of_fractions2 [x]=(x,[])
1.2795 - | step_add_list_of_fractions2 (xs) =
1.2796 - let
1.2797 - val den_list=termlist2denominators (xs); (* list of denominators *)
1.2798 - val (denom,var)=calc_lcm(den_list); (* common denominator *)
1.2799 - val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
1.2800 - in
1.2801 - (
1.2802 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2803 - com_den2(xs,denom, poly2term(denom,var)(*den*),var) $
1.2804 - poly2term(denom,var)
1.2805 - ,
1.2806 - []
1.2807 - )
1.2808 - end;
1.2809 -
1.2810 -(* WN0210???SK brauch ma des überhaupt *)
1.2811 -fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)
1.2812 - | step_add_list_of_fractions2_exp [x]=(x,[])
1.2813 - | step_add_list_of_fractions2_exp (xs) =
1.2814 - let
1.2815 - val den_list=termlist2denominators_exp (xs); (* list of denominators *)
1.2816 - val (denom,var)=calc_lcm(den_list); (* common denominator *)
1.2817 - val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
1.2818 - in
1.2819 - (
1.2820 - Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2821 - com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $
1.2822 - poly2expanded(denom,var)
1.2823 - ,
1.2824 - []
1.2825 - )
1.2826 - end;
1.2827 ----------------------------------------------- *)
1.2828 -
1.2829 -
1.2830 -(* converts a term, which contains several terms seperated by +, into a list of these terms .*)
1.2831 -fun term2list (t as (Const("Fields.inverse_class.divide",_) $ _ $ _)) = [t]
1.2832 - | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =
1.2833 - [Const ("Fields.inverse_class.divide",
1.2834 - [HOLogic.realT,HOLogic.realT] ---> HOLogic.realT) $
1.2835 - t $ Free("1",HOLogic.realT)
1.2836 - ]
1.2837 - | term2list (t as (Free(_,_))) =
1.2838 - [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2839 - t $ Free("1",HOLogic.realT)
1.2840 - ]
1.2841 - | term2list (t as (Const("Groups.times_class.times",_) $ _ $ _)) =
1.2842 - [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1.2843 - t $ Free("1",HOLogic.realT)
1.2844 - ]
1.2845 - | term2list (Const("Groups.plus_class.plus",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)
1.2846 - | term2list (Const("Groups.minus_class.minus",_) $ t1 $ t2) =
1.2847 - error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")
1.2848 - | term2list _ = error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");
1.2849 -
1.2850 -(*.factors out the gcd of nominator and denominator:
1.2851 - a/b = (a' * gcd)/(b' * gcd), a,b,gcd are poly[2].*)
1.2852 -
1.2853 -(*. brings the term into a normal form .*)
1.2854 -fun norm_rational_ (thy:theory) t =
1.2855 - SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;
1.2856 -fun norm_expanded_rat_ (thy:theory) t =
1.2857 - SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;
1.2858 -
1.2859 -
1.2860 -(*.evaluates conditions in calculate_Rational.*)
1.2861 -(*make local with FIXX@ME result:term *term list*)
1.2862 -val calc_rat_erls = prep_rls(
1.2863 - Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.2864 - erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.2865 - rules =
1.2866 - [Calc ("HOL.eq",eval_equal "#equal_"),
1.2867 - Calc ("Atools.is'_const",eval_const "#is_const_"),
1.2868 - Thm ("not_true",num_str @{thm not_true}),
1.2869 - Thm ("not_false",num_str @{thm not_false})
1.2870 - ],
1.2871 - scr = EmptyScr});
1.2872 -
1.2873 -
1.2874 -(*.simplifies expressions with numerals;
1.2875 - does NOT rearrange the term by AC-rewriting; thus terms with variables
1.2876 - need to have constants to be commuted together respectively.*)
1.2877 -val calculate_Rational = prep_rls (merge_rls "calculate_Rational"
1.2878 - (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.2879 - erls = calc_rat_erls, srls = Erls,
1.2880 - calc = [], errpatts = [],
1.2881 - rules =
1.2882 - [Calc ("Fields.inverse_class.divide",eval_cancel "#divide_e"),
1.2883 -
1.2884 - Thm ("minus_divide_left",num_str (@{thm minus_divide_left} RS @{thm sym})),
1.2885 - (*SYM - ?x / ?y = - (?x / ?y) may come from subst*)
1.2886 -
1.2887 - Thm ("rat_add",num_str @{thm rat_add}),
1.2888 - (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \
1.2889 - \a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)
1.2890 - Thm ("rat_add1",num_str @{thm rat_add1}),
1.2891 - (*"[| a is_const; b is_const; c is_const |] ==> a / c + b / c = (a + b) / c"*)
1.2892 - Thm ("rat_add2",num_str @{thm rat_add2}),
1.2893 - (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> ?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)
1.2894 - Thm ("rat_add3",num_str @{thm rat_add3}),
1.2895 - (*"[| a is_const; b is_const; c is_const |] ==> a + b / c = (a * c) / c + b / c"\
1.2896 - .... is_const to be omitted here FIXME*)
1.2897 -
1.2898 - Thm ("rat_mult",num_str @{thm rat_mult}),
1.2899 - (*a / b * (c / d) = a * c / (b * d)*)
1.2900 - Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),
1.2901 - (*?x * (?y / ?z) = ?x * ?y / ?z*)
1.2902 - Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),
1.2903 - (*?y / ?z * ?x = ?y * ?x / ?z*)
1.2904 -
1.2905 - Thm ("real_divide_divide1",num_str @{thm real_divide_divide1}),
1.2906 - (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)
1.2907 - Thm ("divide_divide_eq_left",num_str @{thm divide_divide_eq_left}),
1.2908 - (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
1.2909 -
1.2910 - Thm ("rat_power", num_str @{thm rat_power}),
1.2911 - (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
1.2912 -
1.2913 - Thm ("mult_cross",num_str @{thm mult_cross}),
1.2914 - (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)
1.2915 - Thm ("mult_cross1",num_str @{thm mult_cross1}),
1.2916 - (*" b ~= 0 ==> (a / b = c ) = (a = b * c)*)
1.2917 - Thm ("mult_cross2",num_str @{thm mult_cross2})
1.2918 - (*" d ~= 0 ==> (a = c / d) = (a * d = c)*)
1.2919 - ], scr = EmptyScr})
1.2920 - calculate_Poly);
1.2921 + Thm ("minus_divide_left", num_str (@{thm minus_divide_left} RS @{thm sym})),
1.2922 + (*SYM - ?x / ?y = - (?x / ?y) may come from subst*)
1.2923 + Thm ("rat_add", num_str @{thm rat_add}),
1.2924 + (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \
1.2925 + \a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)
1.2926 + Thm ("rat_add1", num_str @{thm rat_add1}),
1.2927 + (*"[| a is_const; b is_const; c is_const |] ==> a / c + b / c = (a + b) / c"*)
1.2928 + Thm ("rat_add2", num_str @{thm rat_add2}),
1.2929 + (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> ?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)
1.2930 + Thm ("rat_add3", num_str @{thm rat_add3}),
1.2931 + (*"[| a is_const; b is_const; c is_const |] ==> a + b / c = (a * c) / c + b / c"\
1.2932 + .... is_const to be omitted here FIXME*)
1.2933 +
1.2934 + Thm ("rat_mult", num_str @{thm rat_mult}),
1.2935 + (*a / b * (c / d) = a * c / (b * d)*)
1.2936 + Thm ("times_divide_eq_right", num_str @{thm times_divide_eq_right}),
1.2937 + (*?x * (?y / ?z) = ?x * ?y / ?z*)
1.2938 + Thm ("times_divide_eq_left", num_str @{thm times_divide_eq_left}),
1.2939 + (*?y / ?z * ?x = ?y * ?x / ?z*)
1.2940 +
1.2941 + Thm ("real_divide_divide1", num_str @{thm real_divide_divide1}),
1.2942 + (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)
1.2943 + Thm ("divide_divide_eq_left", num_str @{thm divide_divide_eq_left}),
1.2944 + (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
1.2945 +
1.2946 + Thm ("rat_power", num_str @{thm rat_power}),
1.2947 + (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
1.2948 +
1.2949 + Thm ("mult_cross", num_str @{thm mult_cross}),
1.2950 + (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)
1.2951 + Thm ("mult_cross1", num_str @{thm mult_cross1}),
1.2952 + (*" b ~= 0 ==> (a / b = c ) = (a = b * c)*)
1.2953 + Thm ("mult_cross2", num_str @{thm mult_cross2})
1.2954 + (*" d ~= 0 ==> (a = c / d) = (a * d = c)*)],
1.2955 + scr = EmptyScr})
1.2956 + calculate_Poly);
1.2957
1.2958 (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)
1.2959 fun eval_is_expanded (thmid:string) _
1.2960 @@ -2980,309 +461,116 @@
1.2961 Trueprop $ (mk_equality (t, @{term False})))
1.2962 | eval_is_expanded _ _ _ _ = NONE;
1.2963
1.2964 +calclist':= overwritel (!calclist',
1.2965 + [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))]);
1.2966 +
1.2967 val rational_erls =
1.2968 - merge_rls "rational_erls" calculate_Rational
1.2969 - (append_rls "is_expanded" Atools_erls
1.2970 - [Calc ("Rational.is'_expanded", eval_is_expanded "")
1.2971 - ]);
1.2972 + merge_rls "rational_erls" calculate_Rational
1.2973 + (append_rls "is_expanded" Atools_erls
1.2974 + [Calc ("Rational.is'_expanded", eval_is_expanded "")]);
1.2975 +*}
1.2976
1.2977 +subsection {* Embed cancellation into rewriting *}
1.2978 +ML {*
1.2979 +local (* cancel_p *)
1.2980
1.2981 -(*.3 'reverse-rewrite-sets' for symbolic computation on rationals:
1.2982 - =================================================================
1.2983 - A[2] 'cancel_p': .
1.2984 - A[3] 'cancel': .
1.2985 - B[2] 'common_nominator_p': transforms summands in a term [2]
1.2986 - to fractions with the (least) common multiple as nominator.
1.2987 - B[3] 'norm_rational': normalizes arbitrary algebraic terms (without
1.2988 - radicals and transzendental functions) to one canceled fraction,
1.2989 - nominator and denominator in polynomial form.
1.2990 +val {rules = rules, rew_ord = (_, ro), ...} = rep_rls (assoc_rls "rev_rew_p");
1.2991
1.2992 -In order to meet isac's requirements for interactive and stepwise calculation,
1.2993 -each 'reverse-rewerite-set' consists of an initialization for the interpreter
1.2994 -state and of 4 functions, each of which employs rewriting as much as possible.
1.2995 -The signature of these functions are the same in each 'reverse-rewrite-set'
1.2996 -respectively.*)
1.2997 +fun init_state thy eval_rls ro t =
1.2998 + let
1.2999 + val SOME (t', _) = factout_p_ thy t;
1.3000 + val SOME (t'', asm) = cancel_p_ thy t;
1.3001 + val der = reverse_deriv thy eval_rls rules ro NONE t';
1.3002 + val der = der @
1.3003 + [(Thm ("real_mult_div_cancel2", num_str @{thm real_mult_div_cancel2}), (t'', asm))]
1.3004 + val rs = (distinct_Thm o (map #1)) der
1.3005 + val rs = filter_out (eq_Thms
1.3006 + ["sym_real_add_zero_left", "sym_real_mult_0", "sym_real_mult_1"]) rs
1.3007 + in (t, t'', [rs(*one in order to ease locate_rule*)], der) end;
1.3008
1.3009 -(* ************************************************************************* *)
1.3010 +fun locate_rule thy eval_rls ro [rs] t r =
1.3011 + if member op = ((map (id_of_thm)) rs) (id_of_thm r)
1.3012 + then
1.3013 + let val ropt = rewrite_ thy ro eval_rls true (thm_of_thm r) t;
1.3014 + in
1.3015 + case ropt of SOME ta => [(r, ta)]
1.3016 + | NONE => (tracing
1.3017 + ("### locate_rule: rewrite " ^ id_of_thm r ^ " " ^ term2str t ^ " = NONE"); [])
1.3018 + end
1.3019 + else (tracing ("### locate_rule: " ^ id_of_thm r ^ " not mem rrls"); [])
1.3020 + | locate_rule _ _ _ _ _ _ = error "locate_rule: doesnt match rev-sets in istate";
1.3021
1.3022 -local(*. cancel_p
1.3023 -------------------------
1.3024 -cancels a single fraction consisting of two (uni- or multivariate)
1.3025 -polynomials WN0609???SK[2] into another such a fraction; examples:
1.3026 +fun next_rule thy eval_rls ro [rs] t =
1.3027 + let
1.3028 + val der = make_deriv thy eval_rls rs ro NONE t;
1.3029 + in case der of (_, r, _) :: _ => SOME r | _ => NONE end
1.3030 + | next_rule _ _ _ _ _ = error ("next_rule: doesnt match rev-sets in istate");
1.3031
1.3032 - a^2 + -1*b^2 a + b
1.3033 - -------------------- = ---------
1.3034 - a^2 + -2*a*b + b^2 a + -1*b
1.3035 -
1.3036 - a^2 a
1.3037 - --- = ---
1.3038 - a 1
1.3039 -
1.3040 -Remark: the reverse ruleset does _NOT_ work properly with other input !.*)
1.3041 -(*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)
1.3042 -
1.3043 -val {rules, rew_ord=(_,ro),...} =
1.3044 - rep_rls (assoc_rls "make_polynomial");
1.3045 -(*WN060829 ... make_deriv does not terminate with 1st expl above,
1.3046 - see rational.sml --- investigate rulesets for cancel_p ---*)
1.3047 -val {rules, rew_ord=(_,ro),...} =
1.3048 - rep_rls (assoc_rls "rev_rew_p");
1.3049 -
1.3050 -(*.init_state = fn : term -> istate
1.3051 -initialzies the state of the script interpreter. The state is:
1.3052 -
1.3053 -type rrlsstate = (*state for reverse rewriting*)
1.3054 - (term * (*the current formula*)
1.3055 - term * (*the final term*)
1.3056 - rule list (*'reverse rule list' (#)*)
1.3057 - list * (*may be serveral, eg. in norm_rational*)
1.3058 - (rule * (*Thm (+ Thm generated from Calc) resulting in ...*)
1.3059 - (term * (*... rewrite with ...*)
1.3060 - term list)) (*... assumptions*)
1.3061 - list); (*derivation from given term to normalform
1.3062 - in reverse order with sym_thm;
1.3063 - (#) could be extracted from here by (map #1)*).*)
1.3064 -(* val {rules, rew_ord=(_,ro),...} =
1.3065 - rep_rls (assoc_rls "rev_rew_p") (*USE ALWAYS, SEE val cancel_p*);
1.3066 - val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*);
1.3067 - val t = t;
1.3068 - *)
1.3069 -fun init_state thy eval_rls ro t =
1.3070 - let val SOME (t',_) = factout_p_ thy t
1.3071 - val SOME (t'',asm) = cancel_p_ thy t
1.3072 - val der = reverse_deriv thy eval_rls rules ro NONE t'
1.3073 - val der = der @ [(Thm ("real_mult_div_cancel2",
1.3074 - num_str @{thm real_mult_div_cancel2}),
1.3075 - (t'',asm))]
1.3076 - val rs = (distinct_Thm o (map #1)) der
1.3077 - val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
1.3078 - "sym_real_mult_0",
1.3079 - "sym_real_mult_1"
1.3080 - (*..insufficient,eg.make_Polynomial*)])rs
1.3081 - in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end;
1.3082 -
1.3083 -(*.locate_rule = fn : rule list -> term -> rule
1.3084 - -> (rule * (term * term list) option) list.
1.3085 - checks a rule R for being a cancel-rule, and if it is,
1.3086 - then return the list of rules (+ the terms they are rewriting to)
1.3087 - which need to be applied before R should be applied.
1.3088 - precondition: the rule is applicable to the argument-term.
1.3089 -arguments:
1.3090 - rule list: the reverse rule list
1.3091 - -> term : ... to which the rule shall be applied
1.3092 - -> rule : ... to be applied to term
1.3093 -value:
1.3094 - -> (rule : a rule rewriting to ...
1.3095 - * (term : ... the resulting term ...
1.3096 - * term list): ... with the assumptions ( //#0).
1.3097 - ) list : there may be several such rules;
1.3098 - the list is empty, if the rule has nothing to do
1.3099 - with cancelation.*)
1.3100 -(* val () = ();
1.3101 - *)
1.3102 -fun locate_rule thy eval_rls ro [rs] t r =
1.3103 - if (id_of_thm r) mem (map (id_of_thm)) rs
1.3104 - then let val ropt =
1.3105 - rewrite_ thy ro eval_rls true (thm_of_thm r) t;
1.3106 - in case ropt of
1.3107 - SOME ta => [(r, ta)]
1.3108 - | NONE => (tracing("### locate_rule: rewrite "^
1.3109 - (id_of_thm r)^" "^(term2str t)^" = NONE");
1.3110 - []) end
1.3111 - else (tracing("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
1.3112 - | locate_rule _ _ _ _ _ _ =
1.3113 - error ("locate_rule: doesnt match rev-sets in istate");
1.3114 -
1.3115 -(*.next_rule = fn : rule list -> term -> rule option
1.3116 - for a given term return the next rules to be done for cancelling.
1.3117 -arguments:
1.3118 - rule list : the reverse rule list
1.3119 - term : the term for which ...
1.3120 -value:
1.3121 - -> rule option: ... this rule is appropriate for cancellation;
1.3122 - there may be no such rule (if the term is canceled already.*)
1.3123 -(* val thy = thy;
1.3124 - val Rrls {rew_ord=(_,ro),...} = cancel;
1.3125 - val ([rs],t) = (rss,f);
1.3126 - next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
1.3127 -
1.3128 - val (thy, [rs]) = (thy, revsets);
1.3129 - val Rrls {rew_ord=(_,ro),...} = cancel;
1.3130 - nex [rs] t;
1.3131 - *)
1.3132 -fun next_rule thy eval_rls ro [rs] t =
1.3133 - let val der = make_deriv thy eval_rls rs ro NONE t;
1.3134 - in case der of
1.3135 -(* val (_,r,_)::_ = der;
1.3136 - *)
1.3137 - (_,r,_)::_ => SOME r
1.3138 - | _ => NONE
1.3139 - end
1.3140 - | next_rule _ _ _ _ _ =
1.3141 - error ("next_rule: doesnt match rev-sets in istate");
1.3142 -
1.3143 -(*.val attach_form = f : rule list -> term -> term
1.3144 - -> (rule * (term * term list)) list
1.3145 - checks an input term TI, if it may belong to a current cancellation, by
1.3146 - trying to derive it from the given term TG.
1.3147 -arguments:
1.3148 - term : TG, the last one in the cancellation agreed upon by user + math-eng
1.3149 - -> term: TI, the next one input by the user
1.3150 -value:
1.3151 - -> (rule : the rule to be applied in order to reach TI
1.3152 - * (term : ... obtained by applying the rule ...
1.3153 - * term list): ... and the respective assumptions.
1.3154 - ) list : there may be several such rules;
1.3155 - the list is empty, if the users term does not belong
1.3156 - to a cancellation of the term last agreed upon.*)
1.3157 -fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
1.3158 - []:(rule * (term * term list)) list;
1.3159 +fun attach_form (_:rule list list) (_:term) (_:term) =
1.3160 + [(*TODO*)]: (rule * (term * term list)) list;
1.3161
1.3162 in
1.3163
1.3164 -val cancel_p =
1.3165 - Rrls {id = "cancel_p", prepat=[],
1.3166 - rew_ord=("ord_make_polynomial",
1.3167 - ord_make_polynomial false thy),
1.3168 - erls = rational_erls,
1.3169 - calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")),
1.3170 - ("TIMES" ,("Groups.times_class.times" ,eval_binop "#mult_")),
1.3171 - ("DIVIDE" ,("Fields.inverse_class.divide" ,eval_cancel "#divide_e")),
1.3172 - ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
1.3173 - errpatts = [],
1.3174 - scr=Rfuns {init_state = init_state thy Atools_erls ro,
1.3175 - normal_form = cancel_p_ thy,
1.3176 - locate_rule = locate_rule thy Atools_erls ro,
1.3177 - next_rule = next_rule thy Atools_erls ro,
1.3178 - attach_form = attach_form}}
1.3179 -end;(*local*)
1.3180 +val cancel_p =
1.3181 + Rrls {id = "cancel_p", prepat = [],
1.3182 + rew_ord=("ord_make_polynomial", ord_make_polynomial false thy),
1.3183 + erls = rational_erls,
1.3184 + calc =
1.3185 + [("PLUS", ("Groups.plus_class.plus", eval_binop "#add_")),
1.3186 + ("TIMES" , ("Groups.times_class.times", eval_binop "#mult_")),
1.3187 + ("DIVIDE", ("Fields.inverse_class.divide", eval_cancel "#divide_e")),
1.3188 + ("POWER", ("Atools.pow", eval_binop "#power_"))],
1.3189 + errpatts = [],
1.3190 + scr =
1.3191 + Rfuns {init_state = init_state thy Atools_erls ro,
1.3192 + normal_form = cancel_p_ thy,
1.3193 + locate_rule = locate_rule thy Atools_erls ro,
1.3194 + next_rule = next_rule thy Atools_erls ro,
1.3195 + attach_form = attach_form}}
1.3196 +end; (* local cancel_p *)
1.3197 +*}
1.3198
1.3199 -local(*.ad [2] 'common_nominator_p'
1.3200 ----------------------------------
1.3201 -FIXME Beschreibung .*)
1.3202 +subsection {* Embed addition into rewriting *}
1.3203 +ML {*
1.3204 +local (* add_fractions_p *)
1.3205
1.3206 +val {rules = rules, rew_ord = (_, ro), ...} = rep_rls (assoc_rls "make_polynomial");
1.3207 +val {rules, rew_ord=(_,ro),...} = rep_rls (assoc_rls "rev_rew_p");
1.3208
1.3209 -val {rules=rules,rew_ord=(_,ro),...} =
1.3210 - rep_rls (assoc_rls "make_polynomial");
1.3211 -(*WN060829 ... make_deriv does not terminate with 1st expl above,
1.3212 - see rational.sml --- investigate rulesets for cancel_p ---*)
1.3213 -val {rules, rew_ord=(_,ro),...} =
1.3214 - rep_rls (assoc_rls "rev_rew_p");
1.3215 -val thy = thy;
1.3216 +fun init_state thy eval_rls ro t =
1.3217 + let
1.3218 + val SOME (t',_) = common_nominator_p_ thy t;
1.3219 + val SOME (t'', asm) = add_fraction_p_ thy t;
1.3220 + val der = reverse_deriv thy eval_rls rules ro NONE t';
1.3221 + val der = der @
1.3222 + [(Thm ("real_mult_div_cancel2", num_str @{thm real_mult_div_cancel2}), (t'',asm))]
1.3223 + val rs = (distinct_Thm o (map #1)) der;
1.3224 + val rs = filter_out (eq_Thms
1.3225 + ["sym_real_add_zero_left", "sym_real_mult_0", "sym_real_mult_1"]) rs;
1.3226 + in (t, t'', [rs(*here only _ONE_*)], der) end;
1.3227
1.3228 +fun locate_rule thy eval_rls ro [rs] t r =
1.3229 + if member op = ((map (id_of_thm)) rs) (id_of_thm r)
1.3230 + then
1.3231 + let val ropt = rewrite_ thy ro eval_rls true (thm_of_thm r) t;
1.3232 + in
1.3233 + case ropt of
1.3234 + SOME ta => [(r, ta)]
1.3235 + | NONE =>
1.3236 + (tracing ("### locate_rule: rewrite " ^ id_of_thm r ^ " " ^ term2str t ^ " = NONE");
1.3237 + []) end
1.3238 + else (tracing ("### locate_rule: " ^ id_of_thm r ^ " not mem rrls"); [])
1.3239 + | locate_rule _ _ _ _ _ _ = error "locate_rule: doesnt match rev-sets in istate";
1.3240
1.3241 -(*.common_nominator_p_ = fn : theory -> term -> (term * term list) option
1.3242 - as defined above*)
1.3243 -
1.3244 -(*.init_state = fn : term -> istate
1.3245 -initialzies the state of the interactive interpreter. The state is:
1.3246 -
1.3247 -type rrlsstate = (*state for reverse rewriting*)
1.3248 - (term * (*the current formula*)
1.3249 - term * (*the final term*)
1.3250 - rule list (*'reverse rule list' (#)*)
1.3251 - list * (*may be serveral, eg. in norm_rational*)
1.3252 - (rule * (*Thm (+ Thm generated from Calc) resulting in ...*)
1.3253 - (term * (*... rewrite with ...*)
1.3254 - term list)) (*... assumptions*)
1.3255 - list); (*derivation from given term to normalform
1.3256 - in reverse order with sym_thm;
1.3257 - (#) could be extracted from here by (map #1)*).*)
1.3258 -fun init_state thy eval_rls ro t =
1.3259 - let val SOME (t',_) = common_nominator_p_ thy t;
1.3260 - val SOME (t'',asm) = add_fraction_p_ thy t;
1.3261 - val der = reverse_deriv thy eval_rls rules ro NONE t';
1.3262 - val der = der @ [(Thm ("real_mult_div_cancel2",
1.3263 - num_str @{thm real_mult_div_cancel2}),
1.3264 - (t'',asm))]
1.3265 - val rs = (distinct_Thm o (map #1)) der;
1.3266 - val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
1.3267 - "sym_real_mult_0",
1.3268 - "sym_real_mult_1"]) rs;
1.3269 - in (t,t'',[rs(*here only _ONE_*)],der) end;
1.3270 -
1.3271 -(* use"knowledge/Rational.ML";
1.3272 - *)
1.3273 -
1.3274 -(*.locate_rule = fn : rule list -> term -> rule
1.3275 - -> (rule * (term * term list) option) list.
1.3276 - checks a rule R for being a cancel-rule, and if it is,
1.3277 - then return the list of rules (+ the terms they are rewriting to)
1.3278 - which need to be applied before R should be applied.
1.3279 - precondition: the rule is applicable to the argument-term.
1.3280 -arguments:
1.3281 - rule list: the reverse rule list
1.3282 - -> term : ... to which the rule shall be applied
1.3283 - -> rule : ... to be applied to term
1.3284 -value:
1.3285 - -> (rule : a rule rewriting to ...
1.3286 - * (term : ... the resulting term ...
1.3287 - * term list): ... with the assumptions ( //#0).
1.3288 - ) list : there may be several such rules;
1.3289 - the list is empty, if the rule has nothing to do
1.3290 - with cancelation.*)
1.3291 -(* val () = ();
1.3292 - *)
1.3293 -fun locate_rule thy eval_rls ro [rs] t r =
1.3294 - if (id_of_thm r) mem (map (id_of_thm)) rs
1.3295 - then let val ropt =
1.3296 - rewrite_ thy ro eval_rls true (thm_of_thm r) t;
1.3297 - in case ropt of
1.3298 - SOME ta => [(r, ta)]
1.3299 - | NONE => (tracing("### locate_rule: rewrite "^
1.3300 - (id_of_thm r)^" "^(term2str t)^" = NONE");
1.3301 - []) end
1.3302 - else (tracing("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
1.3303 - | locate_rule _ _ _ _ _ _ =
1.3304 - error ("locate_rule: doesnt match rev-sets in istate");
1.3305 -
1.3306 -(*.next_rule = fn : rule list -> term -> rule option
1.3307 - for a given term return the next rules to be done for cancelling.
1.3308 -arguments:
1.3309 - rule list : the reverse rule list
1.3310 - term : the term for which ...
1.3311 -value:
1.3312 - -> rule option: ... this rule is appropriate for cancellation;
1.3313 - there may be no such rule (if the term is canceled already.*)
1.3314 -(* val thy = thy;
1.3315 - val Rrls {rew_ord=(_,ro),...} = cancel;
1.3316 - val ([rs],t) = (rss,f);
1.3317 - next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
1.3318 -
1.3319 - val (thy, [rs]) = (thy, revsets);
1.3320 - val Rrls {rew_ord=(_,ro),...} = cancel;
1.3321 - nex [rs] t;
1.3322 - *)
1.3323 fun next_rule thy eval_rls ro [rs] t =
1.3324 let val der = make_deriv thy eval_rls rs ro NONE t;
1.3325 - in case der of
1.3326 -(* val (_,r,_)::_ = der;
1.3327 - *)
1.3328 - (_,r,_)::_ => SOME r
1.3329 - | _ => NONE
1.3330 + in
1.3331 + case der of
1.3332 + (_,r,_)::_ => SOME r
1.3333 + | _ => NONE
1.3334 end
1.3335 - | next_rule _ _ _ _ _ =
1.3336 - error ("next_rule: doesnt match rev-sets in istate");
1.3337 + | next_rule _ _ _ _ _ = error ("next_rule: doesnt match rev-sets in istate");
1.3338
1.3339 -(*.val attach_form = f : rule list -> term -> term
1.3340 - -> (rule * (term * term list)) list
1.3341 - checks an input term TI, if it may belong to a current cancellation, by
1.3342 - trying to derive it from the given term TG.
1.3343 -arguments:
1.3344 - term : TG, the last one in the cancellation agreed upon by user + math-eng
1.3345 - -> term: TI, the next one input by the user
1.3346 -value:
1.3347 - -> (rule : the rule to be applied in order to reach TI
1.3348 - * (term : ... obtained by applying the rule ...
1.3349 - * term list): ... and the respective assumptions.
1.3350 - ) list : there may be several such rules;
1.3351 - the list is empty, if the users term does not belong
1.3352 - to a cancellation of the term last agreed upon.*)
1.3353 -fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
1.3354 - []:(rule * (term * term list)) list;
1.3355 -
1.3356 -(* if each pat* matches with the current term, the Rrls is applied
1.3357 - (there are no preconditions to be checked, they are @{term True}) *)
1.3358 val pat0 = parse_patt thy "?r/?s+?u/?v :: real";
1.3359 val pat1 = parse_patt thy "?r/?s+?u :: real";
1.3360 val pat2 = parse_patt thy "?r +?u/?v :: real";
1.3361 @@ -3291,98 +579,34 @@
1.3362 ([@{term True}], pat2)];
1.3363 in
1.3364
1.3365 -(*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;
1.3366 - besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht
1.3367 - dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)
1.3368 -val common_nominator_p =
1.3369 - Rrls {id = "common_nominator_p", prepat=prepat,
1.3370 - rew_ord=("ord_make_polynomial",
1.3371 - ord_make_polynomial false thy),
1.3372 - erls = rational_erls,
1.3373 - calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")),
1.3374 - ("TIMES" ,("Groups.times_class.times" ,eval_binop "#mult_")),
1.3375 - ("DIVIDE" ,("Fields.inverse_class.divide" ,eval_cancel "#divide_e")),
1.3376 - ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
1.3377 - errpatts = [],
1.3378 - scr=Rfuns {init_state = init_state thy Atools_erls ro,
1.3379 - normal_form = add_fraction_p_ thy, (*FIXME.WN0211*)
1.3380 - locate_rule = locate_rule thy Atools_erls ro,
1.3381 - next_rule = next_rule thy Atools_erls ro,
1.3382 - attach_form = attach_form}}
1.3383 -end;(*local*)
1.3384 +val add_fractions_p =
1.3385 + Rrls {id = "add_fractions_p", prepat=prepat,
1.3386 + rew_ord = ("ord_make_polynomial", ord_make_polynomial false thy),
1.3387 + erls = rational_erls,
1.3388 + calc = [("PLUS", ("Groups.plus_class.plus", eval_binop "#add_")),
1.3389 + ("TIMES", ("Groups.times_class.times", eval_binop "#mult_")),
1.3390 + ("DIVIDE", ("Fields.inverse_class.divide", eval_cancel "#divide_e")),
1.3391 + ("POWER", ("Atools.pow", eval_binop "#power_"))],
1.3392 + errpatts = [],
1.3393 + scr = Rfuns {init_state = init_state thy Atools_erls ro,
1.3394 + normal_form = add_fraction_p_ thy,
1.3395 + locate_rule = locate_rule thy Atools_erls ro,
1.3396 + next_rule = next_rule thy Atools_erls ro,
1.3397 + attach_form = attach_form}}
1.3398 +end; (*local add_fractions_p *)
1.3399 +*}
1.3400
1.3401 -end; (*struct*)
1.3402 +subsection {* Cancelling and adding all occurrences in a term /////////////////////////////*}
1.3403 +ML {*
1.3404 +(*copying cancel_p_rls + add her caused error in interface.sml*)
1.3405 +*}
1.3406
1.3407 -*}
1.3408 +section {* Rulesets for general simplification *}
1.3409 ML {*
1.3410 -open RationalI;
1.3411 -(*##*)
1.3412
1.3413 -(*.the expression contains + - * ^ / only ?.*)
1.3414 -fun is_ratpolyexp (Free _) = true
1.3415 - | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ Free _ $ Free _) = true
1.3416 - | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ Free _ $ Free _) = true
1.3417 - | is_ratpolyexp (Const ("Groups.times_class.times",_) $ Free _ $ Free _) = true
1.3418 - | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true
1.3419 - | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ Free _ $ Free _) = true
1.3420 - | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ t1 $ t2) =
1.3421 - ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.3422 - | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) =
1.3423 - ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.3424 - | is_ratpolyexp (Const ("Groups.times_class.times",_) $ t1 $ t2) =
1.3425 - ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.3426 - | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =
1.3427 - ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.3428 - | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ t1 $ t2) =
1.3429 - ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
1.3430 - | is_ratpolyexp _ = false;
1.3431 -
1.3432 -(*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)
1.3433 -fun eval_is_ratpolyexp (thmid:string) _
1.3434 - (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =
1.3435 - if is_ratpolyexp arg
1.3436 - then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",
1.3437 - Trueprop $ (mk_equality (t, @{term True})))
1.3438 - else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",
1.3439 - Trueprop $ (mk_equality (t, @{term False})))
1.3440 - | eval_is_ratpolyexp _ _ _ _ = NONE;
1.3441 -
1.3442 -(*("get_denominator", ("Rational.get_denominator", eval_get_denominator ""))*)
1.3443 -fun eval_get_denominator (thmid:string) _
1.3444 - (t as Const ("Rational.get_denominator", _) $
1.3445 - (Const ("Fields.inverse_class.divide", _) $ num $
1.3446 - denom)) thy =
1.3447 - SOME (mk_thmid thmid "" (term_to_string''' thy denom) "",
1.3448 - Trueprop $ (mk_equality (t, denom)))
1.3449 - | eval_get_denominator _ _ _ _ = NONE;
1.3450 -
1.3451 -(*("get_numerator", ("Rational.get_numerator", eval_get_numerator ""))*)
1.3452 -fun eval_get_numerator (thmid:string) _
1.3453 - (t as Const ("Rational.get_numerator", _) $
1.3454 - (Const ("Fields.inverse_class.divide", _) $num
1.3455 - $denom )) thy =
1.3456 - SOME (mk_thmid thmid "" (term_to_string''' thy num) "",
1.3457 - Trueprop $ (mk_equality (t, num)))
1.3458 - | eval_get_numerator _ _ _ _ = NONE;
1.3459 -
1.3460 -(*-------------------18.3.03 --> struct <-----------vvv--*)
1.3461 +(*-------------------18.3.03 --> struct <-----------vvv--
1.3462 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)
1.3463 -
1.3464 -(*WN100906 removed in favour of discard_minus = discard_minus_ formerly
1.3465 - discard binary minus, shift unary minus into -1*;
1.3466 - unary minus before numerals are put into the numeral by parsing;
1.3467 - contains absolute minimum of thms for context in norm_Rational
1.3468 -*** val discard_minus = prep_rls(
1.3469 - Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.3470 - erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3471 - rules =
1.3472 - [Thm ("real_diff_minus", num_str @{thm real_diff_minus}),
1.3473 - (*"a - b = a + -1 * b"*)
1.3474 - Thm ("sym_real_mult_minus1", num_str (@{thm real_mult_minus1} RS @{thm sym}))
1.3475 - (*- ?z = "-1 * ?z"*)],
1.3476 - scr = EmptyScr
1.3477 - }):rls;
1.3478 - *)
1.3479 + -------------------18.3.03 --> struct <-----------vvv--*)
1.3480
1.3481 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)
1.3482 val powers_erls = prep_rls(
1.3483 @@ -3433,7 +657,7 @@
1.3484 (*.contains absolute minimum of thms for context in norm_Rational.*)
1.3485 val rat_mult_divide = prep_rls(
1.3486 Rls {id = "rat_mult_divide", preconds = [],
1.3487 - rew_ord = ("dummy_ord",dummy_ord),
1.3488 + rew_ord = ("dummy_ord", dummy_ord),
1.3489 erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3490 rules = [Thm ("rat_mult",num_str @{thm rat_mult}),
1.3491 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
1.3492 @@ -3493,12 +717,11 @@
1.3493 rules = [Calc ("Atools.is'_const",eval_const "#is_const_")
1.3494 ], scr = EmptyScr}:rls);
1.3495
1.3496 -(*.consists of rls containing the absolute minimum of thms.*)
1.3497 +(* consists of rls containing the absolute minimum of thms *)
1.3498 (*040209: this version has been used by RL for his equations,
1.3499 -which is now replaced by MGs version below
1.3500 -vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)
1.3501 -val norm_Rational = prep_rls(
1.3502 - Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.3503 +which is now replaced by MGs version "norm_Rational" below *)
1.3504 +val norm_Rational_min = prep_rls(
1.3505 + Rls {id = "norm_Rational_min", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.3506 erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],
1.3507 rules = [(*sequence given by operator precedence*)
1.3508 Rls_ discard_minus,
1.3509 @@ -3506,7 +729,6 @@
1.3510 Rls_ rat_mult_divide,
1.3511 Rls_ expand,
1.3512 Rls_ reduce_0_1_2,
1.3513 - (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)
1.3514 Rls_ order_add_mult,
1.3515 Rls_ collect_numerals,
1.3516 Rls_ add_fractions_p,
1.3517 @@ -3519,7 +741,7 @@
1.3518 rew_ord = ("dummy_ord", dummy_ord),
1.3519 erls = Atools_erls, srls = Erls,
1.3520 calc = [], errpatts = [],
1.3521 - rules = [Rls_ norm_Rational, (*from RL -- not the latest one*)
1.3522 + rules = [Rls_ norm_Rational_min,
1.3523 Rls_ discard_parentheses
1.3524 ],
1.3525 scr = EmptyScr}:rls);
1.3526 @@ -3543,149 +765,108 @@
1.3527 (*"?a - ?b + ?b = ?a"*)
1.3528 Thm ("divide_minus1",num_str @{thm divide_minus1})
1.3529 (*"?x / -1 = - ?x"*)
1.3530 -(*
1.3531 -,
1.3532 - Thm ("",num_str @{thm })
1.3533 -*)
1.3534 ]);
1.3535 *}
1.3536 ML {*
1.3537 +val add_fractions_p_rls = prep_rls(
1.3538 + Rls {id = "add_fractions_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.3539 + erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3540 + rules = [Rls_ add_fractions_p],
1.3541 + scr = EmptyScr});
1.3542
1.3543 -(*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)
1.3544 -
1.3545 -(* ------------------------------------------------------------------ *)
1.3546 -(* Simplifier für beliebige Buchterme *)
1.3547 -(* ------------------------------------------------------------------ *)
1.3548 -(*----------------------- norm_Rational_mg ---------------------------*)
1.3549 -(*. description of the simplifier see MG-DA.p.56ff .*)
1.3550 -(* ------------------------------------------------------------------- *)
1.3551 -val common_nominator_p_rls = prep_rls(
1.3552 - Rls {id = "common_nominator_p_rls", preconds = [],
1.3553 - rew_ord = ("dummy_ord",dummy_ord),
1.3554 - erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3555 - rules =
1.3556 - [Rls_ common_nominator_p
1.3557 - (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?
1.3558 - FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)
1.3559 - ],
1.3560 - scr = EmptyScr});
1.3561 -(* ------------------------------------------------------------------- *)
1.3562 (* "Rls" causes repeated application of cancel_p to one and the same term *)
1.3563 val cancel_p_rls = prep_rls(
1.3564 Rls
1.3565 {id = "cancel_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.3566 erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3567 - rules =
1.3568 - [Rls_ cancel_p (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)
1.3569 - ],
1.3570 + rules = [Rls_ cancel_p],
1.3571 scr = EmptyScr});
1.3572 -(* -------------------------------------------------------------------- *)
1.3573 +
1.3574 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
1.3575 used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)
1.3576 val rat_mult_poly = prep_rls(
1.3577 - Rls {id = "rat_mult_poly", preconds = [],
1.3578 - rew_ord = ("dummy_ord",dummy_ord),
1.3579 - erls = append_rls "e_rls-is_polyexp" e_rls
1.3580 - [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
1.3581 - srls = Erls, calc = [], errpatts = [],
1.3582 - rules =
1.3583 - [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),
1.3584 - (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
1.3585 - Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})
1.3586 - (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
1.3587 - ],
1.3588 - scr = EmptyScr});
1.3589 + Rls {id = "rat_mult_poly", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.3590 + erls = append_rls "e_rls-is_polyexp" e_rls [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
1.3591 + srls = Erls, calc = [], errpatts = [],
1.3592 + rules =
1.3593 + [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),
1.3594 + (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
1.3595 + Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})
1.3596 + (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) ],
1.3597 + scr = EmptyScr});
1.3598
1.3599 -(* ------------------------------------------------------------------ *)
1.3600 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
1.3601 used in looping part norm_Rational_rls, see example DA-M02-main.p.60
1.3602 .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls,
1.3603 I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028
1.3604 ... WN0609???MG.*)
1.3605 val rat_mult_div_pow = prep_rls(
1.3606 - Rls {id = "rat_mult_div_pow", preconds = [],
1.3607 - rew_ord = ("dummy_ord",dummy_ord),
1.3608 - erls = e_rls,
1.3609 - (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls
1.3610 - [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
1.3611 - with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get
1.3612 - error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.
1.3613 - thus we decided to go on with this flaw*)
1.3614 - srls = Erls, calc = [], errpatts = [],
1.3615 - rules = [Thm ("rat_mult",num_str @{thm rat_mult}),
1.3616 - (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
1.3617 - Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),
1.3618 - (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
1.3619 - Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}),
1.3620 - (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
1.3621 + Rls {id = "rat_mult_div_pow", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.3622 + erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3623 + rules = [Thm ("rat_mult", num_str @{thm rat_mult}),
1.3624 + (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
1.3625 + Thm ("rat_mult_poly_l", num_str @{thm rat_mult_poly_l}),
1.3626 + (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
1.3627 + Thm ("rat_mult_poly_r", num_str @{thm rat_mult_poly_r}),
1.3628 + (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
1.3629 +
1.3630 + Thm ("real_divide_divide1_mg", num_str @{thm real_divide_divide1_mg}),
1.3631 + (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)
1.3632 + Thm ("divide_divide_eq_right", num_str @{thm divide_divide_eq_right}),
1.3633 + (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
1.3634 + Thm ("divide_divide_eq_left", num_str @{thm divide_divide_eq_left}),
1.3635 + (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
1.3636 + Calc ("Fields.inverse_class.divide", eval_cancel "#divide_e"),
1.3637 +
1.3638 + Thm ("rat_power", num_str @{thm rat_power})
1.3639 + (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
1.3640 + ],
1.3641 + scr = EmptyScr}:rls);
1.3642
1.3643 - Thm ("real_divide_divide1_mg",
1.3644 - num_str @{thm real_divide_divide1_mg}),
1.3645 - (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)
1.3646 - Thm ("divide_divide_eq_right",
1.3647 - num_str @{thm divide_divide_eq_right}),
1.3648 - (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
1.3649 - Thm ("divide_divide_eq_left",
1.3650 - num_str @{thm divide_divide_eq_left}),
1.3651 - (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
1.3652 - Calc ("Fields.inverse_class.divide" ,eval_cancel "#divide_e"),
1.3653 -
1.3654 - Thm ("rat_power", num_str @{thm rat_power})
1.3655 - (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
1.3656 - ],
1.3657 - scr = EmptyScr}:rls);
1.3658 -(* ------------------------------------------------------------------ *)
1.3659 val rat_reduce_1 = prep_rls(
1.3660 - Rls {id = "rat_reduce_1", preconds = [],
1.3661 - rew_ord = ("dummy_ord",dummy_ord),
1.3662 - erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3663 - rules = [Thm ("divide_1",num_str @{thm divide_1}),
1.3664 - (*"?x / 1 = ?x"*)
1.3665 - Thm ("mult_1_left",num_str @{thm mult_1_left})
1.3666 - (*"1 * z = z"*)
1.3667 - ],
1.3668 - scr = EmptyScr}:rls);
1.3669 -(* ------------------------------------------------------------------ *)
1.3670 -(*. looping part of norm_Rational(*_mg*) .*)
1.3671 -val norm_Rational_rls = prep_rls(
1.3672 - Rls {id = "norm_Rational_rls", preconds = [],
1.3673 - rew_ord = ("dummy_ord",dummy_ord),
1.3674 - erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],
1.3675 - rules = [Rls_ common_nominator_p_rls,
1.3676 - Rls_ rat_mult_div_pow,
1.3677 - Rls_ make_rat_poly_with_parentheses,
1.3678 - Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)
1.3679 - Rls_ rat_reduce_1
1.3680 - ],
1.3681 - scr = EmptyScr}:rls);
1.3682 -(* ------------------------------------------------------------------ *)
1.3683 -(*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)
1.3684 - just be renaming:*)
1.3685 -val norm_Rational (*_mg*) = prep_rls(
1.3686 + Rls {id = "rat_reduce_1", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.3687 + erls = e_rls, srls = Erls, calc = [], errpatts = [],
1.3688 + rules =
1.3689 + [Thm ("divide_1", num_str @{thm divide_1}),
1.3690 + (*"?x / 1 = ?x"*)
1.3691 + Thm ("mult_1_left", num_str @{thm mult_1_left})
1.3692 + (*"1 * z = z"*)
1.3693 + ],
1.3694 + scr = EmptyScr}:rls);
1.3695 +
1.3696 +(* looping part of norm_Rational *)
1.3697 +val norm_Rational_rls = prep_rls (
1.3698 + Rls {id = "norm_Rational_rls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
1.3699 + erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],
1.3700 + rules = [Rls_ add_fractions_p_rls,
1.3701 + Rls_ rat_mult_div_pow,
1.3702 + Rls_ make_rat_poly_with_parentheses,
1.3703 + Rls_ cancel_p_rls,
1.3704 + Rls_ rat_reduce_1
1.3705 + ],
1.3706 + scr = EmptyScr}:rls);
1.3707 +
1.3708 +val norm_Rational = prep_rls (
1.3709 Seq
1.3710 - {id = "norm_Rational"(*_mg*), preconds = [],
1.3711 - rew_ord = ("dummy_ord",dummy_ord),
1.3712 + {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
1.3713 erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],
1.3714 rules = [Rls_ discard_minus,
1.3715 - Rls_ rat_mult_poly, (* removes double fractions like a/b/c *)
1.3716 - Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)
1.3717 - Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)
1.3718 - Rls_ norm_Rational_rls, (* the main rls, looping (#) *)
1.3719 - Rls_ discard_parentheses1 (* mult only *)
1.3720 + Rls_ rat_mult_poly, (* removes double fractions like a/b/c *)
1.3721 + Rls_ make_rat_poly_with_parentheses,
1.3722 + Rls_ cancel_p_rls,
1.3723 + Rls_ norm_Rational_rls, (* the main rls, looping (#) *)
1.3724 + Rls_ discard_parentheses1 (* mult only *)
1.3725 ],
1.3726 scr = EmptyScr}:rls);
1.3727 -"-------- rls norm_Rational --------------------------------------------------";
1.3728 -(* ------------------------------------------------------------------ *)
1.3729 -
1.3730 *}
1.3731 -ML {*
1.3732 +ML {*
1.3733 ruleset' := overwritelthy @{theory} (!ruleset',
1.3734 [("calculate_Rational", calculate_Rational),
1.3735 ("calc_rat_erls",calc_rat_erls),
1.3736 ("rational_erls", rational_erls),
1.3737 ("cancel_p", cancel_p),
1.3738 - ("common_nominator_p", common_nominator_p),
1.3739 - ("common_nominator_p_rls", common_nominator_p_rls),
1.3740 + ("add_fractions_p", add_fractions_p),
1.3741 + ("add_fractions_p_rls", add_fractions_p_rls),
1.3742 (*WN120410 ("discard_minus", discard_minus), is defined in Poly.thy*)
1.3743 ("powers_erls", powers_erls),
1.3744 ("powers", powers),
1.3745 @@ -3701,12 +882,10 @@
1.3746 ("cancel_p_rls", cancel_p_rls)
1.3747 ]);
1.3748
1.3749 -calclist':= overwritel (!calclist',
1.3750 - [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))
1.3751 - ]);
1.3752 +*}
1.3753
1.3754 -(** problems **)
1.3755 -
1.3756 +section {* A problem for simplification of rationals *}
1.3757 +ML {*
1.3758 store_pbt
1.3759 (prep_pbt thy "pbl_simp_rat" [] e_pblID
1.3760 (["rational","simplification"],
1.3761 @@ -3717,9 +896,10 @@
1.3762 append_rls "e_rls" e_rls [(*for preds in where_*)],
1.3763 SOME "Simplify t_t",
1.3764 [["simplification","of_rationals"]]));
1.3765 +*}
1.3766
1.3767 -(** methods **)
1.3768 -
1.3769 +section {* A methods for simplification of rationals *}
1.3770 +ML {*
1.3771 (*WN061025 this methods script is copied from (auto-generated) script
1.3772 of norm_Rational in order to ease repair on inform*)
1.3773 store_met (prep_met thy "met_simp_rat" [] e_metID (["simplification","of_rationals"],
1.3774 @@ -3737,14 +917,13 @@
1.3775 " Try (Rewrite_Set make_rat_poly_with_parentheses False) @@ " ^
1.3776 " Try (Rewrite_Set cancel_p_rls False) @@ " ^
1.3777 " (Repeat " ^
1.3778 - " ((Try (Rewrite_Set common_nominator_p_rls False) @@ " ^
1.3779 + " ((Try (Rewrite_Set add_fractions_p_rls False) @@ " ^
1.3780 " Try (Rewrite_Set rat_mult_div_pow False) @@ " ^
1.3781 " Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^
1.3782 " Try (Rewrite_Set cancel_p_rls False) @@ " ^
1.3783 " Try (Rewrite_Set rat_reduce_1 False)))) @@ " ^
1.3784 " Try (Rewrite_Set discard_parentheses1 False)) " ^
1.3785 " t_t)"));
1.3786 +*}
1.3787
1.3788 -*}
1.3789 -ML {*"test"*}
1.3790 -end (* theory*)
1.3791 +end