1 (*.calculate in rationals: gcd, lcm, etc.
3 Institute for Mathematics D and Institute for Software Technology,
5 Use is subject to license terms.
7 use"Knowledge/Rational.ML";
11 use_thy"Knowledge/Isac";
12 ****************************************************************.*)
14 (*.*****************************************************************
15 Remark on notions in the documentation below:
16 referring to the remark on 'polynomials' in Poly.sml we use
17 [2] 'polynomial' normalform (Polynom)
18 [3] 'expanded_term' normalform (Ausmultiplizierter Term),
19 where normalform [2] is a special case of [3], i.e. [3] implies [2].
21 'fraction with numerator and nominator both in normalform [2]'
22 'fraction with numerator and nominator both in normalform [3]'
24 'fraction in normalform [2]'
25 'fraction in normalform [3]'
29 a 'simple fraction' is a term with '/' as outmost operator and
30 numerator and nominator in normalform [2] or [3].
31 ****************************************************************.*)
37 val add_fraction_ : theory -> term -> (term * term list) option
38 val add_fraction_p_ : theory -> term -> (term * term list) option
39 val calculate_Rational : rls
42 val cancel_ : theory -> term -> (term * term list) option
44 val cancel_p_ : theory -> term -> (term * term list) option
45 val common_nominator : rls
46 val common_nominator_ : theory -> term -> (term * term list) option
47 val common_nominator_p : rls
48 val common_nominator_p_ : theory -> term -> (term * term list) option
49 val eval_is_expanded : string -> 'a -> term -> theory ->
50 (string * term) option
51 val expanded2polynomial : term -> term option
52 val factout_ : theory -> term -> (term * term list) option
53 val factout_p_ : theory -> term -> (term * term list) option
54 val is_expanded : term -> bool
55 val is_polynomial : term -> bool
57 val mv_gcd : (int * int list) list -> mv_poly -> mv_poly
58 val mv_lcm : mv_poly -> mv_poly -> mv_poly
60 val norm_expanded_rat_ : theory -> term -> (term * term list) option
61 (*WN0602.2.6.pull into struct !!!
62 val norm_Rational : rls(*.normalizes an arbitrary rational term without
63 roots into a simple and canceled fraction
64 with normalform [2].*)
66 (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME
67 rls (*.normalizes an rational term [2] without
68 roots into a simple and canceled fraction
69 with normalform [2].*)
71 val norm_rational_ : theory -> term -> (term * term list) option
72 val polynomial2expanded : term -> term option
74 rls (*.evaluates an arbitrary rational term with numerals.*)
76 (*WN0210???SK: fehlen Funktionen, die exportiert werden sollen ? *)
79 (*.**************************************************************************
80 survey on the functions
81 ~~~~~~~~~~~~~~~~~~~~~~~
82 [2] 'polynomial' :rls | [3]'expanded_term':rls
83 --------------------:------------------+-------------------:-----------------
84 factout_p_ : | factout_ :
85 cancel_p_ : | cancel_ :
87 --------------------:------------------+-------------------:-----------------
88 common_nominator_p_: | common_nominator_ :
89 :common_nominator_p| :common_nominator
90 add_fraction_p_ : | add_fraction_ :
91 --------------------:------------------+-------------------:-----------------
92 ???SK :norm_rational_p | :norm_rational
94 This survey shows only the principal functions for reuse, and the identifiers
95 of the rls exported. The list below shows some more useful functions.
98 conversion from Isabelle-term to internal representation
99 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101 ... BITTE FORTSETZEN ...
103 polynomial2expanded = ...
104 expanded2polynomial = ...
106 remark: polynomial2expanded o expanded2polynomial = I,
107 where 'o' is function chaining, and 'I' is identity WN0210???SK
109 functions for greatest common divisor and canceling
110 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
117 functions for least common multiple and addition of fractions
118 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122 add_fraction_ (*.add 2 or more fractions.*)
123 add_fraction_p_ (*.add 2 or more fractions.*)
125 functions for normalform of rationals
126 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127 WN0210???SK interne Funktionen f"ur norm_rational:
128 schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?
133 **************************************************************************.*)
137 structure RationalI : RATIONALI =
141 infix mem ins union; (*WN100819 updating to Isabelle2009-2*)
143 | x mem (y :: ys) = x = y orelse x mem ys;
144 fun (x ins xs) = if x mem xs then xs else x :: xs;
147 | (x :: xs) union ys = xs union (x ins ys);
149 (*. gcd of integers .*)
150 (* die gcd Funktion von Isabelle funktioniert nicht richtig !!! *)
151 fun gcd_int a b = if b=0 then a
152 else gcd_int b (a mod b);
154 (*. univariate polynomials (uv) .*)
155 (*. univariate polynomials are represented as a list of the coefficent in reverse maximum degree order .*)
156 (*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*)
157 type uv_poly = int list;
159 (*. adds two uv polynomials .*)
160 fun uv_mod_add_poly ([]:uv_poly,p2:uv_poly) = p2:uv_poly
161 | uv_mod_add_poly (p1,[]) = p1
162 | uv_mod_add_poly (x::p1,y::p2) = (x+y)::(uv_mod_add_poly(p1,p2));
164 (*. multiplies a uv polynomial with a skalar s .*)
165 fun uv_mod_smul_poly ([]:uv_poly,s:int) = []:uv_poly
166 | uv_mod_smul_poly (x::p,s) = (x*s)::(uv_mod_smul_poly(p,s));
168 (*. calculates the remainder of a polynomial divided by a skalar s .*)
169 fun uv_mod_rem_poly ([]:uv_poly,s) = []:uv_poly
170 | uv_mod_rem_poly (x::p,s) = (x mod s)::(uv_mod_smul_poly(p,s));
172 (*. calculates the degree of a uv polynomial .*)
173 fun uv_mod_deg ([]:uv_poly) = 0
174 | uv_mod_deg p = length(p)-1;
176 (*. calculates the remainder of x/p and represents it as value between -p/2 and p/2 .*)
177 fun uv_mod_mod2(x,p)=
181 if (y)>(p div 2) then (y)-p else
183 if (y)<(~p div 2) then p+(y) else (y)
187 (*.calculates the remainder for each element of a integer list divided by p.*)
188 fun uv_mod_list_modp [] p = []
189 | uv_mod_list_modp (x::xs) p = (uv_mod_mod2(x,p))::(uv_mod_list_modp xs p);
191 (*. appends an integer at the end of a integer list .*)
192 fun uv_mod_null (p1:int list,0) = p1
193 | uv_mod_null (p1:int list,n1:int) = uv_mod_null(p1,n1-1) @ [0];
195 (*. uv polynomial division, result is (quotient, remainder) .*)
196 (*. only for uv_mod_divides .*)
197 (* FIXME: Division von x^9+x^5+1 durch x-1000 funktioniert nicht integer zu klein *)
198 fun uv_mod_pdiv (p1:uv_poly) ([]:uv_poly) = raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")
199 | uv_mod_pdiv p1 [x] =
205 xs:=(uv_mod_rem_poly(p1,x));
206 while length(!xs)>0 andalso hd(!xs)=0 do xs:=tl(!xs)
208 else raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero");
209 ([]:uv_poly,!xs:uv_poly)
211 | uv_mod_pdiv p1 p2 =
213 val n= uv_mod_deg(p2);
214 val m= ref (uv_mod_deg(p1));
215 val p1'=ref (rev(p1));
220 val output=ref ([],[]);
223 if (!m)=0 orelse p2=[0] then raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero")
234 c:=hd(!p1') div hd(p2');
237 p1':=uv_mod_add_poly(!p1',uv_mod_null(uv_mod_smul_poly(p2',~(!c)),!m-n));
238 while length(!p1')>0 andalso hd(!p1')=0 do p1':= tl(!p1');
243 output:=(rev(!q),rev(!p1'))
250 (*. divides p1 by p2 in Zp .*)
251 fun uv_mod_pdivp (p1:uv_poly) (p2:uv_poly) p =
253 val n=uv_mod_deg(p2);
254 val m=ref (uv_mod_deg(uv_mod_list_modp p1 p));
255 val p1'=ref (rev(p1));
256 val p2'=(rev(uv_mod_list_modp p2 p));
260 val output=ref ([],[]);
263 if (!m)=0 orelse p2=[0] then raise error ("RATIONALS_UV_MOD_PDIVP_EXCEPTION: Division by zero")
274 c:=uv_mod_mod2(hd(!p1')*(power lc2 1), p);
276 p1':=uv_mod_list_modp(tl(uv_mod_add_poly(uv_mod_smul_poly(!p1',lc2),
277 uv_mod_smul_poly(uv_mod_smul_poly(p2',hd(!p1')),~1)))) p;
281 while !p1'<>[] andalso hd(!p1')=0 do
286 output:=(rev(uv_mod_list_modp (!q) (p)),rev(!p1'))
289 !output:uv_poly * uv_poly
293 (*. calculates the remainder of p1/p2 .*)
294 fun uv_mod_prest (p1:uv_poly) ([]:uv_poly) = raise error("UV_MOD_PREST_EXCEPTION: Division by zero")
295 | uv_mod_prest [] p2 = []:uv_poly
296 | uv_mod_prest p1 p2 = (#2(uv_mod_pdiv p1 p2));
298 (*. calculates the remainder of p1/p2 in Zp .*)
299 fun uv_mod_prestp (p1:uv_poly) ([]:uv_poly) p= raise error("UV_MOD_PRESTP_EXCEPTION: Division by zero")
300 | uv_mod_prestp [] p2 p= []:uv_poly
301 | uv_mod_prestp p1 p2 p = #2(uv_mod_pdivp p1 p2 p);
303 (*. calculates the content of a uv polynomial .*)
304 fun uv_mod_cont ([]:uv_poly) = 0
305 | uv_mod_cont (x::p)= gcd_int x (uv_mod_cont(p));
307 (*. divides each coefficient of a uv polynomial by y .*)
308 fun uv_mod_div_list (p:uv_poly,0) = raise error("UV_MOD_DIV_LIST_EXCEPTION: Division by zero")
309 | uv_mod_div_list ([],y) = []:uv_poly
310 | uv_mod_div_list (x::p,y) = (x div y)::uv_mod_div_list(p,y);
312 (*. calculates the primitiv part of a uv polynomial .*)
313 fun uv_mod_pp ([]:uv_poly) = []:uv_poly
321 if !c=0 then raise error ("RATIONALS_UV_MOD_PP_EXCEPTION: content is 0")
322 else uv_mod_div_list(p,!c)
326 (*. gets the leading coefficient of a uv polynomial .*)
327 fun uv_mod_lc ([]:uv_poly) = 0
328 | uv_mod_lc p = hd(rev(p));
330 (*. calculates the euklidean polynomial remainder sequence in Zp .*)
331 fun uv_mod_prs_euklid_p(p1:uv_poly,p2:uv_poly,p)=
339 while uv_mod_deg(!f')>0 do
341 f':=uv_mod_prestp (hd(tl(!f))) (hd(!f)) p;
354 (*. calculates the gcd of p1 and p2 in Zp .*)
355 fun uv_mod_gcd_modp ([]:uv_poly) (p2:uv_poly) p = p2:uv_poly
356 | uv_mod_gcd_modp p1 [] p= p1
357 | uv_mod_gcd_modp p1 p2 p=
367 if uv_mod_deg(p1)>=uv_mod_deg(p2) then
369 p1':=uv_mod_list_modp (uv_mod_pp(p1)) p;
370 p2':=uv_mod_list_modp (uv_mod_pp(p2)) p
374 p1':=uv_mod_list_modp (uv_mod_pp(p2)) p;
375 p2':=uv_mod_list_modp (uv_mod_pp(p1)) p
377 d:=uv_mod_mod2((gcd_int (uv_mod_cont(p1))) (uv_mod_cont(p2)), p) ;
378 if !d>(p div 2) then d:=(!d)-p else ();
380 prs:=uv_mod_prs_euklid_p(!p1',!p2',p);
382 if hd(!prs)=[] then pc:=hd(tl(!prs))
385 g:=uv_mod_smul_poly(uv_mod_pp(!pc),!d);
390 (*. calculates the minimum of two real values x and y .*)
391 fun uv_mod_r_min(x,y):BasisLibrary.Real.real = if x>y then y else x;
393 (*. calculates the minimum of two integer values x and y .*)
394 fun uv_mod_min(x,y) = if x>y then y else x;
396 (*. adds the squared values of a integer list .*)
397 fun uv_mod_add_qu [] = 0.0
398 | uv_mod_add_qu (x::p) = BasisLibrary.Real.fromInt(x)*BasisLibrary.Real.fromInt(x) + uv_mod_add_qu p;
400 (*. calculates the euklidean norm .*)
401 fun uv_mod_norm ([]:uv_poly) = 0.0
402 | uv_mod_norm p = Math.sqrt(uv_mod_add_qu(p));
404 (*. multipies two values a and b .*)
405 fun uv_mod_multi a b = a * b;
407 (*. decides if x is a prim, the list contains all primes which are lower then x .*)
408 fun uv_mod_prim(x,[])= false
409 | uv_mod_prim(x,[y])=if ((x mod y) <> 0) then true
411 | uv_mod_prim(x,y::ys) = if uv_mod_prim(x,[y])
413 if uv_mod_prim(x,ys) then true
417 (*. gets the first prime, which is greater than p and does not divide g .*)
418 fun uv_mod_nextprime(g,p)=
424 while (!i<p) do (* calculates the primes lower then p *)
426 if uv_mod_prim(!i,!list) then
431 list:= (!i)::(!list);
439 while (!exit=0) do (* calculate next prime which does not divide g *)
441 if uv_mod_prim(!i,!list) then
446 list:= (!i)::(!list);
456 (*. decides if p1 is a factor of p2 in Zp .*)
457 fun uv_mod_dividesp ([]:uv_poly) (p2:uv_poly) p= raise error("UV_MOD_DIVIDESP: Division by zero")
458 | uv_mod_dividesp p1 p2 p= if uv_mod_prestp p2 p1 p = [] then true else false;
460 (*. decides if p1 is a factor of p2 .*)
461 fun uv_mod_divides ([]:uv_poly) (p2:uv_poly) = raise error("UV_MOD_DIVIDES: Division by zero")
462 | uv_mod_divides p1 p2 = if uv_mod_prest p2 p1 = [] then true else false;
464 (*. chinese remainder algorithm .*)
465 fun uv_mod_cra2(r1,r2,m1,m2)=
473 while (uv_mod_mod2((!c)*m1,m2))<>1 do
477 r1':= uv_mod_mod2(r1,m1);
478 d:=uv_mod_mod2(((r2-(!r1'))*(!c)),m2);
483 (*. applies the chinese remainder algorithmen to the coefficients of x1 and x2 .*)
484 fun uv_mod_cra_2 ([],[],m1,m2) = []
485 | uv_mod_cra_2 ([],x2,m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x1")
486 | uv_mod_cra_2 (x1,[],m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x2")
487 | 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));
489 (*. calculates the gcd of two uv polynomials p1' and p2' with the modular algorithm .*)
490 fun uv_mod_gcd (p1':uv_poly) (p2':uv_poly) =
492 val p1=ref (uv_mod_pp(p1'));
493 val p2=ref (uv_mod_pp(p2'));
494 val c=gcd_int (uv_mod_cont(p1')) (uv_mod_cont(p2'));
507 if length(!p1)>length(!p2) then ()
516 d:=gcd_int (uv_mod_lc(!p1)) (uv_mod_lc(!p2));
517 g1:=uv_mod_lc(!p1)*uv_mod_lc(!p2);
520 m:=BasisLibrary.Real.ceil(2.0 *
521 BasisLibrary.Real.fromInt(!d) *
522 BasisLibrary.Real.fromInt(power 2 (uv_mod_min(uv_mod_deg(!p2),uv_mod_deg(!p1)))) *
523 BasisLibrary.Real.fromInt(!d) *
524 uv_mod_r_min(uv_mod_norm(!p1) / BasisLibrary.Real.fromInt(abs(uv_mod_lc(!p1))),
525 uv_mod_norm(!p2) / BasisLibrary.Real.fromInt(abs(uv_mod_lc(!p2)))));
529 p:=uv_mod_nextprime(!d,!p);
530 cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)) ;
531 if abs(uv_mod_lc(!cp))<>1 then (* leading coefficient = 1 ? *)
534 while (!i)<(!p) andalso (abs(uv_mod_mod2((uv_mod_lc(!cp)*(!i)),(!p)))<>1) do
538 cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)
542 qp:= ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp));
544 if uv_mod_deg(!qp)=0 then (q:=[1]; exit:=1) else ();
549 while !pn<= !m andalso !m>(!p) andalso !exit=0 do
551 p:=uv_mod_nextprime(!d,!p);
552 cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p));
553 if uv_mod_lc(!cp)<>1 then (* leading coefficient = 1 ? *)
556 while (!i)<(!p) andalso ((uv_mod_mod2((uv_mod_lc(!q)*(!i)),(!p)))<>1) do
560 cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)
564 qp:=uv_mod_list_modp ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp) ) (!p);
565 if uv_mod_deg(!qp)>uv_mod_deg(!q) then
567 (*print("degree to high!!!\n")*)
571 if uv_mod_deg(!qp)=uv_mod_deg(!q) then
573 q:=uv_mod_cra_2(!q,!qp,!pn,!p);
575 q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn)); (* found already gcd ? *)
576 if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then (exit:=1) else ()
580 if uv_mod_deg(!qp)<uv_mod_deg(!q) then
589 q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn));
590 if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then exit:=1 else ()
592 uv_mod_smul_poly(!q,c):uv_poly
595 (*. multivariate polynomials .*)
596 (*. multivariate polynomials are represented as a list of the pairs,
597 first is the coefficent and the second is a list of the exponents .*)
598 (*. 5 * x^5 * y^3 + 4 * x^3 * z^2 + 2 * x^2 * y * z^3 - z - 19
599 => [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*)
601 (*. global variables .*)
602 (*. order indicators .*)
603 val LEX_=0; (* lexicographical term order *)
604 val GGO_=1; (* greatest degree order *)
606 (*. datatypes for internal representation.*)
607 type mv_monom = (int * (*.coefficient or the monom.*)
608 int list); (*.list of exponents) .*)
609 fun mv_monom2str (i, is) = "("^ int2str i^"," ^ ints2str' is ^ ")";
611 type mv_poly = mv_monom list;
612 fun mv_poly2str p = (strs2str' o (map mv_monom2str)) p;
614 (*. help function for monom_greater and geq .*)
615 fun mv_mg_hlp([]) = EQUAL
616 | mv_mg_hlp(x::list)=if x<0 then LESS
617 else if x>0 then GREATER
618 else mv_mg_hlp(list);
620 (*. adds a list of values .*)
621 fun mv_addlist([]) = 0
622 | mv_addlist(p1) = hd(p1)+mv_addlist(tl(p1));
624 (*. tests if the monomial M1 is greater as the monomial M2 and returns a boolean value .*)
625 (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
626 fun mv_monom_greater((M1x,M1l):mv_monom,(M2x,M2l):mv_monom,order)=
629 if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
630 else if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
635 if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
637 if mv_addlist(M1l)=mv_addlist(M2l) then if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
638 else if mv_addlist(M1l)>mv_addlist(M2l) then true else false
640 else raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");
642 (*. tests if the monomial X is greater as the monomial Y and returns a order value (GREATER,EQUAL,LESS) .*)
643 (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
644 fun mv_geq order ((x1,x):mv_monom,(x2,y):mv_monom) =
650 if length(x)<>length(y) then
651 raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
654 temp:=mv_mg_hlp((map op- (x~~y)));
656 ( if x1=x2 then EQUAL
657 else if x1>x2 then GREATER
666 if length(x)<>length(y) then
667 raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
669 if mv_addlist(x)=mv_addlist(y) then
670 (mv_mg_hlp((map op- (x~~y))))
671 else if mv_addlist(x)>mv_addlist(y) then GREATER else LESS
673 else raise error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")
676 (*. cuts the first variable from a polynomial .*)
677 fun mv_cut([]:mv_poly)=[]:mv_poly
678 | mv_cut((x,[])::list) = raise error ("RATIONALS_MV_CUT_EXCEPTION: Invalid list ")
679 | mv_cut((x,y::ys)::list)=(x,ys)::mv_cut(list);
681 (*. leading power product .*)
682 fun mv_lpp([]:mv_poly,order) = []
683 | mv_lpp([(x,y)],order) = y
684 | mv_lpp(p1,order) = #2(hd(rev(sort (mv_geq order) p1)));
686 (*. leading monomial .*)
687 fun mv_lm([]:mv_poly,order) = (0,[]):mv_monom
688 | mv_lm([x],order) = x
689 | mv_lm(p1,order) = hd(rev(sort (mv_geq order) p1));
691 (*. leading coefficient in term order .*)
692 fun mv_lc2([]:mv_poly,order) = 0
693 | mv_lc2([(x,y)],order) = x
694 | mv_lc2(p1,order) = #1(hd(rev(sort (mv_geq order) p1)));
697 (*. reverse the coefficients in mv polynomial .*)
698 fun mv_rev_to([]:mv_poly) = []:mv_poly
699 | mv_rev_to((c,e)::xs) = (c,rev(e))::mv_rev_to(xs);
701 (*. leading coefficient in reverse term order .*)
702 fun mv_lc([]:mv_poly,order) = []:mv_poly
703 | mv_lc([(x,y)],order) = mv_rev_to(mv_cut(mv_rev_to([(x,y)])))
706 val p1o=ref (rev(sort (mv_geq order) (mv_rev_to(p1))));
707 val lp=hd(#2(hd(!p1o)));
711 while (length(!p1o)>0 andalso hd(#2(hd(!p1o)))=lp) do
713 lc:=hd(mv_cut([hd(!p1o)]))::(!lc);
716 if !lc=[] then raise error ("RATIONALS_MV_LC_EXCEPTION: lc is empty") else ();
721 (*. compares two powerproducts .*)
722 fun mv_monom_equal((_,xlist):mv_monom,(_,ylist):mv_monom) = (foldr and_) (((map op=) (xlist~~ylist)),true);
724 (*. help function for mv_add .*)
725 fun mv_madd([]:mv_poly,[]:mv_poly,order) = []:mv_poly
726 | mv_madd([(0,_)],p2,order) = p2
727 | mv_madd(p1,[(0,_)],order) = p1
728 | mv_madd([],p2,order) = p2
729 | mv_madd(p1,[],order) = p1
730 | mv_madd(p1,p2,order) =
732 if mv_monom_greater(hd(p1),hd(p2),order)
733 then hd(p1)::mv_madd(tl(p1),p2,order)
734 else if mv_monom_equal(hd(p1),hd(p2))
735 then if mv_lc2(p1,order)+mv_lc2(p2,order)<>0
736 then (mv_lc2(p1,order)+mv_lc2(p2,order),mv_lpp(p1,order))::mv_madd(tl(p1),tl(p2),order)
737 else mv_madd(tl(p1),tl(p2),order)
738 else hd(p2)::mv_madd(p1,tl(p2),order)
741 (*. adds two multivariate polynomials .*)
742 fun mv_add([]:mv_poly,p2:mv_poly,order) = p2
743 | mv_add(p1,[],order) = p1
744 | mv_add(p1,p2,order) = mv_madd(rev(sort (mv_geq order) p1),rev(sort (mv_geq order) p2), order);
746 (*. monom multiplication .*)
747 fun mv_mmul((x1,y1):mv_monom,(x2,y2):mv_monom)=(x1*x2,(map op+) (y1~~y2)):mv_monom;
749 (*. deletes all monomials with coefficient 0 .*)
750 fun mv_shorten([]:mv_poly,order) = []:mv_poly
751 | mv_shorten(x::xs,order)=mv_madd([x],mv_shorten(xs,order),order);
755 | mv_null2(x::l)=0::mv_null2(l);
757 (*. multiplies two multivariate polynomials .*)
758 fun mv_mul([]:mv_poly,[]:mv_poly,_) = []:mv_poly
759 | mv_mul([],y::p2,_) = [(0,mv_null2(#2(y)))]
760 | mv_mul(x::p1,[],_) = [(0,mv_null2(#2(x)))]
761 | mv_mul(x::p1,y::p2,order) = mv_shorten(rev(sort (mv_geq order) (mv_mmul(x,y) :: (mv_mul(p1,y::p2,order) @
762 mv_mul([x],p2,order)))),order);
764 (*. gets the maximum value of a list .*)
766 | mv_getmax(x::p1)= let
772 (*. calculates the maximum degree of an multivariate polynomial .*)
773 fun mv_grad([]:mv_poly) = 0
774 | mv_grad(p1:mv_poly)= mv_getmax((map mv_addlist) ((map #2) p1));
776 (*. converts the sign of a value .*)
777 fun mv_minus(x)=(~1) * x;
779 (*. converts the sign of all coefficients of a polynomial .*)
780 fun mv_minus2([]:mv_poly)=[]:mv_poly
781 | mv_minus2(p1)=(mv_minus(#1(hd(p1))),#2(hd(p1)))::(mv_minus2(tl(p1)));
783 (*. searches for a negativ value in a list .*)
784 fun mv_is_negativ([])=false
785 | mv_is_negativ(x::xs)=if x<0 then true else mv_is_negativ(xs);
787 (*. division of monomials .*)
788 fun mv_mdiv((0,[]):mv_monom,_:mv_monom)=(0,[]):mv_monom
789 | mv_mdiv(_,(0,[]))= raise error ("RATIONALS_MV_MDIV_EXCEPTION Division by 0 ")
790 | mv_mdiv(p1:mv_monom,p2:mv_monom)=
796 if !c=0 then raise error("MV_MDIV_EXCEPTION Dividing by zero")
797 else c:=(#1(p1) div #1(p2));
800 pp:=(#2(mv_mmul((1,#2(p1)),(1,(map mv_minus) (#2(p2))))));
801 if mv_is_negativ(!pp) then (0,!pp)
804 else raise error("MV_MDIV_EXCEPTION Dividing by empty Polynom")
808 (*. prints a polynom for (internal use only) .*)
809 fun mv_print_poly([]:mv_poly)=print("[]\n")
810 | mv_print_poly((x,y)::[])= print("("^BasisLibrary.Int.toString(x)^","^ints2str(y)^")\n")
811 | mv_print_poly((x,y)::p1) = (print("("^BasisLibrary.Int.toString(x)^","^ints2str(y)^"),");mv_print_poly(p1));
814 (*. division of two multivariate polynomials .*)
815 fun mv_division([]:mv_poly,g:mv_poly,order)=([]:mv_poly,[]:mv_poly)
816 | mv_division(f,[],order)= raise error ("RATIONALS_MV_DIVISION_EXCEPTION Division by zero")
817 | mv_division(f,g,order)=
826 r := rev(sort (mv_geq order) (mv_shorten(f,order)));
827 g':= rev(sort (mv_geq order) (mv_shorten(g,order)));
828 if #1(hd(!g'))=0 then raise error("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero") else ();
829 if (mv_monom_greater (hd(!g'),hd(!r),order)) then ([(0,mv_null2(#2(hd(f))))],(!r))
833 while (if (!exit)=0 then not(mv_monom_greater (hd(!g'),hd(!r),order)) else false) do
835 if (#1(mv_lm(!g',order)))<>0 then m:=mv_mdiv(mv_lm(!r,order),mv_lm(!g',order))
836 else raise error ("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero");
840 r:=mv_add((!r),mv_minus2(mv_mul(!g',[!m],order)),order)
843 if (if length(!r)<>0 then length(!g')<>0 else false) then ()
850 (*. multiplies a polynomial with an integer .*)
851 fun mv_skalar_mul([]:mv_poly,c) = []:mv_poly
852 | mv_skalar_mul((x,y)::p1,c) = ((x * c),y)::mv_skalar_mul(p1,c);
854 (*. inserts the a first variable into an polynomial with exponent v .*)
855 fun mv_correct([]:mv_poly,v:int)=[]:mv_poly
856 | mv_correct((x,y)::list,v:int)=(x,v::y)::mv_correct(list,v);
858 (*. multivariate case .*)
860 (*. decides if x is a factor of y .*)
861 fun mv_divides([]:mv_poly,[]:mv_poly)= raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
862 | mv_divides(x,[]) = raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
863 | mv_divides(x:mv_poly,y:mv_poly) = #2(mv_division(y,x,LEX_))=[];
865 (*. gets the maximum of a and b .*)
866 fun mv_max(a,b) = if a>b then a else b;
868 (*. gets the maximum exponent of a mv polynomial in the lexicographic term order .*)
869 fun mv_deg([]:mv_poly) = 0
872 val p1'=mv_shorten(p1,LEX_);
874 if length(p1')=0 then 0
875 else mv_max(hd(#2(hd(p1'))),mv_deg(tl(p1')))
878 (*. gets the maximum exponent of a mv polynomial in the reverse lexicographic term order .*)
879 fun mv_deg2([]:mv_poly) = 0
882 val p1'=mv_shorten(p1,LEX_);
884 if length(p1')=0 then 0
885 else mv_max(hd(rev(#2(hd(p1')))),mv_deg2(tl(p1')))
888 (*. evaluates the mv polynomial at the value v of the main variable .*)
889 fun mv_subs([]:mv_poly,v) = []:mv_poly
890 | mv_subs((c,e)::p1:mv_poly,v) = mv_skalar_mul(mv_cut([(c,e)]),power v (hd(e))) @ mv_subs(p1,v);
892 (*. calculates the content of a uv-polynomial in mv-representation .*)
893 fun uv_content2([]:mv_poly) = 0
894 | uv_content2((c,e)::p1) = (gcd_int c (uv_content2(p1)));
896 (*. converts a uv-polynomial from mv-representation to uv-representation .*)
897 fun uv_to_list ([]:mv_poly)=[]:uv_poly
898 | uv_to_list ((c1,e1)::others) =
901 val max=mv_grad((c1,e1)::others);
902 val help=ref ((c1,e1)::others);
905 if length(e1)>1 then raise error ("RATIONALS_TO_LIST_EXCEPTION: not univariate")
906 else if length(e1)=0 then [c1]
910 while (!count)<=max do
912 if length(!help)>0 andalso hd(#2(hd(!help)))=max-(!count) then
914 list:=(#1(hd(!help)))::(!list);
921 count := (!count) + 1
927 (*. converts a uv-polynomial from uv-representation to mv-representation .*)
928 fun uv_to_poly ([]:uv_poly) = []:mv_poly
935 while length(!help)>0 do
937 if hd(!help)=0 then ()
938 else list:=(hd(!help),[!count])::(!list);
945 (*. univariate gcd calculation from polynomials in multivariate representation .*)
946 fun uv_gcd ([]:mv_poly) p2 = p2
947 | uv_gcd p1 ([]:mv_poly) = p1
948 | uv_gcd p1 [(c,[e])] =
950 val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
951 val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
953 [(gcd_int (uv_content2(p1)) c,[min])]
955 | uv_gcd [(c,[e])] p2 =
957 val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p2,LEX_))));
958 val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
960 [(gcd_int (uv_content2(p2)) c,[min])]
962 | uv_gcd p11 p22 = uv_to_poly(uv_mod_gcd (uv_to_list(mv_shorten(p11,LEX_))) (uv_to_list(mv_shorten(p22,LEX_))));
964 (*. help function for the newton interpolation .*)
965 fun mv_newton_help ([]:mv_poly list,k:int) = []:mv_poly list
966 | mv_newton_help (pl:mv_poly list,k) =
975 while length(!x)>1 do
977 if length(hd(!x))>0 then n1:=mv_null2(#2(hd(hd(!x))))
978 else if length(hd(tl(!x)))>0 then n1:=mv_null2(#2(hd(hd(tl(!x)))))
980 t:= #1(mv_division(mv_add(hd(!x),mv_skalar_mul(hd(tl(!x)),~1),LEX_),[(k,!n1)],LEX_));
988 (*. help function for the newton interpolation .*)
989 fun mv_newton_add ([]:mv_poly list) t = []:mv_poly
990 | mv_newton_add [x:mv_poly] t = x
991 | mv_newton_add (pl:mv_poly list) t =
998 while length(!pll)>0 andalso hd(!pll)=[] do
1002 if length(!pll)>0 then expos:= #2(hd(hd(!pll))) else expos:=[];
1005 mv_add(mv_correct(mv_cut([(1,mv_null2(!expos))]),1),[(~t,mv_null2(!expos))],LEX_),
1006 mv_newton_add (tl(pl)) (t+1),
1013 (*. calculates the newton interpolation with polynomial coefficients .*)
1014 (*. step-depth is 1 and if the result is not an integerpolynomial .*)
1015 (*. this function returns [] .*)
1016 fun mv_newton ([]:(mv_poly) list) = []:mv_poly
1017 | mv_newton ([mp]:(mv_poly) list) = mp:mv_poly
1028 c:=mv_newton_help(!c,!k);
1029 c1:=(hd(!c))::(!c1);
1030 while(length(!c)>1 andalso !k<n) do
1033 while length(!c)>0 andalso hd(!c)=[] do c:=tl(!c);
1034 if !c=[] then () else c:=mv_newton_help(!c,!k);
1036 if !c=[] then () else c1:=(hd(!c))::(!c1)
1038 while hd(!c1)=[] do c1:=tl(!c1);
1041 mv_newton_add (!c1) 1
1044 (*. sets the exponents of the first variable to zero .*)
1045 fun mv_null3([]:mv_poly) = []:mv_poly
1046 | mv_null3((x,y)::xs) = (x,0::tl(y))::mv_null3(xs);
1048 (*. calculates the minimum exponents of a multivariate polynomial .*)
1049 fun mv_min_pp([]:mv_poly)=[]
1050 | mv_min_pp((c,e)::xs)=
1057 while length(!y)>0 do
1059 x:=(map uv_mod_min) ((!x) ~~ (#2(hd(!y))));
1066 (*. checks if all elements of the list have value zero .*)
1067 fun list_is_null [] = true
1068 | list_is_null (x::xs) = (x=0 andalso list_is_null(xs));
1070 (* check if main variable is zero*)
1071 fun main_zero (ms : mv_poly) = (list_is_null o (map (hd o #2))) ms;
1073 (*. calculates the content of an polynomial .*)
1074 fun mv_content([]:mv_poly) = []:mv_poly
1077 val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
1078 val test=ref (hd(#2(hd(!list))));
1080 val min=(hd(#2(hd(rev(!list)))));
1083 if length(!list)>1 then
1085 while (if length(!list)>0 then (hd(#2(hd(!list)))=(!test)) else false) do
1087 result:=(#1(hd(!list)),tl(#2(hd(!list))))::(!result);
1089 if length(!list)<1 then list:=[]
1090 else list:=tl(!list)
1093 if length(!list)>0 then
1095 list:=mv_gcd (!result) (mv_cut(mv_content(!list)))
1097 else list:=(!result);
1098 list:=mv_correct(!list,0);
1108 (*. calculates the primitiv part of a polynomial .*)
1109 and mv_pp([]:mv_poly) = []:mv_poly
1114 cont:=mv_content(p1);
1115 pp:=(#1(mv_division(p1,!cont,LEX_)));
1117 then raise error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")
1121 (*. calculates the gcd of two multivariate polynomials with a modular approach .*)
1122 and mv_gcd ([]:mv_poly) ([]:mv_poly) :mv_poly= []:mv_poly
1123 | mv_gcd ([]:mv_poly) (p2) :mv_poly= p2:mv_poly
1124 | mv_gcd (p1:mv_poly) ([]) :mv_poly= p1:mv_poly
1125 | mv_gcd ([(x,xs)]:mv_poly) ([(y,ys)]):mv_poly =
1127 val xpoly:mv_poly = [(x,xs)];
1128 val ypoly:mv_poly = [(y,ys)];
1131 if xs=ys then [((gcd_int x y),xs)]
1132 else [((gcd_int x y),(map uv_mod_min)(xs~~ys))]:mv_poly
1135 | mv_gcd (p1:mv_poly) ([(y,ys)]) :mv_poly=
1137 [(gcd_int (uv_content2(p1)) (y),(map uv_mod_min)(mv_min_pp(p1)~~ys))]:mv_poly
1139 | mv_gcd ([(y,ys)]:mv_poly) (p2):mv_poly =
1141 [(gcd_int (uv_content2(p2)) (y),(map uv_mod_min)(mv_min_pp(p2)~~ys))]:mv_poly
1143 | mv_gcd (p1':mv_poly) (p2':mv_poly):mv_poly=
1145 val vc=length(#2(hd(p1')));
1148 if main_zero(mv_content(p1')) andalso
1149 (main_zero(mv_content(p2'))) then
1150 mv_correct((mv_gcd (mv_cut(mv_content(p1'))) (mv_cut(mv_content(p2')))),0)
1152 mv_gcd (mv_content(p1')) (mv_content(p2'))
1154 val p1= #1(mv_division(p1',mv_content(p1'),LEX_));
1155 val p2= #1(mv_division(p2',mv_content(p2'),LEX_));
1157 val candidate=ref [];
1158 val interpolation_list=ref [];
1169 val current_degree=ref 99999; (*. FIXME: unlimited ! .*)
1172 if vc<2 then (* areUnivariate(p1',p2') *)
1174 gcd:=uv_gcd (mv_shorten(p1',LEX_)) (mv_shorten(p2',LEX_))
1181 p1r := mv_lc(p1,LEX_);
1182 p2r := mv_lc(p2,LEX_);
1183 if main_zero(!p1r) andalso
1187 delta := mv_correct((mv_gcd (mv_cut (!p1r)) (mv_cut (!p2r))),0)
1191 delta := mv_gcd (!p1r) (!p2r)
1193 (*if mv_shorten(mv_subs(!p1r,!r),LEX_)=[] andalso
1194 mv_shorten(mv_subs(!p2r,!r),LEX_)=[] *)
1195 if mv_lc2(mv_shorten(mv_subs(!p1r,!r),LEX_),LEX_)=0 andalso
1196 mv_lc2(mv_shorten(mv_subs(!p2r,!r),LEX_),LEX_)=0
1202 gcd_r:=mv_shorten(mv_gcd (mv_shorten(mv_subs(p1,!r),LEX_))
1203 (mv_shorten(mv_subs(p2,!r),LEX_)) ,LEX_);
1204 gcd_r:= #1(mv_division(mv_mul(mv_correct(mv_subs(!delta,!r),0),!gcd_r,LEX_),
1205 mv_correct(mv_lc(!gcd_r,LEX_),0),LEX_));
1206 d:=mv_deg2(!gcd_r); (* deg(gcd_r,z) *)
1207 if (!d < !current_degree) then
1209 current_degree:= !d;
1210 interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
1214 if (!d = !current_degree) then
1216 interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
1221 if length(!interpolation_list)> uv_mod_min(mv_deg(p1),mv_deg(p2)) then
1223 candidate := mv_newton(rev(!interpolation_list));
1224 if !candidate=[] then ()
1227 candidate:=mv_pp(!candidate);
1228 if mv_divides(!candidate,p1) andalso mv_divides(!candidate,p2) then
1230 gcd:= mv_mul(!candidate,cont,LEX_);
1235 interpolation_list:=[mv_correct(!gcd_r,0)]
1245 (*. calculates the least common divisor of two polynomials .*)
1246 fun mv_lcm (p1:mv_poly) (p2:mv_poly) :mv_poly =
1248 #1(mv_division(mv_mul(p1,p2,LEX_),mv_gcd p1 p2,LEX_))
1251 (*. gets the variables (strings) of a term .*)
1252 fun get_vars(term1) = (map free2str) (vars term1); (*["a","b","c"]; *)
1254 (*. counts the negative coefficents in a polynomial .*)
1255 fun count_neg ([]:mv_poly) = 0
1256 | count_neg ((c,e)::xs) = if c<0 then 1+count_neg xs
1259 (*. help function for is_polynomial
1260 checks the order of the operators .*)
1261 fun test_polynomial (Const ("uminus",_) $ Free (str,_)) _ = true (*WN.13.3.03*)
1262 | test_polynomial (t as Free(str,_)) v = true
1263 | test_polynomial (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false
1264 else (test_polynomial t1 "*") andalso (test_polynomial t2 "*")
1265 | test_polynomial (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1266 else (test_polynomial t1 " ") andalso (test_polynomial t2 " ")
1267 | test_polynomial (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_polynomial t1 "^") andalso (test_polynomial t2 "^")
1268 | test_polynomial _ v = false;
1270 (*. tests if a term is a polynomial .*)
1271 fun is_polynomial t = test_polynomial t " ";
1273 (*. help function for is_expanded
1274 checks the order of the operators .*)
1275 fun test_exp (t as Free(str,_)) v = true
1276 | test_exp (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false
1277 else (test_exp t1 "*") andalso (test_exp t2 "*")
1278 | test_exp (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1279 else (test_exp t1 " ") andalso (test_exp t2 " ")
1280 | test_exp (t as Const ("op -",_) $ t1 $ t2) v = if v="*" orelse v="^" then false
1281 else (test_exp t1 " ") andalso (test_exp t2 " ")
1282 | test_exp (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_exp t1 "^") andalso (test_exp t2 "^")
1283 | test_exp _ v = false;
1286 (*. help function for check_coeff:
1287 converts the term to a list of coefficients .*)
1288 fun term2coef' (t as Free(str,_(*typ*))) v :mv_poly option =
1296 if is_numeral str then
1298 SOME [(((the o int_of_str) str),mv_null2(v))] handle _ => NONE
1304 while ((!len)>(!i)) do
1306 if str=hd((!vh)) then
1317 SOME [(1,rev(!vl))] handle _ => NONE
1320 | term2coef' (Const ("op *",_) $ t1 $ t2) v :mv_poly option=
1328 t1pp:=(#2(hd(the(term2coef' t1 v))));
1329 t2pp:=(#2(hd(the(term2coef' t2 v))));
1330 t1c:=(#1(hd(the(term2coef' t1 v))));
1331 t2c:=(#1(hd(the(term2coef' t2 v))));
1333 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] handle _ => NONE
1337 | term2coef' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ (t2 as Free (str2,_))) v :mv_poly option=
1347 if (not o is_numeral) str1 andalso is_numeral str2 then
1352 while ((!len)>(!i)) do
1354 if str1=hd((!vh)) then
1356 vl:=((the o int_of_str) str2)::(!vl)
1365 SOME [(1,rev(!vl))] handle _ => NONE
1367 else raise error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term")
1370 | term2coef' (Const ("op +",_) $ t1 $ t2) v :mv_poly option=
1372 SOME ((the(term2coef' t1 v)) @ (the(term2coef' t2 v))) handle _ => NONE
1374 | term2coef' (Const ("op -",_) $ t1 $ t2) v :mv_poly option=
1376 SOME ((the(term2coef' t1 v)) @ mv_skalar_mul((the(term2coef' t2 v)),1)) handle _ => NONE
1378 | term2coef' (term) v = raise error ("RATIONALS_TERM2COEF_EXCEPTION 2: Invalid term");
1380 (*. checks if all coefficients of a polynomial are positiv (except the first) .*)
1381 fun check_coeff t = (* erste Koeffizient kann <0 sein !!! *)
1382 if count_neg(tl(the(term2coef' t (get_vars(t)))))=0 then true
1385 (*. checks for expanded term [3] .*)
1386 fun is_expanded t = test_exp t " " andalso check_coeff(t);
1388 (*WN.7.3.03 Hilfsfunktion f"ur term2poly'*)
1389 fun mk_monom v' p vs =
1390 let fun conv p (v: string) = if v'= v then p else 0
1391 in map (conv p) vs end;
1392 (* mk_monom "y" 5 ["a","b","x","y","z"];
1393 val it = [0,0,0,5,0] : int list*)
1395 (*. this function converts the term representation into the internal representation mv_poly .*)
1396 fun term2poly' (Const ("uminus",_) $ Free (str,_)) v = (*WN.7.3.03*)
1398 then SOME [((the o int_of_str) ("-"^str), mk_monom "#" 0 v)]
1399 else SOME [(~1, mk_monom str 1 v)]
1401 | term2poly' (Free(str,_)) v :mv_poly option =
1409 if is_numeral str then
1411 SOME [(((the o int_of_str) str),mv_null2 v)] handle _ => NONE
1417 while ((!len)>(!i)) do
1419 if str=hd((!vh)) then
1430 SOME [(1,rev(!vl))] handle _ => NONE
1433 | term2poly' (Const ("op *",_) $ t1 $ t2) v :mv_poly option=
1441 t1pp:=(#2(hd(the(term2poly' t1 v))));
1442 t2pp:=(#2(hd(the(term2poly' t2 v))));
1443 t1c:=(#1(hd(the(term2poly' t1 v))));
1444 t2c:=(#1(hd(the(term2poly' t2 v))));
1446 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )]
1451 | term2poly' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $
1452 (t2 as Free (str2,_))) v :mv_poly option=
1462 if (not o is_numeral) str1 andalso is_numeral str2 then
1467 while ((!len)>(!i)) do
1469 if str1=hd((!vh)) then
1471 vl:=((the o int_of_str) str2)::(!vl)
1480 SOME [(1,rev(!vl))] handle _ => NONE
1482 else raise error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term")
1485 | term2poly' (Const ("op +",_) $ t1 $ t2) v :mv_poly option =
1487 SOME ((the(term2poly' t1 v)) @ (the(term2poly' t2 v))) handle _ => NONE
1489 | term2poly' (Const ("op -",_) $ t1 $ t2) v :mv_poly option =
1491 SOME ((the(term2poly' t1 v)) @ mv_skalar_mul((the(term2poly' t2 v)),~1)) handle _ => NONE
1493 | term2poly' (term) v = raise error ("RATIONALS_TERM2POLY_EXCEPTION 2: Invalid term");
1495 (*. translates an Isabelle term into internal representation.
1497 fn : term -> (*normalform [2] *)
1498 string list -> (*for ...!!! BITTE DIE ERKLÄRUNG,
1499 DIE DU MIR LETZTES MAL GEGEBEN HAST*)
1500 mv_monom list (*internal representation *)
1501 option (*the translation may fail with NONE*)
1503 fun term2poly (t:term) v =
1504 if is_polynomial t then term2poly' t v
1505 else raise error ("term2poly: invalid = "^(term2str t));
1507 (*. same as term2poly with automatic detection of the variables .*)
1508 fun term2polyx t = term2poly t (((map free2str) o vars) t);
1510 (*. checks if the term is in expanded polynomial form and converts it into the internal representation .*)
1511 fun expanded2poly (t:term) v =
1512 (*if is_expanded t then*) term2poly' t v
1513 (*else raise error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*);
1515 (*. same as expanded2poly with automatic detection of the variables .*)
1516 fun expanded2polyx t = expanded2poly t (((map free2str) o vars) t);
1518 (*. converts a powerproduct into term representation .*)
1519 fun powerproduct2term(xs,v) =
1531 if list_is_null(tl(!xss)) then
1533 if hd(!xss)=1 then Free(hd(!vv), HOLogic.realT)
1536 Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1537 Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT)
1544 Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1545 Free(hd(!vv), HOLogic.realT) $
1546 powerproduct2term(tl(!xss),tl(!vv))
1550 Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1552 Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1553 Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT)
1555 powerproduct2term(tl(!xss),tl(!vv))
1561 (*. converts a monom into term representation .*)
1562 (*fun monom2term ((c,e):mv_monom, v:string list) =
1563 if c=0 then Free(str_of_int 0,HOLogic.realT)
1566 if list_is_null(e) then
1568 Free(str_of_int c,HOLogic.realT)
1574 powerproduct2term(e,v)
1578 Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1579 Free(str_of_int c,HOLogic.realT) $
1580 powerproduct2term(e,v)
1586 (*fun monom2term ((i, is):mv_monom, v) =
1590 then Free (str_of_int i, HOLogic.realT)
1591 else Const ("uminus", HOLogic.realT --> HOLogic.realT) $
1592 Free ((str_of_int o abs) i, HOLogic.realT)
1595 then Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
1596 (Free (str_of_int i, HOLogic.realT)) $
1597 powerproduct2term(is, v)
1598 else Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
1599 (Const ("uminus", HOLogic.realT --> HOLogic.realT) $
1600 Free ((str_of_int o abs) i, HOLogic.realT)) $
1601 powerproduct2term(is, vs);---------------------------*)
1602 fun monom2term ((i, is) : mv_monom, vs) =
1604 then Free (str_of_int i, HOLogic.realT)
1606 then powerproduct2term (is, vs)
1607 else Const ("op *", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
1608 (Free (str_of_int i, HOLogic.realT)) $
1609 powerproduct2term (is, vs);
1611 (*. converts the internal polynomial representation into an Isabelle term.*)
1612 fun poly2term' ([] : mv_poly, vs) = Free(str_of_int 0, HOLogic.realT)
1613 | poly2term' ([(c, e) : mv_monom], vs) = monom2term ((c, e), vs)
1614 | poly2term' ((c, e) :: ces, vs) =
1615 Const("op +", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
1616 poly2term (ces, vs) $ monom2term ((c, e), vs)
1617 and poly2term (xs, vs) = poly2term' (rev (sort (mv_geq LEX_) (xs)), vs);
1620 (*. converts a monom into term representation .*)
1621 (*. ignores the sign of the coefficients => use only for exp-poly functions .*)
1622 fun monom2term2((c,e):mv_monom, v:string list) =
1623 if c=0 then Free(str_of_int 0,HOLogic.realT)
1626 if list_is_null(e) then
1628 Free(str_of_int (abs(c)),HOLogic.realT)
1634 powerproduct2term(e,v)
1638 Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1639 Free(str_of_int (abs(c)),HOLogic.realT) $
1640 powerproduct2term(e,v)
1645 (*. converts the expanded polynomial representation into the term representation .*)
1646 fun exp2term' ([]:mv_poly,vars) = Free(str_of_int 0,HOLogic.realT)
1647 | exp2term' ([(c,e)],vars) = monom2term((c,e),vars)
1648 | exp2term' ((c1,e1)::others,vars) =
1650 Const("op -",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1651 exp2term'(others,vars) $
1653 monom2term2((c1,e1),vars)
1656 Const("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1657 exp2term'(others,vars) $
1659 monom2term2((c1,e1),vars)
1662 (*. sorts the powerproduct by lexicographic termorder and converts them into
1663 a term in polynomial representation .*)
1664 fun poly2expanded (xs,vars) = exp2term'(rev(sort (mv_geq LEX_) (xs)),vars);
1666 (*. converts a polynomial into expanded form .*)
1667 fun polynomial2expanded t =
1669 val vars=(((map free2str) o vars) t);
1671 SOME (poly2expanded (the (term2poly t vars), vars))
1672 end) handle _ => NONE;
1674 (*. converts a polynomial into polynomial form .*)
1675 fun expanded2polynomial t =
1677 val vars=(((map free2str) o vars) t);
1679 SOME (poly2term (the (expanded2poly t vars), vars))
1680 end) handle _ => NONE;
1683 (*. calculates the greatest common divisor of numerator and denominator and seperates it from each .*)
1684 fun step_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) =
1689 val vars = rev(get_vars(p1) union get_vars(p2));
1692 p1':= sort (mv_geq LEX_) (the (term2poly p1 vars ));
1693 p2':= sort (mv_geq LEX_) (the (term2poly p2 vars ));
1694 p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1695 if (!p3)=[(1,mv_null2(vars))] then
1697 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2
1702 p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1703 p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1705 if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1707 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1710 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1711 poly2term(!p1',vars) $
1716 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1717 poly2term(!p2',vars) $
1723 p1':=mv_skalar_mul(!p1',~1);
1724 p2':=mv_skalar_mul(!p2',~1);
1725 p3:=mv_skalar_mul(!p3,~1);
1727 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1730 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1731 poly2term(!p1',vars) $
1736 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1737 poly2term(!p2',vars) $
1745 | step_cancel _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction");
1748 (*. same as step_cancel, this time for expanded forms (input+output) .*)
1749 fun step_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) =
1754 val vars = rev(get_vars(p1) union get_vars(p2));
1757 p1':= sort (mv_geq LEX_) (the (expanded2poly p1 vars ));
1758 p2':= sort (mv_geq LEX_) (the (expanded2poly p2 vars ));
1759 p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1760 if (!p3)=[(1,mv_null2(vars))] then
1762 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2
1767 p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1768 p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1770 if #1(hd(sort (mv_geq LEX_) (!p2')))(* mv_lc2(!p2',LEX_)*)>0 then
1772 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1775 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1776 poly2expanded(!p1',vars) $
1777 poly2expanded(!p3,vars)
1781 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1782 poly2expanded(!p2',vars) $
1783 poly2expanded(!p3,vars)
1788 p1':=mv_skalar_mul(!p1',~1);
1789 p2':=mv_skalar_mul(!p2',~1);
1790 p3:=mv_skalar_mul(!p3,~1);
1792 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1795 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1796 poly2expanded(!p1',vars) $
1797 poly2expanded(!p3,vars)
1801 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
1802 poly2expanded(!p2',vars) $
1803 poly2expanded(!p3,vars)
1810 | step_cancel_expanded _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction");
1812 (*. calculates the greatest common divisor of numerator and denominator and divides each through it .*)
1813 fun direct_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) =
1818 val vars = rev(get_vars(p1) union get_vars(p2));
1821 p1':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p1 vars )),LEX_));
1822 p2':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p2 vars )),LEX_));
1823 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1825 if (!p3)=[(1,mv_null2(vars))] then
1827 (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
1831 p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1832 p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1833 if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1836 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1839 poly2term((!p1'),vars)
1843 poly2term((!p2'),vars)
1847 if mv_grad(!p3)>0 then
1850 Const ("Not",[bool]--->bool) $
1852 Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
1853 poly2term((!p3),vars) $
1854 Free("0",HOLogic.realT)
1863 p1':=mv_skalar_mul(!p1',~1);
1864 p2':=mv_skalar_mul(!p2',~1);
1865 if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);
1867 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1870 poly2term((!p1'),vars)
1874 poly2term((!p2'),vars)
1877 if mv_grad(!p3)>0 then
1880 Const ("Not",[bool]--->bool) $
1882 Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
1883 poly2term((!p3),vars) $
1884 Free("0",HOLogic.realT)
1895 | direct_cancel _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");
1897 (*. same es direct_cancel, this time for expanded forms (input+output).*)
1898 fun direct_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) =
1903 val vars = rev(get_vars(p1) union get_vars(p2));
1906 p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_));
1907 p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_));
1908 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
1910 if (!p3)=[(1,mv_null2(vars))] then
1912 (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
1916 p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
1917 p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
1918 if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
1921 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1924 poly2expanded((!p1'),vars)
1928 poly2expanded((!p2'),vars)
1932 if mv_grad(!p3)>0 then
1935 Const ("Not",[bool]--->bool) $
1937 Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
1938 poly2expanded((!p3),vars) $
1939 Free("0",HOLogic.realT)
1948 p1':=mv_skalar_mul(!p1',~1);
1949 p2':=mv_skalar_mul(!p2',~1);
1950 if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);
1952 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
1955 poly2expanded((!p1'),vars)
1959 poly2expanded((!p2'),vars)
1962 if mv_grad(!p3)>0 then
1965 Const ("Not",[bool]--->bool) $
1967 Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
1968 poly2expanded((!p3),vars) $
1969 Free("0",HOLogic.realT)
1980 | direct_cancel_expanded _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");
1983 (*. adds two fractions .*)
1984 fun add_fract ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =
1986 val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
1987 val t11'=ref (the(term2poly t11 vars));
1988 val _= writeln"### add_fract: done t11"
1989 val t12'=ref (the(term2poly t12 vars));
1990 val _= writeln"### add_fract: done t12"
1991 val t21'=ref (the(term2poly t21 vars));
1992 val _= writeln"### add_fract: done t21"
1993 val t22'=ref (the(term2poly t22 vars));
1994 val _= writeln"### add_fract: done t22"
2002 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
2003 writeln"### add_fract: done sort mv_lcm";
2004 m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
2005 writeln"### add_fract: done sort mv_division t12";
2006 m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
2007 writeln"### add_fract: done sort mv_division t22";
2008 nom :=sort (mv_geq LEX_)
2009 (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),
2010 mv_mul(!t21',!m2,LEX_),
2013 writeln"### add_fract: done sort mv_add";
2015 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2018 poly2term((!nom),vars)
2022 poly2term((!den),vars)
2027 | add_fract (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call");
2029 (*. adds two expanded fractions .*)
2030 fun add_fract_exp ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =
2032 val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
2033 val t11'=ref (the(expanded2poly t11 vars));
2034 val t12'=ref (the(expanded2poly t12 vars));
2035 val t21'=ref (the(expanded2poly t21 vars));
2036 val t22'=ref (the(expanded2poly t22 vars));
2044 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
2045 m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
2046 m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
2047 nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_));
2049 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2052 poly2expanded((!nom),vars)
2056 poly2expanded((!den),vars)
2061 | add_fract_exp (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call");
2063 (*. adds a list of terms .*)
2064 fun add_list_of_fractions []= (Free("0",HOLogic.realT),[])
2065 | add_list_of_fractions [x]= direct_cancel x
2066 | add_list_of_fractions (x::y::xs) =
2068 val (t1a,rest1)=direct_cancel(x);
2069 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(x)";
2070 val (t2a,rest2)=direct_cancel(y);
2071 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(y)";
2072 val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs));
2073 val _= writeln"### add_list_of_fractions xs: has done add_list_of_fraction xs";
2074 val (t4a,rest4)=direct_cancel(t3a);
2075 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(t3a)";
2076 val rest=rest1 union rest2 union rest3 union rest4;
2078 (writeln"### add_list_of_fractions in";
2085 (*. adds a list of expanded terms .*)
2086 fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])
2087 | add_list_of_fractions_exp [x]= direct_cancel_expanded x
2088 | add_list_of_fractions_exp (x::y::xs) =
2090 val (t1a,rest1)=direct_cancel_expanded(x);
2091 val (t2a,rest2)=direct_cancel_expanded(y);
2092 val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));
2093 val (t4a,rest4)=direct_cancel_expanded(t3a);
2094 val rest=rest1 union rest2 union rest3 union rest4;
2101 (*. calculates the lcm of a list of mv_poly .*)
2102 fun calc_lcm ([x],var)= (x,var)
2103 | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
2105 (*. converts a list of terms to a list of mv_poly .*)
2107 | t2d((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);
2109 (*. same as t2d, this time for expanded forms .*)
2110 fun t2d_exp([],_)=[]
2111 | t2d_exp((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
2113 (*. converts a list of fract terms to a list of their denominators .*)
2114 fun termlist2denominators [] = ([],[])
2115 | termlist2denominators xs =
2121 while length(!xxs)>0 do
2124 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
2128 var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
2135 (*. calculates the lcm of a list of mv_poly .*)
2136 fun calc_lcm ([x],var)= (x,var)
2137 | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
2139 (*. converts a list of terms to a list of mv_poly .*)
2141 | t2d((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);
2143 (*. same as t2d, this time for expanded forms .*)
2144 fun t2d_exp([],_)=[]
2145 | t2d_exp((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
2147 (*. converts a list of fract terms to a list of their denominators .*)
2148 fun termlist2denominators [] = ([],[])
2149 | termlist2denominators xs =
2155 while length(!xxs)>0 do
2158 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
2162 var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
2169 (*. same as termlist2denminators, this time for expanded forms .*)
2170 fun termlist2denominators_exp [] = ([],[])
2171 | termlist2denominators_exp xs =
2177 while length(!xxs)>0 do
2180 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
2184 var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
2188 (t2d_exp(xs,!var),!var)
2191 (*. reduces all fractions to the least common denominator .*)
2192 fun com_den(x::xs,denom,den,var)=
2194 val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
2195 val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
2196 val p3= #1(mv_division(denom,p2,LEX_));
2197 val p1var=get_vars(p1');
2199 if length(xs)>0 then
2200 if p3=[(1,mv_null2(var))] then
2202 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2205 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2207 poly2term(the (term2poly p1' p1var),p1var)
2212 #1(com_den(xs,denom,den,var))
2218 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2221 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2224 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2225 poly2term(the (term2poly p1' p1var),p1var) $
2234 #1(com_den(xs,denom,den,var))
2239 if p3=[(1,mv_null2(var))] then
2242 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2244 poly2term(the (term2poly p1' p1var),p1var)
2253 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2256 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2257 poly2term(the (term2poly p1' p1var),p1var) $
2267 (*. same as com_den, this time for expanded forms .*)
2268 fun com_den_exp(x::xs,denom,den,var)=
2270 val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
2271 val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
2272 val p3= #1(mv_division(denom,p2,LEX_));
2273 val p1var=get_vars(p1');
2275 if length(xs)>0 then
2276 if p3=[(1,mv_null2(var))] then
2278 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2281 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2283 poly2expanded(the(expanded2poly p1' p1var),p1var)
2288 #1(com_den_exp(xs,denom,den,var))
2294 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2297 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2300 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2301 poly2expanded(the(expanded2poly p1' p1var),p1var) $
2302 poly2expanded(p3,var)
2310 #1(com_den_exp(xs,denom,den,var))
2315 if p3=[(1,mv_null2(var))] then
2318 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2320 poly2expanded(the(expanded2poly p1' p1var),p1var)
2329 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
2332 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2333 poly2expanded(the(expanded2poly p1' p1var),p1var) $
2334 poly2expanded(p3,var)
2343 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon
2344 -------------------------------------------------------------
2345 (* WN0210???SK brauch ma des überhaupt *)
2346 fun com_den2(x::xs,denom,den,var)=
2348 val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
2349 val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
2350 val p3= #1(mv_division(denom,p2,LEX_));
2351 val p1var=get_vars(p1');
2353 if length(xs)>0 then
2354 if p3=[(1,mv_null2(var))] then
2356 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2357 poly2term(the(term2poly p1' p1var),p1var) $
2358 com_den2(xs,denom,den,var)
2362 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2365 val p3'=poly2term(p3,var);
2366 val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
2368 poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
2371 com_den2(xs,denom,den,var)
2374 if p3=[(1,mv_null2(var))] then
2376 poly2term(the(term2poly p1' p1var),p1var)
2381 val p3'=poly2term(p3,var);
2382 val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
2384 poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
2389 (* WN0210???SK brauch ma des überhaupt *)
2390 fun com_den_exp2(x::xs,denom,den,var)=
2392 val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
2393 val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
2394 val p3= #1(mv_division(denom,p2,LEX_));
2395 val p1var=get_vars p1';
2397 if length(xs)>0 then
2398 if p3=[(1,mv_null2(var))] then
2400 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2401 poly2expanded(the (expanded2poly p1' p1var),p1var) $
2402 com_den_exp2(xs,denom,den,var)
2406 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2409 val p3'=poly2expanded(p3,var);
2410 val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
2412 poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
2415 com_den_exp2(xs,denom,den,var)
2418 if p3=[(1,mv_null2(var))] then
2420 poly2expanded(the (expanded2poly p1' p1var),p1var)
2425 val p3'=poly2expanded(p3,var);
2426 val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
2428 poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
2432 ---------------------------------------------------------*)
2435 (*. searches for an element y of a list ys, which has an gcd not 1 with x .*)
2436 fun exists_gcd (x,[]) = false
2437 | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then exists_gcd (x,ys)
2440 (*. divides each element of the list xs with y .*)
2441 fun list_div ([],y) = []
2442 | list_div (x::xs,y) =
2444 val (d,r)=mv_division(x,y,LEX_);
2448 else x::list_div(xs,y)
2451 (*. checks if x is in the list ys .*)
2452 fun in_list (x,[]) = false
2453 | in_list (x,y::ys) = if x=y then true
2456 (*. deletes all equal elements of the list xs .*)
2457 fun kill_equal [] = []
2458 | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)
2459 else x::kill_equal(xs);
2461 (*. searches for new factors .*)
2462 fun new_factors [] = []
2463 | new_factors (list:mv_poly list):mv_poly list =
2465 val l = kill_equal list;
2466 val len = length(l);
2474 if gcd=[(1,mv_null2(#2(hd(x))))] then
2476 if exists_gcd(x,xs) then new_factors (y::xs @ [x])
2477 else x::new_factors(y::xs)
2479 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))
2483 if len=1 then [hd(l)]
2487 (*. gets the factors of a list .*)
2488 fun get_factors x = new_factors x;
2490 (*. multiplies the elements of the list .*)
2491 fun multi_list [] = []
2492 | multi_list (x::xs) = if xs=[] then x
2493 else mv_mul(x,multi_list xs,LEX_);
2495 (*. makes a term out of the elements of the list (polynomial representation) .*)
2496 fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)
2497 | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)
2500 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2501 poly2term(sort (mv_geq LEX_) (x),vars) $
2505 (*. factorizes the denominator (polynomial representation) .*)
2506 fun factorize_den (l,den,vars) =
2508 val factor_list=kill_equal( (get_factors l));
2509 val mlist=multi_list(factor_list);
2510 val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
2514 if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)
2515 else make_term(last::factor_list,vars)
2517 else raise error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")
2520 (*. makes a term out of the elements of the list (expanded polynomial representation) .*)
2521 fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)
2522 | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)
2525 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2526 poly2expanded(sort (mv_geq LEX_) (x),vars) $
2530 (*. factorizes the denominator (expanded polynomial representation) .*)
2531 fun factorize_den_exp (l,den,vars) =
2533 val factor_list=kill_equal( (get_factors l));
2534 val mlist=multi_list(factor_list);
2535 val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
2539 if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)
2540 else make_exp(last::factor_list,vars)
2542 else raise error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")
2545 (*. calculates the common denominator of all elements of the list and multiplies .*)
2546 (*. the nominators and denominators with the correct factor .*)
2547 (*. (polynomial representation) .*)
2548 fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)
2549 | step_add_list_of_fractions [x]= raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")
2550 | step_add_list_of_fractions (xs) =
2552 val den_list=termlist2denominators (xs); (* list of denominators *)
2553 val (denom,var)=calc_lcm(den_list); (* common denominator *)
2554 val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
2556 com_den(xs,denom,den,var)
2559 (*. calculates the common denominator of all elements of the list and multiplies .*)
2560 (*. the nominators and denominators with the correct factor .*)
2561 (*. (expanded polynomial representation) .*)
2562 fun step_add_list_of_fractions_exp [] = (Free("0",HOLogic.realT),[]:term list)
2563 | step_add_list_of_fractions_exp [x] = raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")
2564 | step_add_list_of_fractions_exp (xs)=
2566 val den_list=termlist2denominators_exp (xs); (* list of denominators *)
2567 val (denom,var)=calc_lcm(den_list); (* common denominator *)
2568 val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
2570 com_den_exp(xs,denom,den,var)
2573 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon
2574 -------------------------------------------------------------
2575 (* WN0210???SK brauch ma des überhaupt *)
2576 fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)
2577 | step_add_list_of_fractions2 [x]=(x,[])
2578 | step_add_list_of_fractions2 (xs) =
2580 val den_list=termlist2denominators (xs); (* list of denominators *)
2581 val (denom,var)=calc_lcm(den_list); (* common denominator *)
2582 val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
2585 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2586 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $
2587 poly2term(denom,var)
2593 (* WN0210???SK brauch ma des überhaupt *)
2594 fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)
2595 | step_add_list_of_fractions2_exp [x]=(x,[])
2596 | step_add_list_of_fractions2_exp (xs) =
2598 val den_list=termlist2denominators_exp (xs); (* list of denominators *)
2599 val (denom,var)=calc_lcm(den_list); (* common denominator *)
2600 val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
2603 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2604 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $
2605 poly2expanded(denom,var)
2610 ---------------------------------------------- *)
2613 (*. converts a term, which contains severel terms seperated by +, into a list of these terms .*)
2614 fun term2list (t as (Const("HOL.divide",_) $ _ $ _)) = [t]
2615 | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =
2616 [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2617 t $ Free("1",HOLogic.realT)
2619 | term2list (t as (Free(_,_))) =
2620 [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2621 t $ Free("1",HOLogic.realT)
2623 | term2list (t as (Const("op *",_) $ _ $ _)) =
2624 [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
2625 t $ Free("1",HOLogic.realT)
2627 | term2list (Const("op +",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)
2628 | term2list (Const("op -",_) $ t1 $ t2) =
2629 raise error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")
2630 | term2list _ = raise error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");
2632 (*.factors out the gcd of nominator and denominator:
2633 a/b = (a' * gcd)/(b' * gcd), a,b,gcd are poly[2].*)
2634 fun factout_p_ (thy:theory) t = SOME (step_cancel t,[]:term list);
2635 fun factout_ (thy:theory) t = SOME (step_cancel_expanded t,[]:term list);
2637 (*.cancels a single fraction with normalform [2]
2638 resulting in a canceled fraction [2], see factout_ .*)
2639 fun cancel_p_ (thy:theory) t = (*WN.2.6.03 no rewrite -> NONE !*)
2640 (let val (t',asm) = direct_cancel(*_expanded ... corrected MG.21.8.03*) t
2641 in if t = t' then NONE else SOME (t',asm)
2642 end) handle _ => NONE;
2643 (*.the same as above with normalform [3]
2645 theory -> (*10.02 unused *)
2646 term -> (*fraction in normalform [3] *)
2647 (term * (*fraction in normalform [3] *)
2648 term list) (*casual asumptions in normalform [3] *)
2649 option (*NONE: the function is not applicable *).*)
2650 fun cancel_ (thy:theory) t = SOME (direct_cancel_expanded t) handle _ => NONE;
2652 (*.transforms sums of at least 2 fractions [3] to
2653 sums with the least common multiple as nominator.*)
2654 fun common_nominator_p_ (thy:theory) t =
2655 ((*writeln("### common_nominator_p_ called");*)
2656 SOME (step_add_list_of_fractions(term2list(t))) handle _ => NONE
2658 fun common_nominator_ (thy:theory) t =
2659 SOME (step_add_list_of_fractions_exp(term2list(t))) handle _ => NONE;
2661 (*.add 2 or more fractions
2662 val add_fraction_p_ :
2663 theory -> (*10.02 unused *)
2664 term -> (*2 or more fractions with normalform [2] *)
2665 (term * (*one fraction with normalform [2] *)
2666 term list) (*casual assumptions in normalform [2] WN0210???SK *)
2667 option (*NONE: the function is not applicable *).*)
2668 fun add_fraction_p_ (thy:theory) t =
2669 (writeln("### add_fraction_p_ called");
2670 (let val ts = term2list t
2672 then SOME (add_list_of_fractions ts)
2673 else NONE (*raise error ("RATIONALS_ADD_EXCEPTION: nothing to add")*)
2674 end) handle _ => NONE
2676 (*.same as add_fraction_p_ but with normalform [3].*)
2677 (*SOME (step_add_list_of_fractions2(term2list(t))); *)
2678 fun add_fraction_ (thy:theory) t =
2679 if length(term2list(t))>1
2680 then SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE
2681 else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)
2683 fun add_fraction_ (thy:theory) t =
2684 (if 1 < length (term2list t)
2685 then SOME (add_list_of_fractions_exp (term2list t))
2686 else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)
2687 NONE) handle _ => NONE;
2689 (*SOME (step_add_list_of_fractions2_exp(term2list(t))); *)
2691 (*. brings the term into a normal form .*)
2692 fun norm_rational_ (thy:theory) t =
2693 SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;
2694 fun norm_expanded_rat_ (thy:theory) t =
2695 SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;
2698 (*.evaluates conditions in calculate_Rational.*)
2699 (*make local with FIXX@ME result:term *term list*)
2700 val calc_rat_erls = prep_rls(
2701 Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
2702 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [], *)
2704 [Calc ("op =",eval_equal "#equal_"),
2705 Calc ("Atools.is'_const",eval_const "#is_const_"),
2706 Thm ("not_true",num_str not_true),
2707 Thm ("not_false",num_str not_false)
2712 (*.simplifies expressions with numerals;
2713 does NOT rearrange the term by AC-rewriting; thus terms with variables
2714 need to have constants to be commuted together respectively.*)
2715 val calculate_Rational = prep_rls(
2716 merge_rls "calculate_Rational"
2717 (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
2718 erls = calc_rat_erls, srls = Erls, (*asm_thm = [],*)
2721 [Calc ("HOL.divide" ,eval_cancel "#divide_"),
2723 Thm ("sym_real_minus_divide_eq",
2724 num_str (real_minus_divide_eq RS sym)),
2725 (*SYM - ?x / ?y = - (?x / ?y) may come from subst*)
2727 Thm ("rat_add",num_str rat_add),
2728 (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \
2729 \"a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)
2730 Thm ("rat_add1",num_str rat_add1),
2731 (*"[| a is_const; b is_const; c is_const |] ==> \
2732 \"a / c + b / c = (a + b) / c"*)
2733 Thm ("rat_add2",num_str rat_add2),
2734 (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> \
2735 \?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)
2736 Thm ("rat_add3",num_str rat_add3),
2737 (*"[| a is_const; b is_const; c is_const |] ==> \
2738 \"a + b / c = (a * c) / c + b / c"\
2739 \.... is_const to be omitted here FIXME*)
2741 Thm ("rat_mult",num_str rat_mult),
2742 (*a / b * (c / d) = a * c / (b * d)*)
2743 Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),
2744 (*?x * (?y / ?z) = ?x * ?y / ?z*)
2745 Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),
2746 (*?y / ?z * ?x = ?y * ?x / ?z*)
2748 Thm ("real_divide_divide1",num_str real_divide_divide1),
2749 (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)
2750 Thm ("real_divide_divide2_eq",num_str real_divide_divide2_eq),
2751 (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
2753 Thm ("rat_power", num_str rat_power),
2754 (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
2756 Thm ("mult_cross",num_str mult_cross),
2757 (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)
2758 Thm ("mult_cross1",num_str mult_cross1),
2759 (*" b ~= 0 ==> (a / b = c ) = (a = b * c)*)
2760 Thm ("mult_cross2",num_str mult_cross2)
2761 (*" d ~= 0 ==> (a = c / d) = (a * d = c)*)
2766 (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)
2767 fun eval_is_expanded (thmid:string) _
2768 (t as (Const("Rational.is'_expanded", _) $ arg)) thy =
2770 then SOME (mk_thmid thmid ""
2771 ((Syntax.string_of_term (thy2ctxt thy)) arg) "",
2772 Trueprop $ (mk_equality (t, HOLogic.true_const)))
2773 else SOME (mk_thmid thmid ""
2774 ((Syntax.string_of_term (thy2ctxt thy)) arg) "",
2775 Trueprop $ (mk_equality (t, HOLogic.false_const)))
2776 | eval_is_expanded _ _ _ _ = NONE;
2779 merge_rls "rational_erls" calculate_Rational
2780 (append_rls "is_expanded" Atools_erls
2781 [Calc ("Rational.is'_expanded", eval_is_expanded "")
2786 (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:
2787 =================================================================
2790 B[2] 'common_nominator_p': transforms summands in a term [2]
2791 to fractions with the (least) common multiple as nominator.
2792 B[3] 'norm_rational': normalizes arbitrary algebraic terms (without
2793 radicals and transzendental functions) to one canceled fraction,
2794 nominator and denominator in polynomial form.
2796 In order to meet isac's requirements for interactive and stepwise calculation,
2797 each 'reverse-rewerite-set' consists of an initialization for the interpreter
2798 state and of 4 functions, each of which employs rewriting as much as possible.
2799 The signature of these functions are the same in each 'reverse-rewrite-set'
2802 (* ************************************************************************* *)
2806 ------------------------
2807 cancels a single fraction consisting of two (uni- or multivariate)
2808 polynomials WN0609???SK[2] into another such a fraction; examples:
2811 -------------------- = ---------
2812 a^2 + -2*a*b + b^2 a + -1*b
2818 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)
2819 (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)
2821 val {rules, rew_ord=(_,ro),...} =
2822 rep_rls (assoc_rls "make_polynomial");
2823 (*WN060829 ... make_deriv does not terminate with 1st expl above,
2824 see rational.sml --- investigate rulesets for cancel_p ---*)
2825 val {rules, rew_ord=(_,ro),...} =
2826 rep_rls (assoc_rls "rev_rew_p");
2828 val thy = Rational.thy;
2830 (*.init_state = fn : term -> istate
2831 initialzies the state of the script interpreter. The state is:
2833 type rrlsstate = (*state for reverse rewriting*)
2834 (term * (*the current formula*)
2835 term * (*the final term*)
2836 rule list (*'reverse rule list' (#)*)
2837 list * (*may be serveral, eg. in norm_rational*)
2838 (rule * (*Thm (+ Thm generated from Calc) resulting in ...*)
2839 (term * (*... rewrite with ...*)
2840 term list)) (*... assumptions*)
2841 list); (*derivation from given term to normalform
2842 in reverse order with sym_thm;
2843 (#) could be extracted from here by (map #1)*).*)
2844 (* val {rules, rew_ord=(_,ro),...} =
2845 rep_rls (assoc_rls "rev_rew_p") (*USE ALWAYS, SEE val cancel_p*);
2846 val (thy, eval_rls, ro) =(Rational.thy, Atools_erls, ro) (*..val cancel_p*);
2849 fun init_state thy eval_rls ro t =
2850 let val SOME (t',_) = factout_p_ thy t
2851 val SOME (t'',asm) = cancel_p_ thy t
2852 val der = reverse_deriv thy eval_rls rules ro NONE t'
2853 val der = der @ [(Thm ("real_mult_div_cancel2",
2854 num_str real_mult_div_cancel2),
2856 val rs = (distinct_Thm o (map #1)) der
2857 val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
2860 (*..insufficient,eg.make_Polynomial*)])rs
2861 in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end;
2863 (*.locate_rule = fn : rule list -> term -> rule
2864 -> (rule * (term * term list) option) list.
2865 checks a rule R for being a cancel-rule, and if it is,
2866 then return the list of rules (+ the terms they are rewriting to)
2867 which need to be applied before R should be applied.
2868 precondition: the rule is applicable to the argument-term.
2870 rule list: the reverse rule list
2871 -> term : ... to which the rule shall be applied
2872 -> rule : ... to be applied to term
2874 -> (rule : a rule rewriting to ...
2875 * (term : ... the resulting term ...
2876 * term list): ... with the assumptions ( //#0).
2877 ) list : there may be several such rules;
2878 the list is empty, if the rule has nothing to do
2882 fun locate_rule thy eval_rls ro [rs] t r =
2883 if (id_of_thm r) mem (map (id_of_thm)) rs
2885 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
2887 SOME ta => [(r, ta)]
2888 | NONE => (writeln("### locate_rule: rewrite "^
2889 (id_of_thm r)^" "^(term2str t)^" = NONE");
2891 else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
2892 | locate_rule _ _ _ _ _ _ =
2893 raise error ("locate_rule: doesnt match rev-sets in istate");
2895 (*.next_rule = fn : rule list -> term -> rule option
2896 for a given term return the next rules to be done for cancelling.
2898 rule list : the reverse rule list
2899 term : the term for which ...
2901 -> rule option: ... this rule is appropriate for cancellation;
2902 there may be no such rule (if the term is canceled already.*)
2903 (* val thy = Rational.thy;
2904 val Rrls {rew_ord=(_,ro),...} = cancel;
2905 val ([rs],t) = (rss,f);
2906 next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
2908 val (thy, [rs]) = (Rational.thy, revsets);
2909 val Rrls {rew_ord=(_,ro),...} = cancel;
2912 fun next_rule thy eval_rls ro [rs] t =
2913 let val der = make_deriv thy eval_rls rs ro NONE t;
2915 (* val (_,r,_)::_ = der;
2917 (_,r,_)::_ => SOME r
2920 | next_rule _ _ _ _ _ =
2921 raise error ("next_rule: doesnt match rev-sets in istate");
2923 (*.val attach_form = f : rule list -> term -> term
2924 -> (rule * (term * term list)) list
2925 checks an input term TI, if it may belong to a current cancellation, by
2926 trying to derive it from the given term TG.
2928 term : TG, the last one in the cancellation agreed upon by user + math-eng
2929 -> term: TI, the next one input by the user
2931 -> (rule : the rule to be applied in order to reach TI
2932 * (term : ... obtained by applying the rule ...
2933 * term list): ... and the respective assumptions.
2934 ) list : there may be several such rules;
2935 the list is empty, if the users term does not belong
2936 to a cancellation of the term last agreed upon.*)
2937 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
2938 []:(rule * (term * term list)) list;
2943 Rrls {id = "cancel_p", prepat=[],
2944 rew_ord=("ord_make_polynomial",
2945 ord_make_polynomial false Rational.thy),
2946 erls = rational_erls,
2947 calc = [("PLUS" ,("op +" ,eval_binop "#add_")),
2948 ("TIMES" ,("op *" ,eval_binop "#mult_")),
2949 ("DIVIDE" ,("HOL.divide" ,eval_cancel "#divide_")),
2950 ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
2951 (*asm_thm=[("real_mult_div_cancel2","")],*)
2952 scr=Rfuns {init_state = init_state thy Atools_erls ro,
2953 normal_form = cancel_p_ thy,
2954 locate_rule = locate_rule thy Atools_erls ro,
2955 next_rule = next_rule thy Atools_erls ro,
2956 attach_form = attach_form}}
2960 local(*.ad (1) 'cancel'
2961 ------------------------------
2962 cancels a single fraction consisting of two (uni- or multivariate)
2963 polynomials WN0609???SK[3] into another such a fraction; examples:
2966 -------------------- = ---------
2967 a^2 - 2*a*b + b^2 a - *b
2969 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)
2970 (*WN 24.8.02: wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)
2973 val SOME (Rls {rules=rules,rew_ord=(_,ro),...}) =
2974 assoc'(!ruleset',"expand_binoms");
2976 val {rules=rules,rew_ord=(_,ro),...} =
2977 rep_rls (assoc_rls "expand_binoms");
2978 val thy = Rational.thy;
2980 fun init_state thy eval_rls ro t =
2981 let val SOME (t',_) = factout_ thy t;
2982 val SOME (t'',asm) = cancel_ thy t;
2983 val der = reverse_deriv thy eval_rls rules ro NONE t';
2984 val der = der @ [(Thm ("real_mult_div_cancel2",
2985 num_str real_mult_div_cancel2),
2987 val rs = map #1 der;
2988 in (t,t'',[rs],der) end;
2990 fun locate_rule thy eval_rls ro [rs] t r =
2991 if (id_of_thm r) mem (map (id_of_thm)) rs
2993 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
2995 SOME ta => [(r, ta)]
2996 | NONE => (writeln("### locate_rule: rewrite "^
2997 (id_of_thm r)^" "^(term2str t)^" = NONE");
2999 else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
3000 | locate_rule _ _ _ _ _ _ =
3001 raise error ("locate_rule: doesnt match rev-sets in istate");
3003 fun next_rule thy eval_rls ro [rs] t =
3004 let val der = make_deriv thy eval_rls rs ro NONE t;
3006 (* val (_,r,_)::_ = der;
3008 (_,r,_)::_ => SOME r
3011 | next_rule _ _ _ _ _ =
3012 raise error ("next_rule: doesnt match rev-sets in istate");
3014 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
3015 []:(rule * (term * term list)) list;
3017 val pat = (term_of o the o (parse thy)) "?r/?s";
3018 val pre1 = (term_of o the o (parse thy)) "?r is_expanded";
3019 val pre2 = (term_of o the o (parse thy)) "?s is_expanded";
3020 val prepat = [([pre1, pre2], pat)];
3026 Rrls {id = "cancel", prepat=prepat,
3027 rew_ord=("ord_make_polynomial",
3028 ord_make_polynomial false Rational.thy),
3029 erls = rational_erls,
3030 calc = [("PLUS" ,("op +" ,eval_binop "#add_")),
3031 ("TIMES" ,("op *" ,eval_binop "#mult_")),
3032 ("DIVIDE" ,("HOL.divide" ,eval_cancel "#divide_")),
3033 ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
3034 scr=Rfuns {init_state = init_state thy Atools_erls ro,
3035 normal_form = cancel_ thy,
3036 locate_rule = locate_rule thy Atools_erls ro,
3037 next_rule = next_rule thy Atools_erls ro,
3038 attach_form = attach_form}}
3043 local(*.ad [2] 'common_nominator_p'
3044 ---------------------------------
3045 FIXME Beschreibung .*)
3048 val {rules=rules,rew_ord=(_,ro),...} =
3049 rep_rls (assoc_rls "make_polynomial");
3050 (*WN060829 ... make_deriv does not terminate with 1st expl above,
3051 see rational.sml --- investigate rulesets for cancel_p ---*)
3052 val {rules, rew_ord=(_,ro),...} =
3053 rep_rls (assoc_rls "rev_rew_p");
3054 val thy = Rational.thy;
3057 (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option
3060 (*.init_state = fn : term -> istate
3061 initialzies the state of the interactive interpreter. The state is:
3063 type rrlsstate = (*state for reverse rewriting*)
3064 (term * (*the current formula*)
3065 term * (*the final term*)
3066 rule list (*'reverse rule list' (#)*)
3067 list * (*may be serveral, eg. in norm_rational*)
3068 (rule * (*Thm (+ Thm generated from Calc) resulting in ...*)
3069 (term * (*... rewrite with ...*)
3070 term list)) (*... assumptions*)
3071 list); (*derivation from given term to normalform
3072 in reverse order with sym_thm;
3073 (#) could be extracted from here by (map #1)*).*)
3074 fun init_state thy eval_rls ro t =
3075 let val SOME (t',_) = common_nominator_p_ thy t;
3076 val SOME (t'',asm) = add_fraction_p_ thy t;
3077 val der = reverse_deriv thy eval_rls rules ro NONE t';
3078 val der = der @ [(Thm ("real_mult_div_cancel2",
3079 num_str real_mult_div_cancel2),
3081 val rs = (distinct_Thm o (map #1)) der;
3082 val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
3084 "sym_real_mult_1"]) rs;
3085 in (t,t'',[rs(*here only _ONE_*)],der) end;
3087 (* use"knowledge/Rational.ML";
3090 (*.locate_rule = fn : rule list -> term -> rule
3091 -> (rule * (term * term list) option) list.
3092 checks a rule R for being a cancel-rule, and if it is,
3093 then return the list of rules (+ the terms they are rewriting to)
3094 which need to be applied before R should be applied.
3095 precondition: the rule is applicable to the argument-term.
3097 rule list: the reverse rule list
3098 -> term : ... to which the rule shall be applied
3099 -> rule : ... to be applied to term
3101 -> (rule : a rule rewriting to ...
3102 * (term : ... the resulting term ...
3103 * term list): ... with the assumptions ( //#0).
3104 ) list : there may be several such rules;
3105 the list is empty, if the rule has nothing to do
3109 fun locate_rule thy eval_rls ro [rs] t r =
3110 if (id_of_thm r) mem (map (id_of_thm)) rs
3112 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
3114 SOME ta => [(r, ta)]
3115 | NONE => (writeln("### locate_rule: rewrite "^
3116 (id_of_thm r)^" "^(term2str t)^" = NONE");
3118 else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
3119 | locate_rule _ _ _ _ _ _ =
3120 raise error ("locate_rule: doesnt match rev-sets in istate");
3122 (*.next_rule = fn : rule list -> term -> rule option
3123 for a given term return the next rules to be done for cancelling.
3125 rule list : the reverse rule list
3126 term : the term for which ...
3128 -> rule option: ... this rule is appropriate for cancellation;
3129 there may be no such rule (if the term is canceled already.*)
3130 (* val thy = Rational.thy;
3131 val Rrls {rew_ord=(_,ro),...} = cancel;
3132 val ([rs],t) = (rss,f);
3133 next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
3135 val (thy, [rs]) = (Rational.thy, revsets);
3136 val Rrls {rew_ord=(_,ro),...} = cancel;
3139 fun next_rule thy eval_rls ro [rs] t =
3140 let val der = make_deriv thy eval_rls rs ro NONE t;
3142 (* val (_,r,_)::_ = der;
3144 (_,r,_)::_ => SOME r
3147 | next_rule _ _ _ _ _ =
3148 raise error ("next_rule: doesnt match rev-sets in istate");
3150 (*.val attach_form = f : rule list -> term -> term
3151 -> (rule * (term * term list)) list
3152 checks an input term TI, if it may belong to a current cancellation, by
3153 trying to derive it from the given term TG.
3155 term : TG, the last one in the cancellation agreed upon by user + math-eng
3156 -> term: TI, the next one input by the user
3158 -> (rule : the rule to be applied in order to reach TI
3159 * (term : ... obtained by applying the rule ...
3160 * term list): ... and the respective assumptions.
3161 ) list : there may be several such rules;
3162 the list is empty, if the users term does not belong
3163 to a cancellation of the term last agreed upon.*)
3164 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
3165 []:(rule * (term * term list)) list;
3167 val pat0 = (term_of o the o (parse thy)) "?r/?s+?u/?v";
3168 val pat1 = (term_of o the o (parse thy)) "?r/?s+?u ";
3169 val pat2 = (term_of o the o (parse thy)) "?r +?u/?v";
3170 val prepat = [([HOLogic.true_const], pat0),
3171 ([HOLogic.true_const], pat1),
3172 ([HOLogic.true_const], pat2)];
3176 (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;
3177 besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht
3178 dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)
3179 val common_nominator_p =
3180 Rrls {id = "common_nominator_p", prepat=prepat,
3181 rew_ord=("ord_make_polynomial",
3182 ord_make_polynomial false Rational.thy),
3183 erls = rational_erls,
3184 calc = [("PLUS" ,("op +" ,eval_binop "#add_")),
3185 ("TIMES" ,("op *" ,eval_binop "#mult_")),
3186 ("DIVIDE" ,("HOL.divide" ,eval_cancel "#divide_")),
3187 ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
3188 scr=Rfuns {init_state = init_state thy Atools_erls ro,
3189 normal_form = add_fraction_p_ thy,(*FIXME.WN0211*)
3190 locate_rule = locate_rule thy Atools_erls ro,
3191 next_rule = next_rule thy Atools_erls ro,
3192 attach_form = attach_form}}
3196 local(*.ad [2] 'common_nominator'
3197 ---------------------------------
3198 FIXME Beschreibung .*)
3201 val {rules=rules,rew_ord=(_,ro),...} =
3202 rep_rls (assoc_rls "make_polynomial");
3203 val thy = Rational.thy;
3206 (*.common_nominator_ = fn : theory -> term -> (term * term list) option
3209 (*.init_state = fn : term -> istate
3210 initialzies the state of the interactive interpreter. The state is:
3212 type rrlsstate = (*state for reverse rewriting*)
3213 (term * (*the current formula*)
3214 term * (*the final term*)
3215 rule list (*'reverse rule list' (#)*)
3216 list * (*may be serveral, eg. in norm_rational*)
3217 (rule * (*Thm (+ Thm generated from Calc) resulting in ...*)
3218 (term * (*... rewrite with ...*)
3219 term list)) (*... assumptions*)
3220 list); (*derivation from given term to normalform
3221 in reverse order with sym_thm;
3222 (#) could be extracted from here by (map #1)*).*)
3223 fun init_state thy eval_rls ro t =
3224 let val SOME (t',_) = common_nominator_ thy t;
3225 val SOME (t'',asm) = add_fraction_ thy t;
3226 val der = reverse_deriv thy eval_rls rules ro NONE t';
3227 val der = der @ [(Thm ("real_mult_div_cancel2",
3228 num_str real_mult_div_cancel2),
3230 val rs = (distinct_Thm o (map #1)) der;
3231 val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
3233 "sym_real_mult_1"]) rs;
3234 in (t,t'',[rs(*here only _ONE_*)],der) end;
3236 (* use"knowledge/Rational.ML";
3239 (*.locate_rule = fn : rule list -> term -> rule
3240 -> (rule * (term * term list) option) list.
3241 checks a rule R for being a cancel-rule, and if it is,
3242 then return the list of rules (+ the terms they are rewriting to)
3243 which need to be applied before R should be applied.
3244 precondition: the rule is applicable to the argument-term.
3246 rule list: the reverse rule list
3247 -> term : ... to which the rule shall be applied
3248 -> rule : ... to be applied to term
3250 -> (rule : a rule rewriting to ...
3251 * (term : ... the resulting term ...
3252 * term list): ... with the assumptions ( //#0).
3253 ) list : there may be several such rules;
3254 the list is empty, if the rule has nothing to do
3258 fun locate_rule thy eval_rls ro [rs] t r =
3259 if (id_of_thm r) mem (map (id_of_thm)) rs
3261 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
3263 SOME ta => [(r, ta)]
3264 | NONE => (writeln("### locate_rule: rewrite "^
3265 (id_of_thm r)^" "^(term2str t)^" = NONE");
3267 else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[])
3268 | locate_rule _ _ _ _ _ _ =
3269 raise error ("locate_rule: doesnt match rev-sets in istate");
3271 (*.next_rule = fn : rule list -> term -> rule option
3272 for a given term return the next rules to be done for cancelling.
3274 rule list : the reverse rule list
3275 term : the term for which ...
3277 -> rule option: ... this rule is appropriate for cancellation;
3278 there may be no such rule (if the term is canceled already.*)
3279 (* val thy = Rational.thy;
3280 val Rrls {rew_ord=(_,ro),...} = cancel;
3281 val ([rs],t) = (rss,f);
3282 next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
3284 val (thy, [rs]) = (Rational.thy, revsets);
3285 val Rrls {rew_ord=(_,ro),...} = cancel_p;
3288 fun next_rule thy eval_rls ro [rs] t =
3289 let val der = make_deriv thy eval_rls rs ro NONE t;
3291 (* val (_,r,_)::_ = der;
3293 (_,r,_)::_ => SOME r
3296 | next_rule _ _ _ _ _ =
3297 raise error ("next_rule: doesnt match rev-sets in istate");
3299 (*.val attach_form = f : rule list -> term -> term
3300 -> (rule * (term * term list)) list
3301 checks an input term TI, if it may belong to a current cancellation, by
3302 trying to derive it from the given term TG.
3304 term : TG, the last one in the cancellation agreed upon by user + math-eng
3305 -> term: TI, the next one input by the user
3307 -> (rule : the rule to be applied in order to reach TI
3308 * (term : ... obtained by applying the rule ...
3309 * term list): ... and the respective assumptions.
3310 ) list : there may be several such rules;
3311 the list is empty, if the users term does not belong
3312 to a cancellation of the term last agreed upon.*)
3313 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
3314 []:(rule * (term * term list)) list;
3316 val pat0 = (term_of o the o (parse thy)) "?r/?s+?u/?v";
3317 val pat01 = (term_of o the o (parse thy)) "?r/?s-?u/?v";
3318 val pat1 = (term_of o the o (parse thy)) "?r/?s+?u ";
3319 val pat11 = (term_of o the o (parse thy)) "?r/?s-?u ";
3320 val pat2 = (term_of o the o (parse thy)) "?r +?u/?v";
3321 val pat21 = (term_of o the o (parse thy)) "?r -?u/?v";
3322 val prepat = [([HOLogic.true_const], pat0),
3323 ([HOLogic.true_const], pat01),
3324 ([HOLogic.true_const], pat1),
3325 ([HOLogic.true_const], pat11),
3326 ([HOLogic.true_const], pat2),
3327 ([HOLogic.true_const], pat21)];
3332 val common_nominator =
3333 Rrls {id = "common_nominator", prepat=prepat,
3334 rew_ord=("ord_make_polynomial",
3335 ord_make_polynomial false Rational.thy),
3336 erls = rational_erls,
3337 calc = [("PLUS" ,("op +" ,eval_binop "#add_")),
3338 ("TIMES" ,("op *" ,eval_binop "#mult_")),
3339 ("DIVIDE" ,("HOL.divide" ,eval_cancel "#divide_")),
3340 ("POWER" ,("Atools.pow" ,eval_binop "#power_"))],
3341 (*asm_thm=[("real_mult_div_cancel2","")],*)
3342 scr=Rfuns {init_state = init_state thy Atools_erls ro,
3343 normal_form = add_fraction_ (*NOT common_nominator_*) thy,
3344 locate_rule = locate_rule thy Atools_erls ro,
3345 next_rule = next_rule thy Atools_erls ro,
3346 attach_form = attach_form}}
3357 (*.the expression contains + - * ^ / only ?.*)
3358 fun is_ratpolyexp (Free _) = true
3359 | is_ratpolyexp (Const ("op +",_) $ Free _ $ Free _) = true
3360 | is_ratpolyexp (Const ("op -",_) $ Free _ $ Free _) = true
3361 | is_ratpolyexp (Const ("op *",_) $ Free _ $ Free _) = true
3362 | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true
3363 | is_ratpolyexp (Const ("HOL.divide",_) $ Free _ $ Free _) = true
3364 | is_ratpolyexp (Const ("op +",_) $ t1 $ t2) =
3365 ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
3366 | is_ratpolyexp (Const ("op -",_) $ t1 $ t2) =
3367 ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
3368 | is_ratpolyexp (Const ("op *",_) $ t1 $ t2) =
3369 ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
3370 | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =
3371 ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
3372 | is_ratpolyexp (Const ("HOL.divide",_) $ t1 $ t2) =
3373 ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
3374 | is_ratpolyexp _ = false;
3376 (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)
3377 fun eval_is_ratpolyexp (thmid:string) _
3378 (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =
3379 if is_ratpolyexp arg
3380 then SOME (mk_thmid thmid ""
3381 ((Syntax.string_of_term (thy2ctxt thy)) arg) "",
3382 Trueprop $ (mk_equality (t, HOLogic.true_const)))
3383 else SOME (mk_thmid thmid ""
3384 ((Syntax.string_of_term (thy2ctxt thy)) arg) "",
3385 Trueprop $ (mk_equality (t, HOLogic.false_const)))
3386 | eval_is_ratpolyexp _ _ _ _ = NONE;
3390 (*-------------------18.3.03 --> struct <-----------vvv--*)
3391 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)
3393 (*.discard binary minus, shift unary minus into -1*;
3394 unary minus before numerals are put into the numeral by parsing;
3395 contains absolute minimum of thms for context in norm_Rational .*)
3396 val discard_minus = prep_rls(
3397 Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
3398 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
3399 rules = [Thm ("real_diff_minus", num_str real_diff_minus),
3400 (*"a - b = a + -1 * b"*)
3401 Thm ("sym_real_mult_minus1",num_str (real_mult_minus1 RS sym))
3402 (*- ?z = "-1 * ?z"*)
3404 scr = Script ((term_of o the o (parse thy))
3407 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)
3408 val powers_erls = prep_rls(
3409 Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
3410 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
3411 rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),
3412 Calc ("Atools.is'_even",eval_is_even "#is_even_"),
3413 Calc ("op <",eval_equ "#less_"),
3414 Thm ("not_false", not_false),
3415 Thm ("not_true", not_true),
3416 Calc ("op +",eval_binop "#add_")
3418 scr = Script ((term_of o the o (parse thy))
3421 (*.all powers over + distributed; atoms over * collected, other distributed
3422 contains absolute minimum of thms for context in norm_Rational .*)
3423 val powers = prep_rls(
3424 Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
3425 erls = powers_erls, srls = Erls, calc = [], (*asm_thm = [],*)
3426 rules = [Thm ("realpow_multI", num_str realpow_multI),
3427 (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)
3428 Thm ("realpow_pow",num_str realpow_pow),
3429 (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)
3430 Thm ("realpow_oneI",num_str realpow_oneI),
3432 Thm ("realpow_minus_even",num_str realpow_minus_even),
3433 (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)
3434 Thm ("realpow_minus_odd",num_str realpow_minus_odd),
3435 (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)
3437 (*----- collect atoms over * -----*)
3438 Thm ("realpow_two_atom",num_str realpow_two_atom),
3439 (*"r is_atom ==> r * r = r ^^^ 2"*)
3440 Thm ("realpow_plus_1",num_str realpow_plus_1),
3441 (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)
3442 Thm ("realpow_addI_atom",num_str realpow_addI_atom),
3443 (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)
3445 (*----- distribute none-atoms -----*)
3446 Thm ("realpow_def_atom",num_str realpow_def_atom),
3447 (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)
3448 Thm ("realpow_eq_oneI",num_str realpow_eq_oneI),
3450 Calc ("op +",eval_binop "#add_")
3452 scr = Script ((term_of o the o (parse thy))
3455 (*.contains absolute minimum of thms for context in norm_Rational.*)
3456 val rat_mult_divide = prep_rls(
3457 Rls {id = "rat_mult_divide", preconds = [],
3458 rew_ord = ("dummy_ord",dummy_ord),
3459 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
3460 rules = [Thm ("rat_mult",num_str rat_mult),
3461 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
3462 Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),
3463 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],
3464 otherwise inv.to a / b / c = ...*)
3465 Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),
3466 (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much
3467 and does not commute a / b * c ^^^ 2 !*)
3469 Thm ("real_divide_divide1_eq", real_divide_divide1_eq),
3470 (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
3471 Thm ("real_divide_divide2_eq", real_divide_divide2_eq),
3472 (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
3473 Calc ("HOL.divide" ,eval_cancel "#divide_")
3475 scr = Script ((term_of o the o (parse thy)) "empty_script")
3477 (*.contains absolute minimum of thms for context in norm_Rational.*)
3478 val reduce_0_1_2 = prep_rls(
3479 Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
3480 erls = e_rls,srls = Erls,calc = [],(*asm_thm = [],*)
3481 rules = [(*Thm ("real_divide_1",num_str real_divide_1),
3482 "?x / 1 = ?x" unnecess.for normalform*)
3483 Thm ("real_mult_1",num_str real_mult_1),
3485 (*Thm ("real_mult_minus1",num_str real_mult_minus1),
3487 (*Thm ("real_minus_mult_cancel",num_str real_minus_mult_cancel),
3488 "- ?x * - ?y = ?x * ?y"*)
3490 Thm ("real_mult_0",num_str real_mult_0),
3492 Thm ("real_add_zero_left",num_str real_add_zero_left),
3494 (*Thm ("real_add_minus",num_str real_add_minus),
3497 Thm ("sym_real_mult_2",num_str (real_mult_2 RS sym)),
3498 (*"z1 + z1 = 2 * z1"*)
3499 Thm ("real_mult_2_assoc",num_str real_mult_2_assoc),
3500 (*"z1 + (z1 + k) = 2 * z1 + k"*)
3502 Thm ("real_0_divide",num_str real_0_divide)
3504 ], scr = EmptyScr}:rls);
3506 (*erls for calculate_Rational;
3507 make local with FIXX@ME result:term *term list WN0609???SKMG*)
3508 val norm_rat_erls = prep_rls(
3509 Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
3510 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
3511 rules = [Calc ("Atools.is'_const",eval_const "#is_const_")
3513 scr = Script ((term_of o the o (parse thy))
3516 (*.consists of rls containing the absolute minimum of thms.*)
3517 (*040209: this version has been used by RL for his equations,
3518 which is now replaced by MGs version below
3519 vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)
3520 val norm_Rational = prep_rls(
3521 Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord),
3522 erls = norm_rat_erls, srls = Erls, calc = [], (*asm_thm = [],*)
3523 rules = [(*sequence given by operator precedence*)
3526 Rls_ rat_mult_divide,
3529 (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)
3530 Rls_ order_add_mult,
3531 Rls_ collect_numerals,
3532 Rls_ add_fractions_p,
3535 scr = Script ((term_of o the o (parse thy))
3538 val norm_Rational_parenthesized = prep_rls(
3539 Seq {id = "norm_Rational_parenthesized", preconds = []:term list,
3540 rew_ord = ("dummy_ord", dummy_ord),
3541 erls = Atools_erls, srls = Erls,
3542 calc = [], (*asm_thm = [],*)
3543 rules = [Rls_ norm_Rational, (*from RL -- not the latest one*)
3544 Rls_ discard_parentheses
3550 (*-------------------18.3.03 --> struct <-----------^^^--*)
3554 theory' := overwritel (!theory', [("Rational.thy",Rational.thy)]);
3557 (*WN030318???SK: simplifies all but cancel and common_nominator*)
3558 val simplify_rational =
3559 merge_rls "simplify_rational" expand_binoms
3560 (append_rls "divide" calculate_Rational
3561 [Thm ("real_divide_1",num_str real_divide_1),
3563 Thm ("rat_mult",num_str rat_mult),
3564 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
3565 Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),
3566 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],
3567 otherwise inv.to a / b / c = ...*)
3568 Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),
3569 (*"?a / ?b * ?c = ?a * ?c / ?b"*)
3570 Thm ("add_minus",num_str add_minus),
3571 (*"?a + ?b - ?b = ?a"*)
3572 Thm ("add_minus1",num_str add_minus1),
3573 (*"?a - ?b + ?b = ?a"*)
3574 Thm ("real_divide_minus1",num_str real_divide_minus1)
3575 (*"?x / -1 = - ?x"*)
3582 (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)
3584 (* ------------------------------------------------------------------ *)
3585 (* Simplifier für beliebige Buchterme *)
3586 (* ------------------------------------------------------------------ *)
3587 (*----------------------- norm_Rational_mg ---------------------------*)
3588 (*. description of the simplifier see MG-DA.p.56ff .*)
3589 (* ------------------------------------------------------------------- *)
3590 val common_nominator_p_rls = prep_rls(
3591 Rls {id = "common_nominator_p_rls", preconds = [],
3592 rew_ord = ("dummy_ord",dummy_ord),
3593 erls = e_rls, srls = Erls, calc = [],
3595 [Rls_ common_nominator_p
3596 (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?
3597 FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)
3600 (* ------------------------------------------------------------------- *)
3601 val cancel_p_rls = prep_rls(
3602 Rls {id = "cancel_p_rls", preconds = [],
3603 rew_ord = ("dummy_ord",dummy_ord),
3604 erls = e_rls, srls = Erls, calc = [],
3607 (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)
3610 (* -------------------------------------------------------------------- *)
3611 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
3612 used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)
3613 val rat_mult_poly = prep_rls(
3614 Rls {id = "rat_mult_poly", preconds = [],
3615 rew_ord = ("dummy_ord",dummy_ord),
3616 erls = append_rls "e_rls-is_polyexp" e_rls
3617 [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
3618 srls = Erls, calc = [],
3620 [Thm ("rat_mult_poly_l",num_str rat_mult_poly_l),
3621 (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
3622 Thm ("rat_mult_poly_r",num_str rat_mult_poly_r)
3623 (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
3626 (* ------------------------------------------------------------------ *)
3627 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
3628 used in looping part norm_Rational_rls, see example DA-M02-main.p.60
3629 .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls,
3630 I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028
3632 val rat_mult_div_pow = prep_rls(
3633 Rls {id = "rat_mult_div_pow", preconds = [],
3634 rew_ord = ("dummy_ord",dummy_ord),
3636 (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls
3637 [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
3638 with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get
3639 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.
3640 thus we decided to go on with this flaw*)
3641 srls = Erls, calc = [],
3642 rules = [Thm ("rat_mult",num_str rat_mult),
3643 (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
3644 Thm ("rat_mult_poly_l",num_str rat_mult_poly_l),
3645 (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
3646 Thm ("rat_mult_poly_r",num_str rat_mult_poly_r),
3647 (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
3649 Thm ("real_divide_divide1_mg", real_divide_divide1_mg),
3650 (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)
3651 Thm ("real_divide_divide1_eq", real_divide_divide1_eq),
3652 (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
3653 Thm ("real_divide_divide2_eq", real_divide_divide2_eq),
3654 (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
3655 Calc ("HOL.divide" ,eval_cancel "#divide_"),
3657 Thm ("rat_power", num_str rat_power)
3658 (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
3660 scr = Script ((term_of o the o (parse thy)) "empty_script")
3662 (* ------------------------------------------------------------------ *)
3663 val rat_reduce_1 = prep_rls(
3664 Rls {id = "rat_reduce_1", preconds = [],
3665 rew_ord = ("dummy_ord",dummy_ord),
3666 erls = e_rls, srls = Erls, calc = [],
3667 rules = [Thm ("real_divide_1",num_str real_divide_1),
3669 Thm ("real_mult_1",num_str real_mult_1)
3672 scr = Script ((term_of o the o (parse thy)) "empty_script")
3674 (* ------------------------------------------------------------------ *)
3675 (*. looping part of norm_Rational(*_mg*) .*)
3676 val norm_Rational_rls = prep_rls(
3677 Rls {id = "norm_Rational_rls", preconds = [],
3678 rew_ord = ("dummy_ord",dummy_ord),
3679 erls = norm_rat_erls, srls = Erls, calc = [],
3680 rules = [Rls_ common_nominator_p_rls,
3681 Rls_ rat_mult_div_pow,
3682 Rls_ make_rat_poly_with_parentheses,
3683 Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)
3686 scr = Script ((term_of o the o (parse thy)) "empty_script")
3688 (* ------------------------------------------------------------------ *)
3689 (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)
3691 val norm_Rational(*_mg*) = prep_rls(
3692 Seq {id = "norm_Rational"(*_mg*), preconds = [],
3693 rew_ord = ("dummy_ord",dummy_ord),
3694 erls = norm_rat_erls, srls = Erls, calc = [],
3695 rules = [Rls_ discard_minus_,
3696 Rls_ rat_mult_poly,(* removes double fractions like a/b/c *)
3697 Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)
3698 Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)
3699 Rls_ norm_Rational_rls, (* the main rls, looping (#) *)
3700 Rls_ discard_parentheses_ (* mult only *)
3702 scr = Script ((term_of o the o (parse thy)) "empty_script")
3704 (* ------------------------------------------------------------------ *)
3707 ruleset' := overwritelthy thy (!ruleset',
3708 [("calculate_Rational", calculate_Rational),
3709 ("calc_rat_erls",calc_rat_erls),
3710 ("rational_erls", rational_erls),
3711 ("cancel_p", cancel_p),
3713 ("common_nominator_p", common_nominator_p),
3714 ("common_nominator_p_rls", common_nominator_p_rls),
3715 ("common_nominator" , common_nominator),
3716 ("discard_minus", discard_minus),
3717 ("powers_erls", powers_erls),
3719 ("rat_mult_divide", rat_mult_divide),
3720 ("reduce_0_1_2", reduce_0_1_2),
3721 ("rat_reduce_1", rat_reduce_1),
3722 ("norm_rat_erls", norm_rat_erls),
3723 ("norm_Rational", norm_Rational),
3724 ("norm_Rational_rls", norm_Rational_rls),
3725 ("norm_Rational_parenthesized", norm_Rational_parenthesized),
3726 ("rat_mult_poly", rat_mult_poly),
3727 ("rat_mult_div_pow", rat_mult_div_pow),
3728 ("cancel_p_rls", cancel_p_rls)
3731 calclist':= overwritel (!calclist',
3732 [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))
3738 (prep_pbt Rational.thy "pbl_simp_rat" [] e_pblID
3739 (["rational","simplification"],
3740 [("#Given" ,["term t_"]),
3741 ("#Where" ,["t_ is_ratpolyexp"]),
3742 ("#Find" ,["normalform n_"])
3744 append_rls "e_rls" e_rls [(*for preds in where_*)],
3746 [["simplification","of_rationals"]]));
3750 (*WN061025 this methods script is copied from (auto-generated) script
3751 of norm_Rational in order to ease repair on inform*)
3753 (prep_met Rational.thy "met_simp_rat" [] e_metID
3754 (["simplification","of_rationals"],
3755 [("#Given" ,["term t_"]),
3756 ("#Where" ,["t_ is_ratpolyexp"]),
3757 ("#Find" ,["normalform n_"])
3759 {rew_ord'="tless_true",
3761 calc = [], srls = e_rls,
3762 prls = append_rls "simplification_of_rationals_prls" e_rls
3763 [(*for preds in where_*)
3764 Calc ("Rational.is'_ratpolyexp",
3765 eval_is_ratpolyexp "")],
3766 crls = e_rls, nrls = norm_Rational_rls},
3767 "Script SimplifyScript (t_::real) = \
3768 \ ((Try (Rewrite_Set discard_minus_ False) @@ \
3769 \ Try (Rewrite_Set rat_mult_poly False) @@ \
3770 \ Try (Rewrite_Set make_rat_poly_with_parentheses False) @@ \
3771 \ Try (Rewrite_Set cancel_p_rls False) @@ \
3773 \ ((Try (Rewrite_Set common_nominator_p_rls False) @@ \
3774 \ Try (Rewrite_Set rat_mult_div_pow False) @@ \
3775 \ Try (Rewrite_Set make_rat_poly_with_parentheses False) @@\
3776 \ Try (Rewrite_Set cancel_p_rls False) @@ \
3777 \ Try (Rewrite_Set rat_reduce_1 False)))) @@ \
3778 \ Try (Rewrite_Set discard_parentheses_ False)) \
3782 (* use"../Knowledge/Rational.ML";
3783 use"Knowledge/Rational.ML";