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