
neuper@37906: (* rationals, i.e. fractions of multivariate polynomials over the real field
neuper@37906:    author: isac team
neuper@37906:    Copyright (c) isac team 2002
neuper@37906:    Use is subject to license terms.
neuper@37906: 
neuper@37906:    depends on Poly (and not on Atools), because 
neuper@37906:    fractions with _normalized_ polynomials are canceled, added, etc.
neuper@37906: *)
neuper@37906: 
neuper@37950: theory Rational imports Poly begin
neuper@37906: 
neuper@37906: consts
neuper@37906: 
neuper@37906:   is'_expanded   :: "real => bool" ("_ is'_expanded")     (*RL->Poly.thy*)
neuper@37906:   is'_ratpolyexp :: "real => bool" ("_ is'_ratpolyexp") 
neuper@37906: 
neuper@37950: axioms (*.not contained in Isabelle2002,
neuper@37950:           stated as axioms, TODO?: prove as theorems*)
neuper@37950: 
neuper@37978:   mult_cross:      "[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)"
neuper@37978:   mult_cross1:     "   b ~= 0            ==> (a / b = c    ) = (a     = b * c)"
neuper@37978:   mult_cross2:     "           d ~= 0    ==> (a     = c / d) = (a * d =     c)"
neuper@37950:                   
neuper@37978:   add_minus:       "a + b - b = a"(*RL->Poly.thy*)
neuper@37978:   add_minus1:      "a - b + b = a"(*RL->Poly.thy*)
neuper@37950:                   
neuper@37978:   rat_mult:        "a / b * (c / d) = a * c / (b * d)"(*?Isa02*) 
neuper@37978:   rat_mult2:       "a / b *  c      = a * c /  b     "(*?Isa02*)
neuper@37978: 
neuper@37978:   rat_mult_poly_l: "c is_polyexp ==> c * (a / b) = c * a /  b"
neuper@37978:   rat_mult_poly_r: "c is_polyexp ==> (a / b) * c = a * c /  b"
neuper@37906: 
neuper@37906: (*real_times_divide1_eq .. Isa02*) 
neuper@37978:   real_times_divide_1_eq:  "-1    * (c / d) =-1 * c /      d "
neuper@37978:   real_times_divide_num:   "a is_const ==> 
neuper@37950: 	          	   a     * (c / d) = a * c /      d "
neuper@37906: 
neuper@37978:   real_mult_div_cancel2:   "k ~= 0 ==> m * k / (n * k) = m / n"
neuper@37978: (*real_mult_div_cancel1:   "k ~= 0 ==> k * m / (k * n) = m / n"..Isa02*)
neuper@37906: 			  
neuper@37978:   real_divide_divide1:     "y ~= 0 ==> (u / v) / (y / z) = (u / v) * (z / y)"
neuper@37978:   real_divide_divide1_mg:  "y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"
neuper@37978: (*real_divide_divide2_eq:  "x / y / z = x / (y * z)"..Isa02*)
neuper@37906: 			  
neuper@37978:   rat_power:               "(a / b)^^^n = (a^^^n) / (b^^^n)"
neuper@37978: 
neuper@37978: 
neuper@37978:   rat_add:         "[| a is_const; b is_const; c is_const; d is_const |] ==> 
neuper@37950: 	           a / c + b / d = (a * d + b * c) / (c * d)"
neuper@37978:   rat_add_assoc:   "[| a is_const; b is_const; c is_const; d is_const |] ==> 
neuper@37950: 	           a / c +(b / d + e) = (a * d + b * c)/(d * c) + e"
neuper@37978:   rat_add1:        "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a / c + b / c = (a + b) / c"
neuper@37978:   rat_add1_assoc:   "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a / c + (b / c + e) = (a + b) / c + e"
neuper@37978:   rat_add2:        "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a / c + b = (a + b * c) / c"
neuper@37978:   rat_add2_assoc:  "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a / c + (b + e) = (a + b * c) / c + e"
neuper@37978:   rat_add3:        "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a + b / c = (a * c + b) / c"
neuper@37978:   rat_add3_assoc:   "[| a is_const; b is_const; c is_const |] ==> 
neuper@37950: 	           a + (b / c + e) = (a * c + b) / c + e"
neuper@37950: 
neuper@37950: text {*calculate in rationals: gcd, lcm, etc.
neuper@37950:       (c) Stefan Karnel 2002
neuper@37950:       Institute for Mathematics D and Institute for Software Technology, 
neuper@37950:       TU-Graz SS 2002 *}
neuper@37950: 
neuper@37950: text {* Remark on notions in the documentation below:
neuper@37950:     referring to the remark on 'polynomials' in Poly.sml we use
neuper@37950:     [2] 'polynomial' normalform (Polynom)
neuper@37950:     [3] 'expanded_term' normalform (Ausmultiplizierter Term),
neuper@37950:     where normalform [2] is a special case of [3], i.e. [3] implies [2].
neuper@37950:     Instead of 
neuper@37950:       'fraction with numerator and nominator both in normalform [2]'
neuper@37950:       'fraction with numerator and nominator both in normalform [3]' 
neuper@37950:     we say: 
neuper@37950:       'fraction in normalform [2]'
neuper@37950:       'fraction in normalform [3]' 
neuper@37950:     or
neuper@37950:       'fraction [2]'
neuper@37950:       'fraction [3]'.
neuper@37950:     a 'simple fraction' is a term with '/' as outmost operator and
neuper@37950:     numerator and nominator in normalform [2] or [3].
neuper@37950: *}
neuper@37950: 
neuper@37950: ML {*
neuper@37972: val thy = @{theory};
neuper@37978: val ------------------------------------------------------+ = "11111";
neuper@37972: 
neuper@37950: signature RATIONALI =
neuper@37950: sig
neuper@37950:   type mv_monom
neuper@37950:   type mv_poly 
neuper@37950:   val add_fraction_ : theory -> term -> (term * term list) option      
neuper@37950:   val add_fraction_p_ : theory -> term -> (term * term list) option       
neuper@37950:   val calculate_Rational : rls
neuper@37950:   val calc_rat_erls:rls
neuper@37950:   val cancel : rls
neuper@37950:   val cancel_ : theory -> term -> (term * term list) option    
neuper@37950:   val cancel_p : rls   
neuper@37950:   val cancel_p_ : theory -> term -> (term * term list) option
neuper@37950:   val common_nominator : rls              
neuper@37950:   val common_nominator_ : theory -> term -> (term * term list) option
neuper@37950:   val common_nominator_p : rls              
neuper@37950:   val common_nominator_p_ : theory -> term -> (term * term list) option
neuper@37950:   val eval_is_expanded : string -> 'a -> term -> theory -> 
neuper@37950: 			 (string * term) option                    
neuper@37950:   val expanded2polynomial : term -> term option
neuper@37950:   val factout_ : theory -> term -> (term * term list) option
neuper@37950:   val factout_p_ : theory -> term -> (term * term list) option
neuper@37950:   val is_expanded : term -> bool
neuper@37950:   val is_polynomial : term -> bool
neuper@37950: 
neuper@37950:   val mv_gcd : (int * int list) list -> mv_poly -> mv_poly
neuper@37950:   val mv_lcm : mv_poly -> mv_poly -> mv_poly
neuper@37950: 
neuper@37950:   val norm_expanded_rat_ : theory -> term -> (term * term list) option
neuper@37950: (*WN0602.2.6.pull into struct !!!
neuper@37950:   val norm_Rational : rls(*.normalizes an arbitrary rational term without
neuper@37950:                            roots into a simple and canceled fraction
neuper@37950:                            with normalform [2].*)
neuper@37950: *)
neuper@37950: (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME
neuper@37950:       rls               (*.normalizes an rational term [2] without
neuper@37950:                            roots into a simple and canceled fraction
neuper@37950:                            with normalform [2].*)
neuper@37950: *)
neuper@37950:   val norm_rational_ : theory -> term -> (term * term list) option
neuper@37950:   val polynomial2expanded : term -> term option
neuper@37950:   val rational_erls : 
neuper@37950:       rls             (*.evaluates an arbitrary rational term with numerals.*)
neuper@37950: 
neuper@37950: (*WN0210???SK: fehlen Funktionen, die exportiert werden sollen ? *)
neuper@37906: end
neuper@37950: 
neuper@37978: val ------------------------------------------------------++ = "22222";
neuper@37978: 
neuper@37950: (*.**************************************************************************
neuper@37950: survey on the functions
neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~
neuper@37950:  [2] 'polynomial'   :rls               | [3]'expanded_term':rls
neuper@37950: --------------------:------------------+-------------------:-----------------
neuper@37950:  factout_p_         :                  | factout_          :
neuper@37950:  cancel_p_          :                  | cancel_           :
neuper@37950:                     :cancel_p          |                   :cancel
neuper@37950: --------------------:------------------+-------------------:-----------------
neuper@37950:  common_nominator_p_:                  | common_nominator_ :
neuper@37950:                     :common_nominator_p|                   :common_nominator
neuper@37950:  add_fraction_p_    :                  | add_fraction_     :
neuper@37950: --------------------:------------------+-------------------:-----------------
neuper@37950: ???SK                 :norm_rational_p   |                   :norm_rational
neuper@37950: 
neuper@37950: This survey shows only the principal functions for reuse, and the identifiers 
neuper@37950: of the rls exported. The list below shows some more useful functions.
neuper@37950: 
neuper@37950: 
neuper@37950: conversion from Isabelle-term to internal representation
neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neuper@37950: 
neuper@37950: ... BITTE FORTSETZEN ...
neuper@37950: 
neuper@37950: polynomial2expanded = ...
neuper@37950: expanded2polynomial = ...
neuper@37950: 
neuper@37950: remark: polynomial2expanded o expanded2polynomial = I, 
neuper@37950:         where 'o' is function chaining, and 'I' is identity WN0210???SK
neuper@37950: 
neuper@37950: functions for greatest common divisor and canceling
neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neuper@37950: mv_gcd
neuper@37950: factout_
neuper@37950: factout_p_
neuper@37950: cancel_
neuper@37950: cancel_p_
neuper@37950: 
neuper@37950: functions for least common multiple and addition of fractions
neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neuper@37950: mv_lcm
neuper@37950: common_nominator_
neuper@37950: common_nominator_p_
neuper@37950: add_fraction_       (*.add 2 or more fractions.*)
neuper@37950: add_fraction_p_     (*.add 2 or more fractions.*)
neuper@37950: 
neuper@37950: functions for normalform of rationals
neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neuper@37950: WN0210???SK interne Funktionen f"ur norm_rational: 
neuper@37950:           schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?
neuper@37950: 
neuper@37950: norm_rational_
neuper@37950: norm_expanded_rat_
neuper@37950: 
neuper@37950: **************************************************************************.*)
neuper@37950: 
neuper@37950: 
neuper@37950: (*##*)
neuper@37950: structure RationalI : RATIONALI = 
neuper@37950: struct 
neuper@37950: (*##*)
neuper@37950: 
neuper@37950: infix mem ins union; (*WN100819 updating to Isabelle2009-2*)
neuper@37950: fun x mem [] = false
neuper@37950:   | x mem (y :: ys) = x = y orelse x mem ys;
neuper@37950: fun (x ins xs) = if x mem xs then xs else x :: xs;
neuper@37950: fun xs union [] = xs
neuper@37950:   | [] union ys = ys
neuper@37950:   | (x :: xs) union ys = xs union (x ins ys);
neuper@37950: 
neuper@37950: (*. gcd of integers .*)
neuper@37950: (* die gcd Funktion von Isabelle funktioniert nicht richtig !!! *)
neuper@37950: fun gcd_int a b = if b=0 then a
neuper@37950: 		  else gcd_int b (a mod b);
neuper@37950: 
neuper@37950: (*. univariate polynomials (uv) .*)
neuper@37950: (*. univariate polynomials are represented as a list 
neuper@37950:     of the coefficent in reverse maximum degree order .*)
neuper@37950: (*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*)
neuper@37950: type uv_poly = int list;
neuper@37950: 
neuper@37950: (*. adds two uv polynomials .*)
neuper@37950: fun uv_mod_add_poly ([]:uv_poly,p2:uv_poly) = p2:uv_poly 
neuper@37950:   | uv_mod_add_poly (p1,[]) = p1
neuper@37950:   | uv_mod_add_poly (x::p1,y::p2) = (x+y)::(uv_mod_add_poly(p1,p2)); 
neuper@37950: 
neuper@37950: (*. multiplies a uv polynomial with a skalar s .*)
neuper@37950: fun uv_mod_smul_poly ([]:uv_poly,s:int) = []:uv_poly 
neuper@37950:   | uv_mod_smul_poly (x::p,s) = (x*s)::(uv_mod_smul_poly(p,s)); 
neuper@37950: 
neuper@37950: (*. calculates the remainder of a polynomial divided by a skalar s .*)
neuper@37950: fun uv_mod_rem_poly ([]:uv_poly,s) = []:uv_poly 
neuper@37950:   | uv_mod_rem_poly (x::p,s) = (x mod s)::(uv_mod_smul_poly(p,s)); 
neuper@37950: 
neuper@37950: (*. calculates the degree of a uv polynomial .*)
neuper@37950: fun uv_mod_deg ([]:uv_poly) = 0  
neuper@37950:   | uv_mod_deg p = length(p)-1;
neuper@37950: 
neuper@37950: (*. calculates the remainder of x/p and represents it as 
neuper@37950:     value between -p/2 and p/2 .*)
neuper@37950: fun uv_mod_mod2(x,p)=
neuper@37950:     let
neuper@37950: 	val y=(x mod p);
neuper@37950:     in
neuper@37950: 	if (y)>(p div 2) then (y)-p else 
neuper@37950: 	    (
neuper@37950: 	     if (y)<(~p div 2) then p+(y) else (y)
neuper@37950: 	     )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*.calculates the remainder for each element of a integer list divided by p.*)  
neuper@37950: fun uv_mod_list_modp [] p = [] 
neuper@37950:   | uv_mod_list_modp (x::xs) p = (uv_mod_mod2(x,p))::(uv_mod_list_modp xs p);
neuper@37950: 
neuper@37950: (*. appends an integer at the end of a integer list .*)
neuper@37950: fun uv_mod_null (p1:int list,0) = p1 
neuper@37950:   | uv_mod_null (p1:int list,n1:int) = uv_mod_null(p1,n1-1) @ [0];
neuper@37950: 
neuper@37950: (*. uv polynomial division, result is (quotient, remainder) .*)
neuper@37950: (*. only for uv_mod_divides .*)
neuper@37950: (* FIXME: Division von x^9+x^5+1 durch x-1000 funktioniert nicht,
neuper@37950:    integer zu klein  *)
neuper@37950: fun uv_mod_pdiv (p1:uv_poly) ([]:uv_poly) = 
neuper@37950:     raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")
neuper@37950:   | uv_mod_pdiv p1 [x] = 
neuper@37950:     let
neuper@37950: 	val xs=ref [];
neuper@37950:     in
neuper@37950: 	if x<>0 then 
neuper@37950: 	    (
neuper@37950: 	     xs:=(uv_mod_rem_poly(p1,x));
neuper@37950: 	     while length(!xs)>0 andalso hd(!xs)=0 do xs:=tl(!xs)
neuper@37950: 	     )
neuper@37950: 	else raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero");
neuper@37950: 	([]:uv_poly,!xs:uv_poly)
neuper@37950:     end
neuper@37950:   | uv_mod_pdiv p1 p2 =  
neuper@37950:     let
neuper@37950: 	val n= uv_mod_deg(p2);
neuper@37950: 	val m= ref (uv_mod_deg(p1));
neuper@37950: 	val p1'=ref (rev(p1));
neuper@37950: 	val p2'=(rev(p2));
neuper@37950: 	val lc2=hd(p2');
neuper@37950: 	val q=ref [];
neuper@37950: 	val c=ref 0;
neuper@37950: 	val output=ref ([],[]);
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 if (!m)=0 orelse p2=[0] 
neuper@37950:          then raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero") 
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      if (!m)<n then 
neuper@37950: 		  (
neuper@37950: 		   output:=([0],p1) 
neuper@37950: 		   ) 
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   while (!m)>=n do
neuper@37950: 		       (
neuper@37950: 			c:=hd(!p1') div hd(p2');
neuper@37950: 			if !c<>0 then
neuper@37950: 			    (
neuper@37950: 			     p1':=uv_mod_add_poly(!p1',uv_mod_null(uv_mod_smul_poly(p2',~(!c)),!m-n));
neuper@37950: 			     while length(!p1')>0 andalso hd(!p1')=0  do p1':= tl(!p1');
neuper@37950: 			     m:=uv_mod_deg(!p1')
neuper@37950: 			     )
neuper@37950: 			else m:=0
neuper@37950: 			);
neuper@37950:     		   output:=(rev(!q),rev(!p1'))
neuper@37950: 		   )
neuper@37950: 	      );
neuper@37950: 	     !output
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. divides p1 by p2 in Zp .*)
neuper@37950: fun uv_mod_pdivp (p1:uv_poly) (p2:uv_poly) p =  
neuper@37950:     let
neuper@37950: 	val n=uv_mod_deg(p2);
neuper@37950: 	val m=ref (uv_mod_deg(uv_mod_list_modp p1 p));
neuper@37950: 	val p1'=ref (rev(p1));
neuper@37950: 	val p2'=(rev(uv_mod_list_modp p2 p));
neuper@37950: 	val lc2=hd(p2');
neuper@37950: 	val q=ref [];
neuper@37950: 	val c=ref 0;
neuper@37950: 	val output=ref ([],[]);
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 if (!m)=0 orelse p2=[0] then raise error ("RATIONALS_UV_MOD_PDIVP_EXCEPTION: Division by zero") 
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      if (!m)<n then 
neuper@37950: 		  (
neuper@37950: 		   output:=([0],p1) 
neuper@37950: 		   ) 
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   while (!m)>=n do
neuper@37950: 		       (
neuper@37950: 			c:=uv_mod_mod2(hd(!p1')*(power lc2 1), p);
neuper@37950: 			q:=(!c)::(!q);
neuper@37950: 			p1':=uv_mod_list_modp(tl(uv_mod_add_poly(uv_mod_smul_poly(!p1',lc2),
neuper@37950: 								  uv_mod_smul_poly(uv_mod_smul_poly(p2',hd(!p1')),~1)))) p;
neuper@37950: 			m:=(!m)-1
neuper@37950: 			);
neuper@37950: 		   
neuper@37950: 		   while !p1'<>[] andalso hd(!p1')=0 do
neuper@37950: 		       (
neuper@37950: 			p1':=tl(!p1')
neuper@37950: 			);
neuper@37950: 
neuper@37950:     		   output:=(rev(uv_mod_list_modp (!q) (p)),rev(!p1'))
neuper@37950: 		   )
neuper@37950: 	      );
neuper@37950: 	     !output:uv_poly * uv_poly
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the remainder of p1/p2 .*)
neuper@37950: fun uv_mod_prest (p1:uv_poly) ([]:uv_poly) = raise error("UV_MOD_PREST_EXCEPTION: Division by zero") 
neuper@37950:   | uv_mod_prest [] p2 = []:uv_poly
neuper@37950:   | uv_mod_prest p1 p2 = (#2(uv_mod_pdiv p1 p2));
neuper@37950: 
neuper@37950: (*. calculates the remainder of p1/p2 in Zp .*)
neuper@37950: fun uv_mod_prestp (p1:uv_poly) ([]:uv_poly) p= raise error("UV_MOD_PRESTP_EXCEPTION: Division by zero") 
neuper@37950:   | uv_mod_prestp [] p2 p= []:uv_poly 
neuper@37950:   | uv_mod_prestp p1 p2 p = #2(uv_mod_pdivp p1 p2 p); 
neuper@37950: 
neuper@37950: (*. calculates the content of a uv polynomial .*)
neuper@37950: fun uv_mod_cont ([]:uv_poly) = 0  
neuper@37950:   | uv_mod_cont (x::p)= gcd_int x (uv_mod_cont(p));
neuper@37950: 
neuper@37950: (*. divides each coefficient of a uv polynomial by y .*)
neuper@37950: fun uv_mod_div_list (p:uv_poly,0) = raise error("UV_MOD_DIV_LIST_EXCEPTION: Division by zero") 
neuper@37950:   | uv_mod_div_list ([],y)   = []:uv_poly
neuper@37950:   | uv_mod_div_list (x::p,y) = (x div y)::uv_mod_div_list(p,y); 
neuper@37950: 
neuper@37950: (*. calculates the primitiv part of a uv polynomial .*)
neuper@37950: fun uv_mod_pp ([]:uv_poly) = []:uv_poly
neuper@37950:   | uv_mod_pp p =  
neuper@37950:     let
neuper@37950: 	val c=ref 0;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 c:=uv_mod_cont(p);
neuper@37950: 	 
neuper@37950: 	 if !c=0 then raise error ("RATIONALS_UV_MOD_PP_EXCEPTION: content is 0")
neuper@37950: 	 else uv_mod_div_list(p,!c)
neuper@37950: 	)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. gets the leading coefficient of a uv polynomial .*)
neuper@37950: fun uv_mod_lc ([]:uv_poly) = 0 
neuper@37950:   | uv_mod_lc p  = hd(rev(p)); 
neuper@37950: 
neuper@37950: (*. calculates the euklidean polynomial remainder sequence in Zp .*)
neuper@37950: fun uv_mod_prs_euklid_p(p1:uv_poly,p2:uv_poly,p)= 
neuper@37950:     let
neuper@37950: 	val f =ref [];
neuper@37950: 	val f'=ref p2;
neuper@37950: 	val fi=ref [];
neuper@37950:     in
neuper@37950: 	( 
neuper@37950: 	 f:=p2::p1::[]; 
neuper@37950:  	 while uv_mod_deg(!f')>0 do
neuper@37950: 	     (
neuper@37950: 	      f':=uv_mod_prestp (hd(tl(!f))) (hd(!f)) p;
neuper@37950: 	      if (!f')<>[] then 
neuper@37950: 		  (
neuper@37950: 		   fi:=(!f');
neuper@37950: 		   f:=(!fi)::(!f)
neuper@37950: 		   )
neuper@37950: 	      else ()
neuper@37950: 	      );
neuper@37950: 	      (!f)
neuper@37950: 	 
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the gcd of p1 and p2 in Zp .*)
neuper@37950: fun uv_mod_gcd_modp ([]:uv_poly) (p2:uv_poly) p = p2:uv_poly 
neuper@37950:   | uv_mod_gcd_modp p1 [] p= p1
neuper@37950:   | uv_mod_gcd_modp p1 p2 p=
neuper@37950:     let
neuper@37950: 	val p1'=ref[];
neuper@37950: 	val p2'=ref[];
neuper@37950: 	val pc=ref[];
neuper@37950: 	val g=ref [];
neuper@37950: 	val d=ref 0;
neuper@37950: 	val prs=ref [];
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 if uv_mod_deg(p1)>=uv_mod_deg(p2) then
neuper@37950: 	     (
neuper@37950: 	      p1':=uv_mod_list_modp (uv_mod_pp(p1)) p;
neuper@37950: 	      p2':=uv_mod_list_modp (uv_mod_pp(p2)) p
neuper@37950: 	      )
neuper@37950: 	 else 
neuper@37950: 	     (
neuper@37950: 	      p1':=uv_mod_list_modp (uv_mod_pp(p2)) p;
neuper@37950: 	      p2':=uv_mod_list_modp (uv_mod_pp(p1)) p
neuper@37950: 	      );
neuper@37950: 	 d:=uv_mod_mod2((gcd_int (uv_mod_cont(p1))) (uv_mod_cont(p2)), p) ;
neuper@37950: 	 if !d>(p div 2) then d:=(!d)-p else ();
neuper@37950: 	 
neuper@37950: 	 prs:=uv_mod_prs_euklid_p(!p1',!p2',p);
neuper@37950: 
neuper@37950: 	 if hd(!prs)=[] then pc:=hd(tl(!prs))
neuper@37950: 	 else pc:=hd(!prs);
neuper@37950: 
neuper@37950: 	 g:=uv_mod_smul_poly(uv_mod_pp(!pc),!d);
neuper@37950: 	 !g
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the minimum of two real values x and y .*)
neuper@37978: fun uv_mod_r_min(x,y):Real.real = if x>y then y else x;
neuper@37950: 
neuper@37950: (*. calculates the minimum of two integer values x and y .*)
neuper@37950: fun uv_mod_min(x,y) = if x>y then y else x;
neuper@37950: 
neuper@37950: (*. adds the squared values of a integer list .*)
neuper@37950: fun uv_mod_add_qu [] = 0.0 
neuper@37978:   | uv_mod_add_qu (x::p) =  Real.fromInt(x)*Real.fromInt(x) + uv_mod_add_qu p;
neuper@37950: 
neuper@37950: (*. calculates the euklidean norm .*)
neuper@37950: fun uv_mod_norm ([]:uv_poly) = 0.0
neuper@37950:   | uv_mod_norm p = Math.sqrt(uv_mod_add_qu(p));
neuper@37950: 
neuper@37950: (*. multipies two values a and b .*)
neuper@37950: fun uv_mod_multi a b = a * b;
neuper@37950: 
neuper@37950: (*. decides if x is a prim, the list contains all primes which are lower then x .*)
neuper@37950: fun uv_mod_prim(x,[])= false 
neuper@37950:   | uv_mod_prim(x,[y])=if ((x mod y) <> 0) then true
neuper@37950: 		else false
neuper@37950:   | uv_mod_prim(x,y::ys) = if uv_mod_prim(x,[y])
neuper@37950: 			then 
neuper@37950: 			    if uv_mod_prim(x,ys) then true 
neuper@37950: 			    else false
neuper@37950: 		    else false;
neuper@37950: 
neuper@37950: (*. gets the first prime, which is greater than p and does not divide g .*)
neuper@37950: fun uv_mod_nextprime(g,p)= 
neuper@37950:     let
neuper@37950: 	val list=ref [2];
neuper@37950: 	val exit=ref 0;
neuper@37950: 	val i = ref 2
neuper@37950:     in
neuper@37950: 	while (!i<p) do (* calculates the primes lower then p *)
neuper@37950: 	    (
neuper@37950: 	     if uv_mod_prim(!i,!list) then
neuper@37950: 		 (
neuper@37950: 		  if (p mod !i <> 0)
neuper@37950: 		      then
neuper@37950: 			  (
neuper@37950: 			   list:= (!i)::(!list);
neuper@37950: 			   i:= (!i)+1
neuper@37950: 			   )
neuper@37950: 		  else i:=(!i)+1
neuper@37950: 		  )
neuper@37950: 	     else i:= (!i)+1
neuper@37950: 		 );
neuper@37950: 	    i:=(p+1);
neuper@37950: 	    while (!exit=0) do   (* calculate next prime which does not divide g *)
neuper@37950: 	    (
neuper@37950: 	     if uv_mod_prim(!i,!list) then
neuper@37950: 		 (
neuper@37950: 		  if (g mod !i <> 0)
neuper@37950: 		      then
neuper@37950: 			  (
neuper@37950: 			   list:= (!i)::(!list);
neuper@37950: 			   exit:= (!i)
neuper@37950: 			   )
neuper@37950: 		  else i:=(!i)+1
neuper@37950: 		      )
neuper@37950: 	     else i:= (!i)+1
neuper@37950: 		 ); 
neuper@37950: 	    !exit
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. decides if p1 is a factor of p2 in Zp .*)
neuper@37950: fun uv_mod_dividesp ([]:uv_poly) (p2:uv_poly) p= raise error("UV_MOD_DIVIDESP: Division by zero") 
neuper@37950:   | uv_mod_dividesp p1 p2 p= if uv_mod_prestp p2 p1 p = [] then true else false;
neuper@37950: 
neuper@37950: (*. decides if p1 is a factor of p2 .*)
neuper@37950: fun uv_mod_divides ([]:uv_poly) (p2:uv_poly) = raise error("UV_MOD_DIVIDES: Division by zero")
neuper@37950:   | uv_mod_divides p1 p2 = if uv_mod_prest p2 p1  = [] then true else false;
neuper@37950: 
neuper@37950: (*. chinese remainder algorithm .*)
neuper@37950: fun uv_mod_cra2(r1,r2,m1,m2)=     
neuper@37950:     let 
neuper@37950: 	val c=ref 0;
neuper@37950: 	val r1'=ref 0;
neuper@37950: 	val d=ref 0;
neuper@37950: 	val a=ref 0;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 while (uv_mod_mod2((!c)*m1,m2))<>1 do 
neuper@37950: 	     (
neuper@37950: 	      c:=(!c)+1
neuper@37950: 	      );
neuper@37950: 	 r1':= uv_mod_mod2(r1,m1);
neuper@37950: 	 d:=uv_mod_mod2(((r2-(!r1'))*(!c)),m2);
neuper@37950: 	 !r1'+(!d)*m1    
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. applies the chinese remainder algorithmen to the coefficients of x1 and x2 .*)
neuper@37950: fun uv_mod_cra_2 ([],[],m1,m2) = [] 
neuper@37950:   | uv_mod_cra_2 ([],x2,m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x1")
neuper@37950:   | uv_mod_cra_2 (x1,[],m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x2")
neuper@37950:   | 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));
neuper@37950: 
neuper@37950: (*. calculates the gcd of two uv polynomials p1' and p2' with the modular algorithm .*)
neuper@37950: fun uv_mod_gcd (p1':uv_poly) (p2':uv_poly) =
neuper@37950:     let 
neuper@37950: 	val p1=ref (uv_mod_pp(p1'));
neuper@37950: 	val p2=ref (uv_mod_pp(p2'));
neuper@37950: 	val c=gcd_int (uv_mod_cont(p1')) (uv_mod_cont(p2'));
neuper@37950: 	val temp=ref [];
neuper@37950: 	val cp=ref [];
neuper@37950: 	val qp=ref [];
neuper@37950: 	val q=ref[];
neuper@37950: 	val pn=ref 0;
neuper@37950: 	val d=ref 0;
neuper@37950: 	val g1=ref 0;
neuper@37950: 	val p=ref 0;    
neuper@37950: 	val m=ref 0;
neuper@37950: 	val exit=ref 0;
neuper@37950: 	val i=ref 1;
neuper@37950:     in
neuper@37950: 	if length(!p1)>length(!p2) then ()
neuper@37950: 	else 
neuper@37950: 	    (
neuper@37950: 	     temp:= !p1;
neuper@37950: 	     p1:= !p2;
neuper@37950: 	     p2:= !temp
neuper@37950: 	     );
neuper@37950: 
neuper@37950: 	 
neuper@37950: 	d:=gcd_int (uv_mod_lc(!p1)) (uv_mod_lc(!p2));
neuper@37950: 	g1:=uv_mod_lc(!p1)*uv_mod_lc(!p2);
neuper@37950: 	p:=4;
neuper@37950: 	
neuper@37978: 	m:=Real.ceil(2.0 * Real.fromInt(!d) *
neuper@37978: 	  Real.fromInt(power 2 (uv_mod_min(uv_mod_deg(!p2),uv_mod_deg(!p1)))) *
neuper@37978: 	  Real.fromInt(!d) * 
neuper@37978: 	  uv_mod_r_min(uv_mod_norm(!p1) / Real.fromInt(abs(uv_mod_lc(!p1))),
neuper@37978: 	  uv_mod_norm(!p2) / Real.fromInt(abs(uv_mod_lc(!p2))))); 
neuper@37950: 
neuper@37950: 	while (!exit=0) do  
neuper@37950: 	    (
neuper@37950: 	     p:=uv_mod_nextprime(!d,!p);
neuper@37950: 	     cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)) ;
neuper@37950: 	     if abs(uv_mod_lc(!cp))<>1 then  (* leading coefficient = 1 ? *)
neuper@37950: 		 (
neuper@37950: 		  i:=1;
neuper@37950: 		  while (!i)<(!p) andalso (abs(uv_mod_mod2((uv_mod_lc(!cp)*(!i)),(!p)))<>1) do
neuper@37950: 		      (
neuper@37950: 		       i:=(!i)+1
neuper@37950: 		       );
neuper@37950: 		      cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p) 
neuper@37950: 		  )
neuper@37950: 	     else ();
neuper@37950: 
neuper@37950: 	     qp:= ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp));
neuper@37950: 
neuper@37950: 	     if uv_mod_deg(!qp)=0 then (q:=[1]; exit:=1) else ();
neuper@37950: 
neuper@37950: 	     pn:=(!p);
neuper@37950: 	     q:=(!qp);
neuper@37950: 
neuper@37950: 	     while !pn<= !m andalso !m>(!p) andalso !exit=0 do
neuper@37950: 		 (
neuper@37950: 		  p:=uv_mod_nextprime(!d,!p);
neuper@37950:  		  cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)); 
neuper@37950:  		  if uv_mod_lc(!cp)<>1 then  (* leading coefficient = 1 ? *)
neuper@37950:  		      (
neuper@37950:  		       i:=1;
neuper@37950:  		       while (!i)<(!p) andalso ((uv_mod_mod2((uv_mod_lc(!q)*(!i)),(!p)))<>1) do
neuper@37950:  			   (
neuper@37950:  			    i:=(!i)+1
neuper@37950: 		           );
neuper@37950: 		       cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)
neuper@37950:  		      )
neuper@37950:  		  else ();    
neuper@37950:  		 
neuper@37950: 		  qp:=uv_mod_list_modp ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp)  ) (!p);
neuper@37950:  		  if uv_mod_deg(!qp)>uv_mod_deg(!q) then
neuper@37950:  		      (
neuper@37950:  		       (*print("degree to high!!!\n")*)
neuper@37950:  		       )
neuper@37950:  		  else
neuper@37950:  		      (
neuper@37950:  		       if uv_mod_deg(!qp)=uv_mod_deg(!q) then
neuper@37950:  			   (
neuper@37950:  			    q:=uv_mod_cra_2(!q,!qp,!pn,!p);
neuper@37950: 			    pn:=(!pn) * !p;
neuper@37950: 			    q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn)); (* found already gcd ? *)
neuper@37950: 			    if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then (exit:=1) else ()
neuper@37950: 		 	    )
neuper@37950: 		       else
neuper@37950: 			   (
neuper@37950: 			    if  uv_mod_deg(!qp)<uv_mod_deg(!q) then
neuper@37950: 				(
neuper@37950: 				 pn:= !p;
neuper@37950: 				 q:= !qp
neuper@37950: 				 )
neuper@37950: 			    else ()
neuper@37950: 			    )
neuper@37950: 		       )
neuper@37950: 		  );
neuper@37950:  	     q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn));
neuper@37950: 	     if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then exit:=1 else ()
neuper@37950: 	     );
neuper@37950: 	    uv_mod_smul_poly(!q,c):uv_poly
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. multivariate polynomials .*)
neuper@37950: (*. multivariate polynomials are represented as a list of the pairs, 
neuper@37950:  first is the coefficent and the second is a list of the exponents .*)
neuper@37950: (*. 5 * x^5 * y^3 + 4 * x^3 * z^2 + 2 * x^2 * y * z^3 - z - 19 
neuper@37950:  => [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*)
neuper@37950: 
neuper@37950: (*. global variables .*)
neuper@37950: (*. order indicators .*)
neuper@37950: val LEX_=0; (* lexicographical term order *)
neuper@37950: val GGO_=1; (* greatest degree order *)
neuper@37950: 
neuper@37950: (*. datatypes for internal representation.*)
neuper@37950: type mv_monom = (int *      (*.coefficient or the monom.*)
neuper@37950: 		 int list); (*.list of exponents)      .*)
neuper@37950: fun mv_monom2str (i, is) = "("^ int2str i^"," ^ ints2str' is ^ ")";
neuper@37950: 
neuper@37950: type mv_poly = mv_monom list; 
neuper@37950: fun mv_poly2str p = (strs2str' o (map mv_monom2str)) p;
neuper@37950: 
neuper@37950: (*. help function for monom_greater and geq .*)
neuper@37950: fun mv_mg_hlp([]) = EQUAL 
neuper@37950:   | mv_mg_hlp(x::list)=if x<0 then LESS
neuper@37950: 		    else if x>0 then GREATER
neuper@37950: 			 else mv_mg_hlp(list);
neuper@37950: 
neuper@37950: (*. adds a list of values .*)
neuper@37950: fun mv_addlist([]) = 0
neuper@37950:   | mv_addlist(p1) = hd(p1)+mv_addlist(tl(p1));
neuper@37950: 			   
neuper@37950: (*. tests if the monomial M1 is greater as the monomial M2 and returns a boolean value .*)
neuper@37950: (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
neuper@37950: fun mv_monom_greater((M1x,M1l):mv_monom,(M2x,M2l):mv_monom,order)=
neuper@37950:     if order=LEX_ then
neuper@37950: 	( 
neuper@37950: 	 if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
neuper@37950: 	 else if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
neuper@37950: 	     )
neuper@37950:     else
neuper@37950: 	if order=GGO_ then
neuper@37950: 	    ( 
neuper@37950: 	     if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")
neuper@37950: 	     else 
neuper@37950: 		 if mv_addlist(M1l)=mv_addlist(M2l)  then if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true
neuper@37950: 		 else if mv_addlist(M1l)>mv_addlist(M2l) then true else false
neuper@37950: 	     )
neuper@37950: 	else raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");
neuper@37950: 		   
neuper@37950: (*. tests if the monomial X is greater as the monomial Y and returns a order value (GREATER,EQUAL,LESS) .*)
neuper@37950: (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)
neuper@37950: fun mv_geq order ((x1,x):mv_monom,(x2,y):mv_monom) =
neuper@37950: let 
neuper@37950:     val temp=ref EQUAL;
neuper@37950: in
neuper@37950:     if order=LEX_ then
neuper@37950: 	(
neuper@37950: 	 if length(x)<>length(y) then 
neuper@37950: 	     raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
neuper@37950: 	 else 
neuper@37950: 	     (
neuper@37950: 	      temp:=mv_mg_hlp((map op- (x~~y)));
neuper@37950: 	      if !temp=EQUAL then 
neuper@37950: 		  ( if x1=x2 then EQUAL 
neuper@37950: 		    else if x1>x2 then GREATER
neuper@37950: 			 else LESS
neuper@37950: 			     )
neuper@37950: 	      else (!temp)
neuper@37950: 	      )
neuper@37950: 	     )
neuper@37950:     else 
neuper@37950: 	if order=GGO_ then 
neuper@37950: 	    (
neuper@37950: 	     if length(x)<>length(y) then 
neuper@37950: 	      raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")
neuper@37950: 	     else 
neuper@37950: 		 if mv_addlist(x)=mv_addlist(y) then 
neuper@37950: 		     (mv_mg_hlp((map op- (x~~y))))
neuper@37950: 		 else if mv_addlist(x)>mv_addlist(y) then GREATER else LESS
neuper@37950: 		     )
neuper@37950: 	else raise error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")
neuper@37950: end;
neuper@37950: 
neuper@37950: (*. cuts the first variable from a polynomial .*)
neuper@37950: fun mv_cut([]:mv_poly)=[]:mv_poly
neuper@37950:   | mv_cut((x,[])::list) = raise error ("RATIONALS_MV_CUT_EXCEPTION: Invalid list ")
neuper@37950:   | mv_cut((x,y::ys)::list)=(x,ys)::mv_cut(list);
neuper@37950: 	    
neuper@37950: (*. leading power product .*)
neuper@37950: fun mv_lpp([]:mv_poly,order)  = []
neuper@37950:   | mv_lpp([(x,y)],order) = y
neuper@37950:   | mv_lpp(p1,order)  = #2(hd(rev(sort (mv_geq order) p1)));
neuper@37950:     
neuper@37950: (*. leading monomial .*)
neuper@37950: fun mv_lm([]:mv_poly,order)  = (0,[]):mv_monom
neuper@37950:   | mv_lm([x],order) = x 
neuper@37950:   | mv_lm(p1,order)  = hd(rev(sort (mv_geq order) p1));
neuper@37950:     
neuper@37950: (*. leading coefficient in term order .*)
neuper@37950: fun mv_lc2([]:mv_poly,order)  = 0
neuper@37950:   | mv_lc2([(x,y)],order) = x
neuper@37950:   | mv_lc2(p1,order)  = #1(hd(rev(sort (mv_geq order) p1)));
neuper@37950: 
neuper@37950: 
neuper@37950: (*. reverse the coefficients in mv polynomial .*)
neuper@37950: fun mv_rev_to([]:mv_poly) = []:mv_poly
neuper@37950:   | mv_rev_to((c,e)::xs) = (c,rev(e))::mv_rev_to(xs);
neuper@37950: 
neuper@37950: (*. leading coefficient in reverse term order .*)
neuper@37950: fun mv_lc([]:mv_poly,order)  = []:mv_poly 
neuper@37950:   | mv_lc([(x,y)],order) = mv_rev_to(mv_cut(mv_rev_to([(x,y)])))
neuper@37950:   | mv_lc(p1,order)  = 
neuper@37950:     let
neuper@37950: 	val p1o=ref (rev(sort (mv_geq order) (mv_rev_to(p1))));
neuper@37950: 	val lp=hd(#2(hd(!p1o)));
neuper@37950: 	val lc=ref [];
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 while (length(!p1o)>0 andalso hd(#2(hd(!p1o)))=lp) do
neuper@37950: 	     (
neuper@37950: 	      lc:=hd(mv_cut([hd(!p1o)]))::(!lc);
neuper@37950: 	      p1o:=tl(!p1o)
neuper@37950: 	      );
neuper@37950: 	 if !lc=[] then raise error ("RATIONALS_MV_LC_EXCEPTION: lc is empty") else ();
neuper@37950: 	 mv_rev_to(!lc)
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. compares two powerproducts .*)
neuper@37950: fun mv_monom_equal((_,xlist):mv_monom,(_,ylist):mv_monom) = (foldr and_) (((map op=) (xlist~~ylist)),true);
neuper@37950:     
neuper@37950: (*. help function for mv_add .*)
neuper@37950: fun mv_madd([]:mv_poly,[]:mv_poly,order) = []:mv_poly
neuper@37950:   | mv_madd([(0,_)],p2,order) = p2
neuper@37950:   | mv_madd(p1,[(0,_)],order) = p1  
neuper@37950:   | mv_madd([],p2,order) = p2
neuper@37950:   | mv_madd(p1,[],order) = p1
neuper@37950:   | mv_madd(p1,p2,order) = 
neuper@37950:     (
neuper@37950:      if mv_monom_greater(hd(p1),hd(p2),order) 
neuper@37950: 	 then hd(p1)::mv_madd(tl(p1),p2,order)
neuper@37950:      else if mv_monom_equal(hd(p1),hd(p2)) 
neuper@37950: 	      then if mv_lc2(p1,order)+mv_lc2(p2,order)<>0 
neuper@37950: 		       then (mv_lc2(p1,order)+mv_lc2(p2,order),mv_lpp(p1,order))::mv_madd(tl(p1),tl(p2),order)
neuper@37950: 		   else mv_madd(tl(p1),tl(p2),order)
neuper@37950: 	  else hd(p2)::mv_madd(p1,tl(p2),order)
neuper@37950: 	      )
neuper@37950: 	      
neuper@37950: (*. adds two multivariate polynomials .*)
neuper@37950: fun mv_add([]:mv_poly,p2:mv_poly,order) = p2
neuper@37950:   | mv_add(p1,[],order) = p1
neuper@37950:   | mv_add(p1,p2,order) = mv_madd(rev(sort (mv_geq order) p1),rev(sort (mv_geq order) p2), order);
neuper@37950: 
neuper@37950: (*. monom multiplication .*)
neuper@37950: fun mv_mmul((x1,y1):mv_monom,(x2,y2):mv_monom)=(x1*x2,(map op+) (y1~~y2)):mv_monom;
neuper@37950: 
neuper@37950: (*. deletes all monomials with coefficient 0 .*)
neuper@37950: fun mv_shorten([]:mv_poly,order) = []:mv_poly
neuper@37950:   | mv_shorten(x::xs,order)=mv_madd([x],mv_shorten(xs,order),order);
neuper@37950: 
neuper@37950: (*. zeros a list .*)
neuper@37950: fun mv_null2([])=[]
neuper@37950:   | mv_null2(x::l)=0::mv_null2(l);
neuper@37950: 
neuper@37950: (*. multiplies two multivariate polynomials .*)
neuper@37950: fun mv_mul([]:mv_poly,[]:mv_poly,_) = []:mv_poly
neuper@37950:   | mv_mul([],y::p2,_) = [(0,mv_null2(#2(y)))]
neuper@37950:   | mv_mul(x::p1,[],_) = [(0,mv_null2(#2(x)))] 
neuper@37950:   | mv_mul(x::p1,y::p2,order) = mv_shorten(rev(sort (mv_geq order) (mv_mmul(x,y) :: (mv_mul(p1,y::p2,order) @
neuper@37950: 									    mv_mul([x],p2,order)))),order);
neuper@37950: 
neuper@37950: (*. gets the maximum value of a list .*)
neuper@37950: fun mv_getmax([])=0
neuper@37950:   | mv_getmax(x::p1)= let 
neuper@37950: 		       val m=mv_getmax(p1);
neuper@37950: 		   in
neuper@37950: 		       if m>x then m
neuper@37950: 		       else x
neuper@37950: 		   end;
neuper@37950: (*. calculates the maximum degree of an multivariate polynomial .*)
neuper@37950: fun mv_grad([]:mv_poly) = 0 
neuper@37950:   | mv_grad(p1:mv_poly)= mv_getmax((map mv_addlist) ((map #2) p1));
neuper@37950: 
neuper@37950: (*. converts the sign of a value .*)
neuper@37950: fun mv_minus(x)=(~1) * x;
neuper@37950: 
neuper@37950: (*. converts the sign of all coefficients of a polynomial .*)
neuper@37950: fun mv_minus2([]:mv_poly)=[]:mv_poly
neuper@37950:   | mv_minus2(p1)=(mv_minus(#1(hd(p1))),#2(hd(p1)))::(mv_minus2(tl(p1)));
neuper@37950: 
neuper@37950: (*. searches for a negativ value in a list .*)
neuper@37950: fun mv_is_negativ([])=false
neuper@37950:   | mv_is_negativ(x::xs)=if x<0 then true else mv_is_negativ(xs);
neuper@37950: 
neuper@37950: (*. division of monomials .*)
neuper@37950: fun mv_mdiv((0,[]):mv_monom,_:mv_monom)=(0,[]):mv_monom
neuper@37950:   | mv_mdiv(_,(0,[]))= raise error ("RATIONALS_MV_MDIV_EXCEPTION Division by 0 ")
neuper@37950:   | mv_mdiv(p1:mv_monom,p2:mv_monom)= 
neuper@37950:     let
neuper@37950: 	val c=ref (#1(p2));
neuper@37950: 	val pp=ref [];
neuper@37950:     in 
neuper@37950: 	(
neuper@37950: 	 if !c=0 then raise error("MV_MDIV_EXCEPTION Dividing by zero")
neuper@37950: 	 else c:=(#1(p1) div #1(p2));
neuper@37950: 	     if #1(p2)<>0 then 
neuper@37950: 		 (
neuper@37950: 		  pp:=(#2(mv_mmul((1,#2(p1)),(1,(map mv_minus) (#2(p2))))));
neuper@37950: 		  if mv_is_negativ(!pp) then (0,!pp)
neuper@37950: 		  else (!c,!pp) 
neuper@37950: 		      )
neuper@37950: 	     else raise error("MV_MDIV_EXCEPTION Dividing by empty Polynom")
neuper@37950: 		 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. prints a polynom for (internal use only) .*)
neuper@37978: fun mv_print_poly([]:mv_poly)=writeln("[]\n")
neuper@37978:   | mv_print_poly((x,y)::[])= writeln("("^Int.toString(x)^","^ints2str(y)^")\n")
neuper@37978:   | mv_print_poly((x,y)::p1) = (writeln("("^Int.toString(x)^","^ints2str(y)^"),");mv_print_poly(p1));
neuper@37950: 
neuper@37950: 
neuper@37950: (*. division of two multivariate polynomials .*) 
neuper@37950: fun mv_division([]:mv_poly,g:mv_poly,order)=([]:mv_poly,[]:mv_poly)
neuper@37950:   | mv_division(f,[],order)= raise error ("RATIONALS_MV_DIVISION_EXCEPTION Division by zero")
neuper@37950:   | mv_division(f,g,order)=
neuper@37950:     let 
neuper@37950: 	val r=ref [];
neuper@37950: 	val q=ref [];
neuper@37950: 	val g'=ref [];
neuper@37950: 	val k=ref 0;
neuper@37950: 	val m=ref (0,[0]);
neuper@37950: 	val exit=ref 0;
neuper@37950:     in
neuper@37950: 	r := rev(sort (mv_geq order) (mv_shorten(f,order)));
neuper@37950: 	g':= rev(sort (mv_geq order) (mv_shorten(g,order)));
neuper@37950: 	if #1(hd(!g'))=0 then raise error("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero") else ();
neuper@37950: 	if  (mv_monom_greater (hd(!g'),hd(!r),order)) then ([(0,mv_null2(#2(hd(f))))],(!r))
neuper@37950: 	else
neuper@37950: 	    (
neuper@37950: 	     exit:=0;
neuper@37950: 	     while (if (!exit)=0 then not(mv_monom_greater (hd(!g'),hd(!r),order)) else false) do
neuper@37950: 		 (
neuper@37950: 		  if (#1(mv_lm(!g',order)))<>0 then m:=mv_mdiv(mv_lm(!r,order),mv_lm(!g',order))
neuper@37950: 		  else raise error ("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero");	  
neuper@37950: 		  if #1(!m)<>0 then
neuper@37950: 		      ( 
neuper@37950: 		       q:=(!m)::(!q);
neuper@37950: 		       r:=mv_add((!r),mv_minus2(mv_mul(!g',[!m],order)),order)
neuper@37950: 		       )
neuper@37950: 		  else exit:=1;
neuper@37950: 		  if (if length(!r)<>0 then length(!g')<>0 else false) then ()
neuper@37950: 		  else (exit:=1)
neuper@37950: 		  );
neuper@37950: 		 (rev(!q),!r)
neuper@37950: 		 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. multiplies a polynomial with an integer .*)
neuper@37950: fun mv_skalar_mul([]:mv_poly,c) = []:mv_poly
neuper@37950:   | mv_skalar_mul((x,y)::p1,c) = ((x * c),y)::mv_skalar_mul(p1,c); 
neuper@37950: 
neuper@37950: (*. inserts the a first variable into an polynomial with exponent v .*)
neuper@37950: fun mv_correct([]:mv_poly,v:int)=[]:mv_poly
neuper@37950:   | mv_correct((x,y)::list,v:int)=(x,v::y)::mv_correct(list,v);
neuper@37950: 
neuper@37950: (*. multivariate case .*)
neuper@37950: 
neuper@37950: (*. decides if x is a factor of y .*)
neuper@37950: fun mv_divides([]:mv_poly,[]:mv_poly)=  raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
neuper@37950:   | mv_divides(x,[]) =  raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")
neuper@37950:   | mv_divides(x:mv_poly,y:mv_poly) = #2(mv_division(y,x,LEX_))=[];
neuper@37950: 
neuper@37950: (*. gets the maximum of a and b .*)
neuper@37950: fun mv_max(a,b) = if a>b then a else b;
neuper@37950: 
neuper@37950: (*. gets the maximum exponent of a mv polynomial in the lexicographic term order .*)
neuper@37950: fun mv_deg([]:mv_poly) = 0  
neuper@37950:   | mv_deg(p1)=
neuper@37950:     let
neuper@37950: 	val p1'=mv_shorten(p1,LEX_);
neuper@37950:     in
neuper@37950: 	if length(p1')=0 then 0 
neuper@37950: 	else mv_max(hd(#2(hd(p1'))),mv_deg(tl(p1')))
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. gets the maximum exponent of a mv polynomial in the reverse lexicographic term order .*)
neuper@37950: fun mv_deg2([]:mv_poly) = 0
neuper@37950:   | mv_deg2(p1)=
neuper@37950:     let
neuper@37950: 	val p1'=mv_shorten(p1,LEX_);
neuper@37950:     in
neuper@37950: 	if length(p1')=0 then 0 
neuper@37950: 	else mv_max(hd(rev(#2(hd(p1')))),mv_deg2(tl(p1')))
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. evaluates the mv polynomial at the value v of the main variable .*)
neuper@37950: fun mv_subs([]:mv_poly,v) = []:mv_poly
neuper@37950:   | mv_subs((c,e)::p1:mv_poly,v) = mv_skalar_mul(mv_cut([(c,e)]),power v (hd(e))) @ mv_subs(p1,v);
neuper@37950: 
neuper@37950: (*. calculates the content of a uv-polynomial in mv-representation .*)
neuper@37950: fun uv_content2([]:mv_poly) = 0
neuper@37950:   | uv_content2((c,e)::p1) = (gcd_int c (uv_content2(p1)));
neuper@37950: 
neuper@37950: (*. converts a uv-polynomial from mv-representation to  uv-representation .*)
neuper@37950: fun uv_to_list ([]:mv_poly)=[]:uv_poly
neuper@37950:   | uv_to_list ((c1,e1)::others) = 
neuper@37950:     let
neuper@37950: 	val count=ref 0;
neuper@37950: 	val max=mv_grad((c1,e1)::others); 
neuper@37950: 	val help=ref ((c1,e1)::others);
neuper@37950: 	val list=ref [];
neuper@37950:     in
neuper@37950: 	if length(e1)>1 then raise error ("RATIONALS_TO_LIST_EXCEPTION: not univariate")
neuper@37950: 	else if length(e1)=0 then [c1]
neuper@37950: 	     else
neuper@37950: 		 (
neuper@37950: 		  count:=0;
neuper@37950: 		  while (!count)<=max do
neuper@37950: 		      (
neuper@37950: 		       if length(!help)>0 andalso hd(#2(hd(!help)))=max-(!count) then 
neuper@37950: 			   (
neuper@37950: 			    list:=(#1(hd(!help)))::(!list);		       
neuper@37950: 			    help:=tl(!help) 
neuper@37950: 			    )
neuper@37950: 		       else 
neuper@37950: 			   (
neuper@37950: 			    list:= 0::(!list)
neuper@37950: 			    );
neuper@37950: 		       count := (!count) + 1
neuper@37950: 		       );
neuper@37950: 		      (!list)
neuper@37950: 		      )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. converts a uv-polynomial from uv-representation to mv-representation .*)
neuper@37950: fun uv_to_poly ([]:uv_poly) = []:mv_poly
neuper@37950:   | uv_to_poly p1 = 
neuper@37950:     let
neuper@37950: 	val count=ref 0;
neuper@37950: 	val help=ref p1;
neuper@37950: 	val list=ref [];
neuper@37950:     in
neuper@37950: 	while length(!help)>0 do
neuper@37950: 	    (
neuper@37950: 	     if hd(!help)=0 then ()
neuper@37950: 	     else list:=(hd(!help),[!count])::(!list);
neuper@37950: 	     count:=(!count)+1;
neuper@37950: 	     help:=tl(!help)
neuper@37950: 	     );
neuper@37950: 	    (!list)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. univariate gcd calculation from polynomials in multivariate representation .*)
neuper@37950: fun uv_gcd ([]:mv_poly) p2 = p2
neuper@37950:   | uv_gcd p1 ([]:mv_poly) = p1
neuper@37950:   | uv_gcd p1 [(c,[e])] = 
neuper@37950:     let 
neuper@37950: 	val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
neuper@37950: 	val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
neuper@37950:     in
neuper@37950: 	[(gcd_int (uv_content2(p1)) c,[min])]
neuper@37950:     end
neuper@37950:   | uv_gcd [(c,[e])] p2 = 
neuper@37950:     let 
neuper@37950: 	val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p2,LEX_))));
neuper@37950: 	val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));
neuper@37950:     in
neuper@37950: 	[(gcd_int (uv_content2(p2)) c,[min])]
neuper@37950:     end 
neuper@37950:   | uv_gcd p11 p22 = uv_to_poly(uv_mod_gcd (uv_to_list(mv_shorten(p11,LEX_))) (uv_to_list(mv_shorten(p22,LEX_))));
neuper@37950: 
neuper@37950: (*. help function for the newton interpolation .*)
neuper@37950: fun mv_newton_help ([]:mv_poly list,k:int) = []:mv_poly list
neuper@37950:   | mv_newton_help (pl:mv_poly list,k) = 
neuper@37950:     let
neuper@37950: 	val x=ref (rev(pl));
neuper@37950: 	val t=ref [];
neuper@37950: 	val y=ref [];
neuper@37950: 	val n=ref 1;
neuper@37950: 	val n1=ref[];
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 while length(!x)>1 do 
neuper@37950: 	     (
neuper@37950: 	      if length(hd(!x))>0 then n1:=mv_null2(#2(hd(hd(!x))))
neuper@37950: 	      else if length(hd(tl(!x)))>0 then n1:=mv_null2(#2(hd(hd(tl(!x)))))
neuper@37950: 		   else n1:=[]; 
neuper@37950: 	      t:= #1(mv_division(mv_add(hd(!x),mv_skalar_mul(hd(tl(!x)),~1),LEX_),[(k,!n1)],LEX_)); 
neuper@37950: 	      y:=(!t)::(!y);
neuper@37950: 	      x:=tl(!x)
neuper@37950: 	      );
neuper@37950: 	 (!y)
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950:     
neuper@37950: (*. help function for the newton interpolation .*)
neuper@37950: fun mv_newton_add ([]:mv_poly list) t = []:mv_poly
neuper@37950:   | mv_newton_add [x:mv_poly] t = x 
neuper@37950:   | mv_newton_add (pl:mv_poly list) t = 
neuper@37950:     let
neuper@37950: 	val expos=ref [];
neuper@37950: 	val pll=ref pl;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 
neuper@37950: 	 while length(!pll)>0 andalso hd(!pll)=[]  do 
neuper@37950: 	     ( 
neuper@37950: 	      pll:=tl(!pll)
neuper@37950: 	      ); 
neuper@37950: 	 if length(!pll)>0 then expos:= #2(hd(hd(!pll))) else expos:=[];
neuper@37950: 	 mv_add(hd(pl),
neuper@37950: 		mv_mul(
neuper@37950: 		       mv_add(mv_correct(mv_cut([(1,mv_null2(!expos))]),1),[(~t,mv_null2(!expos))],LEX_),
neuper@37950: 		       mv_newton_add (tl(pl)) (t+1),
neuper@37950: 		       LEX_
neuper@37950: 		       ),
neuper@37950: 		LEX_)
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the newton interpolation with polynomial coefficients .*)
neuper@37950: (*. step-depth is 1 and if the result is not an integerpolynomial .*)
neuper@37950: (*. this function returns [] .*)
neuper@37950: fun mv_newton ([]:(mv_poly) list) = []:mv_poly 
neuper@37950:   | mv_newton ([mp]:(mv_poly) list) = mp:mv_poly
neuper@37950:   | mv_newton pl =
neuper@37950:     let
neuper@37950: 	val c=ref pl;
neuper@37950: 	val c1=ref [];
neuper@37950: 	val n=length(pl);
neuper@37950: 	val k=ref 1;
neuper@37950: 	val i=ref n;
neuper@37950: 	val ppl=ref [];
neuper@37950:     in
neuper@37950: 	c1:=hd(pl)::[];
neuper@37950: 	c:=mv_newton_help(!c,!k);
neuper@37950: 	c1:=(hd(!c))::(!c1);
neuper@37950: 	while(length(!c)>1 andalso !k<n) do
neuper@37950: 	    (	 
neuper@37950: 	     k:=(!k)+1; 
neuper@37950: 	     while  length(!c)>0 andalso hd(!c)=[] do c:=tl(!c); 	  
neuper@37950: 	     if !c=[] then () else c:=mv_newton_help(!c,!k);
neuper@37950: 	     ppl:= !c;
neuper@37950: 	     if !c=[] then () else  c1:=(hd(!c))::(!c1)
neuper@37950: 	     );
neuper@37950: 	while  hd(!c1)=[] do c1:=tl(!c1);
neuper@37950: 	c1:=rev(!c1);
neuper@37950: 	ppl:= !c1;
neuper@37950: 	mv_newton_add (!c1) 1
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. sets the exponents of the first variable to zero .*)
neuper@37950: fun mv_null3([]:mv_poly)    = []:mv_poly
neuper@37950:   | mv_null3((x,y)::xs) = (x,0::tl(y))::mv_null3(xs);
neuper@37950: 
neuper@37950: (*. calculates the minimum exponents of a multivariate polynomial .*)
neuper@37950: fun mv_min_pp([]:mv_poly)=[]
neuper@37950:   | mv_min_pp((c,e)::xs)=
neuper@37950:     let
neuper@37950: 	val y=ref xs;
neuper@37950: 	val x=ref [];
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 x:=e;
neuper@37950: 	 while length(!y)>0 do
neuper@37950: 	     (
neuper@37950: 	      x:=(map uv_mod_min) ((!x) ~~ (#2(hd(!y))));
neuper@37950: 	      y:=tl(!y)
neuper@37950: 	      );
neuper@37950: 	 !x
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. checks if all elements of the list have value zero .*)
neuper@37950: fun list_is_null [] = true 
neuper@37950:   | list_is_null (x::xs) = (x=0 andalso list_is_null(xs)); 
neuper@37950: 
neuper@37950: (* check if main variable is zero*)
neuper@37950: fun main_zero (ms : mv_poly) = (list_is_null o (map (hd o #2))) ms;
neuper@37950: 
neuper@37950: (*. calculates the content of an polynomial .*)
neuper@37950: fun mv_content([]:mv_poly) = []:mv_poly
neuper@37950:   | mv_content(p1) = 
neuper@37950:     let
neuper@37950: 	val list=ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));
neuper@37950: 	val test=ref (hd(#2(hd(!list))));
neuper@37950: 	val result=ref []; 
neuper@37950: 	val min=(hd(#2(hd(rev(!list)))));
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 if length(!list)>1 then
neuper@37950: 	     (
neuper@37950: 	      while (if length(!list)>0 then (hd(#2(hd(!list)))=(!test)) else false) do
neuper@37950: 		  (
neuper@37950: 		   result:=(#1(hd(!list)),tl(#2(hd(!list))))::(!result);
neuper@37950: 		   
neuper@37950: 		   if length(!list)<1 then list:=[]
neuper@37950: 		   else list:=tl(!list) 
neuper@37950: 		       
neuper@37950: 		       );		  
neuper@37950: 		  if length(!list)>0 then  
neuper@37950: 		   ( 
neuper@37950: 		    list:=mv_gcd (!result) (mv_cut(mv_content(!list))) 
neuper@37950: 		    ) 
neuper@37950: 		  else list:=(!result); 
neuper@37950: 		  list:=mv_correct(!list,0);  
neuper@37950: 		  (!list) 
neuper@37950: 		  )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      mv_null3(!list) 
neuper@37950: 	      )
neuper@37950: 	     )
neuper@37950:     end
neuper@37950: 
neuper@37950: (*. calculates the primitiv part of a polynomial .*)
neuper@37950: and mv_pp([]:mv_poly) = []:mv_poly
neuper@37950:   | mv_pp(p1) = let
neuper@37950: 		    val cont=ref []; 
neuper@37950: 		    val pp=ref[];
neuper@37950: 		in
neuper@37950: 		    cont:=mv_content(p1);
neuper@37950: 		    pp:=(#1(mv_division(p1,!cont,LEX_)));
neuper@37950: 		    if !pp=[] 
neuper@37950: 			then raise error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")
neuper@37950: 		    else (!pp)
neuper@37950: 		end
neuper@37950: 
neuper@37950: (*. calculates the gcd of two multivariate polynomials with a modular approach .*)
neuper@37950: and mv_gcd ([]:mv_poly) ([]:mv_poly) :mv_poly= []:mv_poly
neuper@37950:   | mv_gcd ([]:mv_poly) (p2) :mv_poly= p2:mv_poly
neuper@37950:   | mv_gcd (p1:mv_poly) ([]) :mv_poly= p1:mv_poly
neuper@37950:   | mv_gcd ([(x,xs)]:mv_poly) ([(y,ys)]):mv_poly = 
neuper@37950:      let
neuper@37950:       val xpoly:mv_poly = [(x,xs)];
neuper@37950:       val ypoly:mv_poly = [(y,ys)];
neuper@37950:      in 
neuper@37950: 	(
neuper@37950: 	 if xs=ys then [((gcd_int x y),xs)]
neuper@37950: 	 else [((gcd_int x y),(map uv_mod_min)(xs~~ys))]:mv_poly
neuper@37950:         )
neuper@37950:     end 
neuper@37950:   | mv_gcd (p1:mv_poly) ([(y,ys)]) :mv_poly= 
neuper@37950: 	(
neuper@37950: 	 [(gcd_int (uv_content2(p1)) (y),(map  uv_mod_min)(mv_min_pp(p1)~~ys))]:mv_poly
neuper@37950: 	)
neuper@37950:   | mv_gcd ([(y,ys)]:mv_poly) (p2):mv_poly = 
neuper@37950: 	(
neuper@37950:          [(gcd_int (uv_content2(p2)) (y),(map  uv_mod_min)(mv_min_pp(p2)~~ys))]:mv_poly
neuper@37950:         )
neuper@37950:   | mv_gcd (p1':mv_poly) (p2':mv_poly):mv_poly=
neuper@37950:     let
neuper@37950: 	val vc=length(#2(hd(p1')));
neuper@37950: 	val cont = 
neuper@37950: 		  (
neuper@37950:                    if main_zero(mv_content(p1')) andalso 
neuper@37950:                      (main_zero(mv_content(p2'))) then
neuper@37950:                      mv_correct((mv_gcd (mv_cut(mv_content(p1'))) (mv_cut(mv_content(p2')))),0)
neuper@37950:                    else 
neuper@37950:                      mv_gcd (mv_content(p1')) (mv_content(p2'))
neuper@37950:                   );
neuper@37950: 	val p1= #1(mv_division(p1',mv_content(p1'),LEX_));
neuper@37950: 	val p2= #1(mv_division(p2',mv_content(p2'),LEX_)); 
neuper@37950: 	val gcd=ref [];
neuper@37950: 	val candidate=ref [];
neuper@37950: 	val interpolation_list=ref [];
neuper@37950: 	val delta=ref [];
neuper@37950:         val p1r = ref [];
neuper@37950:         val p2r = ref [];
neuper@37950:         val p1r' = ref [];
neuper@37950:         val p2r' = ref [];
neuper@37950: 	val factor=ref [];
neuper@37950: 	val r=ref 0;
neuper@37950: 	val gcd_r=ref [];
neuper@37950: 	val d=ref 0;
neuper@37950: 	val exit=ref 0;
neuper@37950: 	val current_degree=ref 99999; (*. FIXME: unlimited ! .*)
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 if vc<2 then (* areUnivariate(p1',p2') *)
neuper@37950: 	     (
neuper@37950: 	      gcd:=uv_gcd (mv_shorten(p1',LEX_)) (mv_shorten(p2',LEX_))
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      while !exit=0 do
neuper@37950: 		  (
neuper@37950: 		   r:=(!r)+1;
neuper@37950:                    p1r := mv_lc(p1,LEX_);
neuper@37950: 		   p2r := mv_lc(p2,LEX_);
neuper@37950:                    if main_zero(!p1r) andalso
neuper@37950:                       main_zero(!p2r) 
neuper@37950:                    then
neuper@37950:                        (
neuper@37950:                         delta := mv_correct((mv_gcd (mv_cut (!p1r)) (mv_cut (!p2r))),0)
neuper@37950:                        )
neuper@37950:                    else
neuper@37950:                        (
neuper@37950: 		        delta := mv_gcd (!p1r) (!p2r)
neuper@37950:                        );
neuper@37950: 		   (*if mv_shorten(mv_subs(!p1r,!r),LEX_)=[] andalso 
neuper@37950: 		      mv_shorten(mv_subs(!p2r,!r),LEX_)=[] *)
neuper@37950: 		   if mv_lc2(mv_shorten(mv_subs(!p1r,!r),LEX_),LEX_)=0 andalso 
neuper@37950: 		      mv_lc2(mv_shorten(mv_subs(!p2r,!r),LEX_),LEX_)=0 
neuper@37950:                    then 
neuper@37950:                        (
neuper@37950: 		       )
neuper@37950: 		   else 
neuper@37950: 		       (
neuper@37950: 			gcd_r:=mv_shorten(mv_gcd (mv_shorten(mv_subs(p1,!r),LEX_)) 
neuper@37950: 					         (mv_shorten(mv_subs(p2,!r),LEX_)) ,LEX_);
neuper@37950: 			gcd_r:= #1(mv_division(mv_mul(mv_correct(mv_subs(!delta,!r),0),!gcd_r,LEX_),
neuper@37950: 					       mv_correct(mv_lc(!gcd_r,LEX_),0),LEX_));
neuper@37950: 			d:=mv_deg2(!gcd_r); (* deg(gcd_r,z) *)
neuper@37950: 			if (!d < !current_degree) then 
neuper@37950: 			    (
neuper@37950: 			     current_degree:= !d;
neuper@37950: 			     interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
neuper@37950: 			     )
neuper@37950: 			else
neuper@37950: 			    (
neuper@37950: 			     if (!d = !current_degree) then
neuper@37950: 				 (
neuper@37950: 				  interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)
neuper@37950: 				  )
neuper@37950: 			     else () 
neuper@37950: 				 )
neuper@37950: 			    );
neuper@37950: 		      if length(!interpolation_list)> uv_mod_min(mv_deg(p1),mv_deg(p2)) then 
neuper@37950: 			  (
neuper@37950: 			   candidate := mv_newton(rev(!interpolation_list));
neuper@37950: 			   if !candidate=[] then ()
neuper@37950: 			   else
neuper@37950: 			       (
neuper@37950: 				candidate:=mv_pp(!candidate);
neuper@37950: 				if mv_divides(!candidate,p1) andalso mv_divides(!candidate,p2) then
neuper@37950: 				    (
neuper@37950: 				     gcd:= mv_mul(!candidate,cont,LEX_);
neuper@37950: 				     exit:=1
neuper@37950: 				     )
neuper@37950: 				else ()
neuper@37950: 				    );
neuper@37950: 			       interpolation_list:=[mv_correct(!gcd_r,0)]
neuper@37950: 			       )
neuper@37950: 		      else ()
neuper@37950: 			  )
neuper@37950: 	     );
neuper@37950: 	     (!gcd):mv_poly
neuper@37950: 	     )
neuper@37950:     end;	
neuper@37950: 
neuper@37950: 
neuper@37950: (*. calculates the least common divisor of two polynomials .*)
neuper@37950: fun mv_lcm (p1:mv_poly) (p2:mv_poly) :mv_poly = 
neuper@37950:     (
neuper@37950:      #1(mv_division(mv_mul(p1,p2,LEX_),mv_gcd p1 p2,LEX_))
neuper@37950:      );
neuper@37950: 
neuper@37950: (*. gets the variables (strings) of a term .*)
neuper@37950: fun get_vars(term1) = (map free2str) (vars term1); (*["a","b","c"]; *)
neuper@37950: 
neuper@37950: (*. counts the negative coefficents in a polynomial .*)
neuper@37950: fun count_neg ([]:mv_poly) = 0 
neuper@37950:   | count_neg ((c,e)::xs) = if c<0 then 1+count_neg xs
neuper@37950: 			  else count_neg xs;
neuper@37950: 
neuper@37950: (*. help function for is_polynomial  
neuper@37950:     checks the order of the operators .*)
neuper@37950: fun test_polynomial (Const ("uminus",_) $ Free (str,_)) _ = true (*WN.13.3.03*)
neuper@37950:   | test_polynomial (t as Free(str,_)) v = true
neuper@37950:   | test_polynomial (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false
neuper@37950: 						     else (test_polynomial t1 "*") andalso (test_polynomial t2 "*")
neuper@37950:   | test_polynomial (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false 
neuper@37950: 							  else (test_polynomial t1 " ") andalso (test_polynomial t2 " ")
neuper@37950:   | test_polynomial (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_polynomial t1 "^") andalso (test_polynomial t2 "^")
neuper@37950:   | test_polynomial _ v = false;  
neuper@37950: 
neuper@37950: (*. tests if a term is a polynomial .*)  
neuper@37950: fun is_polynomial t = test_polynomial t " ";
neuper@37950: 
neuper@37950: (*. help function for is_expanded 
neuper@37950:     checks the order of the operators .*)
neuper@37950: fun test_exp (t as Free(str,_)) v = true 
neuper@37950:   | test_exp (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false
neuper@37950: 						     else (test_exp t1 "*") andalso (test_exp t2 "*")
neuper@37950:   | test_exp (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false 
neuper@37950: 							  else (test_exp t1 " ") andalso (test_exp t2 " ") 
neuper@37950:   | test_exp (t as Const ("op -",_) $ t1 $ t2) v = if v="*" orelse v="^" then false 
neuper@37950: 							  else (test_exp t1 " ") andalso (test_exp t2 " ")
neuper@37950:   | test_exp (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_exp t1 "^") andalso (test_exp t2 "^")
neuper@37950:   | test_exp  _ v = false;
neuper@37950: 
neuper@37950: 
neuper@37950: (*. help function for check_coeff: 
neuper@37950:     converts the term to a list of coefficients .*) 
neuper@37950: fun term2coef' (t as Free(str,_(*typ*))) v :mv_poly option = 
neuper@37950:     let
neuper@37950: 	val x=ref NONE;
neuper@37950: 	val len=ref 0;
neuper@37950: 	val vl=ref [];
neuper@37950: 	val vh=ref [];
neuper@37950: 	val i=ref 0;
neuper@37950:     in 
neuper@37950: 	if is_numeral str then
neuper@37950: 	    (
neuper@37950: 	     SOME [(((the o int_of_str) str),mv_null2(v))] handle _ => NONE
neuper@37950: 		 )
neuper@37950: 	else (* variable *)
neuper@37950: 	    (
neuper@37950: 	     len:=length(v);
neuper@37950: 	     vh:=v;
neuper@37950: 	     while ((!len)>(!i)) do
neuper@37950: 		 (
neuper@37950: 		  if str=hd((!vh)) then
neuper@37950: 		      (
neuper@37950: 		       vl:=1::(!vl)
neuper@37950: 		       )
neuper@37950: 		  else 
neuper@37950: 		      (
neuper@37950: 		       vl:=0::(!vl)
neuper@37950: 		       );
neuper@37950: 		      vh:=tl(!vh);
neuper@37950: 		      i:=(!i)+1    
neuper@37950: 		      );		
neuper@37950: 		 SOME [(1,rev(!vl))] handle _ => NONE
neuper@37950: 	    )
neuper@37950:     end
neuper@37950:   | term2coef' (Const ("op *",_) $ t1 $ t2) v :mv_poly option= 
neuper@37950:     let
neuper@37950: 	val t1pp=ref [];
neuper@37950: 	val t2pp=ref [];
neuper@37950: 	val t1c=ref 0;
neuper@37950: 	val t2c=ref 0;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 t1pp:=(#2(hd(the(term2coef' t1 v))));
neuper@37950: 	 t2pp:=(#2(hd(the(term2coef' t2 v))));
neuper@37950: 	 t1c:=(#1(hd(the(term2coef' t1 v))));
neuper@37950: 	 t2c:=(#1(hd(the(term2coef' t2 v))));
neuper@37950: 	
neuper@37950: 	 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] handle _ => NONE 
neuper@37950: 		
neuper@37950: 	 )
neuper@37950:     end
neuper@37950:   | term2coef' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ (t2 as Free (str2,_))) v :mv_poly option= 
neuper@37950:     let
neuper@37950: 	val x=ref NONE;
neuper@37950: 	val len=ref 0;
neuper@37950: 	val vl=ref [];
neuper@37950: 	val vh=ref [];
neuper@37950: 	val vtemp=ref [];
neuper@37950: 	val i=ref 0;	 
neuper@37950:     in
neuper@37950:     (
neuper@37950:      if (not o is_numeral) str1 andalso is_numeral str2 then
neuper@37950: 	 (
neuper@37950: 	  len:=length(v);
neuper@37950: 	  vh:=v;
neuper@37950: 
neuper@37950: 	  while ((!len)>(!i)) do
neuper@37950: 	      (
neuper@37950: 	       if str1=hd((!vh)) then
neuper@37950: 		   (
neuper@37950: 		    vl:=((the o int_of_str) str2)::(!vl)
neuper@37950: 		    )
neuper@37950: 	       else 
neuper@37950: 		   (
neuper@37950: 		    vl:=0::(!vl)
neuper@37950: 		    );
neuper@37950: 		   vh:=tl(!vh);
neuper@37950: 		   i:=(!i)+1     
neuper@37950: 		   );
neuper@37950: 	      SOME [(1,rev(!vl))] handle _ => NONE
neuper@37950: 	      )
neuper@37950:      else raise error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term")
neuper@37950: 	 )
neuper@37950:     end
neuper@37950:   | term2coef' (Const ("op +",_) $ t1 $ t2) v :mv_poly option= 
neuper@37950:     (
neuper@37950:      SOME ((the(term2coef' t1 v)) @ (the(term2coef' t2 v))) handle _ => NONE
neuper@37950: 	 )
neuper@37950:   | term2coef' (Const ("op -",_) $ t1 $ t2) v :mv_poly option= 
neuper@37950:     (
neuper@37950:      SOME ((the(term2coef' t1 v)) @ mv_skalar_mul((the(term2coef' t2 v)),1)) handle _ => NONE
neuper@37950: 	 )
neuper@37950:   | term2coef' (term) v = raise error ("RATIONALS_TERM2COEF_EXCEPTION 2: Invalid term");
neuper@37950: 
neuper@37950: (*. checks if all coefficients of a polynomial are positiv (except the first) .*)
neuper@37950: fun check_coeff t = (* erste Koeffizient kann <0 sein !!! *)
neuper@37950:     if count_neg(tl(the(term2coef' t (get_vars(t)))))=0 then true 
neuper@37950:     else false;
neuper@37950: 
neuper@37950: (*. checks for expanded term [3] .*)
neuper@37950: fun is_expanded t = test_exp t " " andalso check_coeff(t); 
neuper@37950: 
neuper@37950: (*WN.7.3.03 Hilfsfunktion f"ur term2poly'*)
neuper@37950: fun mk_monom v' p vs = 
neuper@37950:     let fun conv p (v: string) = if v'= v then p else 0
neuper@37950:     in map (conv p) vs end;
neuper@37950: (* mk_monom "y" 5 ["a","b","x","y","z"];
neuper@37950: val it = [0,0,0,5,0] : int list*)
neuper@37950: 
neuper@37950: (*. this function converts the term representation into the internal representation mv_poly .*)
neuper@37950: fun term2poly' (Const ("uminus",_) $ Free (str,_)) v = (*WN.7.3.03*)
neuper@37950:     if is_numeral str 
neuper@37950:     then SOME [((the o int_of_str) ("-"^str), mk_monom "#" 0 v)]
neuper@37950:     else SOME [(~1, mk_monom str 1 v)]
neuper@37950: 
neuper@37950:   | term2poly' (Free(str,_)) v :mv_poly option = 
neuper@37950:     let
neuper@37950: 	val x=ref NONE;
neuper@37950: 	val len=ref 0;
neuper@37950: 	val vl=ref [];
neuper@37950: 	val vh=ref [];
neuper@37950: 	val i=ref 0;
neuper@37950:     in 
neuper@37950: 	if is_numeral str then
neuper@37950: 	    (
neuper@37950: 	     SOME [(((the o int_of_str) str),mv_null2 v)] handle _ => NONE
neuper@37950: 		 )
neuper@37950: 	else (* variable *)
neuper@37950: 	    (
neuper@37950: 	     len:=length v;
neuper@37950: 	     vh:= v;
neuper@37950: 	     while ((!len)>(!i)) do
neuper@37950: 		 (
neuper@37950: 		  if str=hd((!vh)) then
neuper@37950: 		      (
neuper@37950: 		       vl:=1::(!vl)
neuper@37950: 		       )
neuper@37950: 		  else 
neuper@37950: 		      (
neuper@37950: 		       vl:=0::(!vl)
neuper@37950: 		       );
neuper@37950: 		      vh:=tl(!vh);
neuper@37950: 		      i:=(!i)+1    
neuper@37950: 		      );		
neuper@37950: 		 SOME [(1,rev(!vl))] handle _ => NONE
neuper@37950: 	    )
neuper@37950:     end
neuper@37950:   | term2poly' (Const ("op *",_) $ t1 $ t2) v :mv_poly option= 
neuper@37950:     let
neuper@37950: 	val t1pp=ref [];
neuper@37950: 	val t2pp=ref [];
neuper@37950: 	val t1c=ref 0;
neuper@37950: 	val t2c=ref 0;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 t1pp:=(#2(hd(the(term2poly' t1 v))));
neuper@37950: 	 t2pp:=(#2(hd(the(term2poly' t2 v))));
neuper@37950: 	 t1c:=(#1(hd(the(term2poly' t1 v))));
neuper@37950: 	 t2c:=(#1(hd(the(term2poly' t2 v))));
neuper@37950: 	
neuper@37950: 	 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] 
neuper@37950: 	 handle _ => NONE 
neuper@37950: 		
neuper@37950: 	 )
neuper@37950:     end
neuper@37950:   | term2poly' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ 
neuper@37950: 		      (t2 as Free (str2,_))) v :mv_poly option= 
neuper@37950:     let
neuper@37950: 	val x=ref NONE;
neuper@37950: 	val len=ref 0;
neuper@37950: 	val vl=ref [];
neuper@37950: 	val vh=ref [];
neuper@37950: 	val vtemp=ref [];
neuper@37950: 	val i=ref 0;	 
neuper@37950:     in
neuper@37950:     (
neuper@37950:      if (not o is_numeral) str1 andalso is_numeral str2 then
neuper@37950: 	 (
neuper@37950: 	  len:=length(v);
neuper@37950: 	  vh:=v;
neuper@37950: 
neuper@37950: 	  while ((!len)>(!i)) do
neuper@37950: 	      (
neuper@37950: 	       if str1=hd((!vh)) then
neuper@37950: 		   (
neuper@37950: 		    vl:=((the o int_of_str) str2)::(!vl)
neuper@37950: 		    )
neuper@37950: 	       else 
neuper@37950: 		   (
neuper@37950: 		    vl:=0::(!vl)
neuper@37950: 		    );
neuper@37950: 		   vh:=tl(!vh);
neuper@37950: 		   i:=(!i)+1     
neuper@37950: 		   );
neuper@37950: 	      SOME [(1,rev(!vl))] handle _ => NONE
neuper@37950: 	      )
neuper@37950:      else raise error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term")
neuper@37950: 	 )
neuper@37950:     end
neuper@37950:   | term2poly' (Const ("op +",_) $ t1 $ t2) v :mv_poly option = 
neuper@37950:     (
neuper@37950:      SOME ((the(term2poly' t1 v)) @ (the(term2poly' t2 v))) handle _ => NONE
neuper@37950: 	 )
neuper@37950:   | term2poly' (Const ("op -",_) $ t1 $ t2) v :mv_poly option = 
neuper@37950:     (
neuper@37950:      SOME ((the(term2poly' t1 v)) @ mv_skalar_mul((the(term2poly' t2 v)),~1)) handle _ => NONE
neuper@37950: 	 )
neuper@37950:   | term2poly' (term) v = raise error ("RATIONALS_TERM2POLY_EXCEPTION 2: Invalid term");
neuper@37950: 
neuper@37950: (*. translates an Isabelle term into internal representation.
neuper@37950:     term2poly
neuper@37950:     fn : term ->              (*normalform [2]                    *)
neuper@37950:     	 string list ->       (*for ...!!! BITTE DIE ERKLÄRUNG, 
neuper@37950:     			       DIE DU MIR LETZTES MAL GEGEBEN HAST*)
neuper@37950:     	 mv_monom list        (*internal representation           *)
neuper@37950:     		  option      (*the translation may fail with NONE*)
neuper@37950: .*)
neuper@37950: fun term2poly (t:term) v = 
neuper@37950:      if is_polynomial t then term2poly' t v
neuper@37950:      else raise error ("term2poly: invalid = "^(term2str t));
neuper@37950: 
neuper@37950: (*. same as term2poly with automatic detection of the variables .*)
neuper@37950: fun term2polyx t = term2poly t (((map free2str) o vars) t); 
neuper@37950: 
neuper@37950: (*. checks if the term is in expanded polynomial form and converts it into the internal representation .*)
neuper@37950: fun expanded2poly (t:term) v = 
neuper@37950:     (*if is_expanded t then*) term2poly' t v
neuper@37950:     (*else raise error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*);
neuper@37950: 
neuper@37950: (*. same as expanded2poly with automatic detection of the variables .*)
neuper@37950: fun expanded2polyx t = expanded2poly t (((map free2str) o vars) t);
neuper@37950: 
neuper@37950: (*. converts a powerproduct into term representation .*)
neuper@37950: fun powerproduct2term(xs,v) =  
neuper@37950:     let
neuper@37950: 	val xss=ref xs;
neuper@37950: 	val vv=ref v;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 while hd(!xss)=0 do 
neuper@37950: 	     (
neuper@37950: 	      xss:=tl(!xss);
neuper@37950: 	      vv:=tl(!vv)
neuper@37950: 	      );
neuper@37950: 	     
neuper@37950: 	 if list_is_null(tl(!xss)) then 
neuper@37950: 	     (
neuper@37950: 	      if hd(!xss)=1 then Free(hd(!vv), HOLogic.realT)
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   Free(hd(!vv), HOLogic.realT) $  Free(str_of_int (hd(!xss)),HOLogic.realT)
neuper@37950: 		   )
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      if hd(!xss)=1 then 
neuper@37950: 		  ( 
neuper@37950: 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   Free(hd(!vv), HOLogic.realT) $
neuper@37950: 		   powerproduct2term(tl(!xss),tl(!vv))
neuper@37950: 		   )
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   (
neuper@37950: 		    Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		    Free(hd(!vv), HOLogic.realT) $  Free(str_of_int (hd(!xss)),HOLogic.realT)
neuper@37950: 		    ) $
neuper@37950: 		    powerproduct2term(tl(!xss),tl(!vv))
neuper@37950: 		   )
neuper@37950: 	      )
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. converts a monom into term representation .*)
neuper@37950: (*fun monom2term ((c,e):mv_monom, v:string list) = 
neuper@37950:     if c=0 then Free(str_of_int 0,HOLogic.realT)  
neuper@37950:     else
neuper@37950: 	(
neuper@37950: 	 if list_is_null(e) then
neuper@37950: 	     ( 
neuper@37950: 	      Free(str_of_int c,HOLogic.realT)  
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      if c=1 then 
neuper@37950: 		  (
neuper@37950: 		   powerproduct2term(e,v)
neuper@37950: 		   )
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
neuper@37950: 		   Free(str_of_int c,HOLogic.realT)  $
neuper@37950: 		   powerproduct2term(e,v)
neuper@37950: 		   )
neuper@37950: 		  )
neuper@37950: 	     );*)
neuper@37950: 
neuper@37950: 
neuper@37950: (*fun monom2term ((i, is):mv_monom, v) = 
neuper@37950:     if list_is_null is 
neuper@37950:     then 
neuper@37950: 	if i >= 0 
neuper@37950: 	then Free (str_of_int i, HOLogic.realT)
neuper@37950: 	else Const ("uminus", HOLogic.realT --> HOLogic.realT) $
neuper@37950: 		   Free ((str_of_int o abs) i, HOLogic.realT)
neuper@37950:     else
neuper@37950: 	if i > 0 
neuper@37950: 	then Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
neuper@37950: 		   (Free (str_of_int i, HOLogic.realT)) $
neuper@37950: 		   powerproduct2term(is, v)
neuper@37950: 	else Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $
neuper@37950: 		   (Const ("uminus", HOLogic.realT --> HOLogic.realT) $
neuper@37950: 		   Free ((str_of_int o abs) i, HOLogic.realT)) $
neuper@37950: 		   powerproduct2term(is, vs);---------------------------*)
neuper@37950: fun monom2term ((i, is) : mv_monom, vs) = 
neuper@37950:     if list_is_null is 
neuper@37950:     then Free (str_of_int i, HOLogic.realT)
neuper@37950:     else if i = 1
neuper@37950:     then powerproduct2term (is, vs)
neuper@37950:     else Const ("op *", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
neuper@37950: 	       (Free (str_of_int i, HOLogic.realT)) $
neuper@37950: 	       powerproduct2term (is, vs);
neuper@37950:     
neuper@37950: (*. converts the internal polynomial representation into an Isabelle term.*)
neuper@37950: fun poly2term' ([] : mv_poly, vs) = Free(str_of_int 0, HOLogic.realT)  
neuper@37950:   | poly2term' ([(c, e) : mv_monom], vs) = monom2term ((c, e), vs)
neuper@37950:   | poly2term' ((c, e) :: ces, vs) =  
neuper@37950:     Const("op +", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $
neuper@37950:          poly2term (ces, vs) $ monom2term ((c, e), vs)
neuper@37950: and poly2term (xs, vs) = poly2term' (rev (sort (mv_geq LEX_) (xs)), vs);
neuper@37950: 
neuper@37950: 
neuper@37950: (*. converts a monom into term representation .*)
neuper@37950: (*. ignores the sign of the coefficients => use only for exp-poly functions .*)
neuper@37950: fun monom2term2((c,e):mv_monom, v:string list) =  
neuper@37950:     if c=0 then Free(str_of_int 0,HOLogic.realT)  
neuper@37950:     else
neuper@37950: 	(
neuper@37950: 	 if list_is_null(e) then
neuper@37950: 	     ( 
neuper@37950: 	      Free(str_of_int (abs(c)),HOLogic.realT)  
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      if abs(c)=1 then 
neuper@37950: 		  (
neuper@37950: 		   powerproduct2term(e,v)
neuper@37950: 		   )
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
neuper@37950: 		   Free(str_of_int (abs(c)),HOLogic.realT)  $
neuper@37950: 		   powerproduct2term(e,v)
neuper@37950: 		   )
neuper@37950: 		  )
neuper@37950: 	     );
neuper@37950: 
neuper@37950: (*. converts the expanded polynomial representation into the term representation .*)
neuper@37950: fun exp2term' ([]:mv_poly,vars) =  Free(str_of_int 0,HOLogic.realT)  
neuper@37950:   | exp2term' ([(c,e)],vars) =     monom2term((c,e),vars) 			     
neuper@37950:   | exp2term' ((c1,e1)::others,vars) =  
neuper@37950:     if c1<0 then 	
neuper@37950: 	Const("op -",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
neuper@37950: 	exp2term'(others,vars) $
neuper@37950: 	( 
neuper@37950: 	 monom2term2((c1,e1),vars)
neuper@37950: 	 ) 
neuper@37950:     else
neuper@37950: 	Const("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $
neuper@37950: 	exp2term'(others,vars) $
neuper@37950: 	( 
neuper@37950: 	 monom2term2((c1,e1),vars)
neuper@37950: 	 );
neuper@37950: 	
neuper@37950: (*. sorts the powerproduct by lexicographic termorder and converts them into 
neuper@37950:     a term in polynomial representation .*)
neuper@37950: fun poly2expanded (xs,vars) = exp2term'(rev(sort (mv_geq LEX_) (xs)),vars);
neuper@37950: 
neuper@37950: (*. converts a polynomial into expanded form .*)
neuper@37950: fun polynomial2expanded t =  
neuper@37950:     (let 
neuper@37950: 	val vars=(((map free2str) o vars) t);
neuper@37950:     in
neuper@37950: 	SOME (poly2expanded (the (term2poly t vars), vars))
neuper@37950:     end) handle _ => NONE;
neuper@37950: 
neuper@37950: (*. converts a polynomial into polynomial form .*)
neuper@37950: fun expanded2polynomial t =  
neuper@37950:     (let 
neuper@37950: 	val vars=(((map free2str) o vars) t);
neuper@37950:     in
neuper@37950: 	SOME (poly2term (the (expanded2poly t vars), vars))
neuper@37950:     end) handle _ => NONE;
neuper@37950: 
neuper@37950: 
neuper@37950: (*. calculates the greatest common divisor of numerator and denominator and seperates it from each .*)
neuper@37950: fun step_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) = 
neuper@37950:     let
neuper@37950: 	val p1' = ref [];
neuper@37950: 	val p2' = ref [];
neuper@37950: 	val p3  = ref []
neuper@37950: 	val vars = rev(get_vars(p1) union get_vars(p2));
neuper@37950:     in
neuper@37950: 	(
neuper@37950:          p1':= sort (mv_geq LEX_) (the (term2poly p1 vars ));
neuper@37950:        	 p2':= sort (mv_geq LEX_) (the (term2poly p2 vars ));
neuper@37950: 	 p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
neuper@37950: 	 if (!p3)=[(1,mv_null2(vars))] then 
neuper@37950: 	     (
neuper@37950: 	      Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 
neuper@37950: 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
neuper@37950: 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
neuper@37950: 	      
neuper@37950: 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then
neuper@37950: 	      (
neuper@37950: 	       Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 	       $ 
neuper@37950: 	       (
neuper@37950: 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		poly2term(!p1',vars) $ 
neuper@37950: 		poly2term(!p3,vars) 
neuper@37950: 		) 
neuper@37950: 	       $ 
neuper@37950: 	       (
neuper@37950: 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		poly2term(!p2',vars) $ 
neuper@37950: 		poly2term(!p3,vars)
neuper@37950: 		) 	
neuper@37950: 	       )	
neuper@37950: 	      else
neuper@37950: 	      (
neuper@37950: 	       p1':=mv_skalar_mul(!p1',~1);
neuper@37950: 	       p2':=mv_skalar_mul(!p2',~1);
neuper@37950: 	       p3:=mv_skalar_mul(!p3,~1);
neuper@37950: 	       (
neuper@37950: 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2term(!p1',vars) $ 
neuper@37950: 		 poly2term(!p3,vars) 
neuper@37950: 		 ) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2term(!p2',vars) $ 
neuper@37950: 		 poly2term(!p3,vars)
neuper@37950: 		 ) 	
neuper@37950: 		)	
neuper@37950: 	       )	  
neuper@37950: 	      )
neuper@37950: 	     )
neuper@37950:     end
neuper@37950: | step_cancel _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction"); 
neuper@37950: 
neuper@37950: 
neuper@37950: (*. same as step_cancel, this time for expanded forms (input+output) .*)
neuper@37950: fun step_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) = 
neuper@37950:     let
neuper@37950: 	val p1' = ref [];
neuper@37950: 	val p2' = ref [];
neuper@37950: 	val p3  = ref []
neuper@37950: 	val vars = rev(get_vars(p1) union get_vars(p2));
neuper@37950:     in
neuper@37950: 	(
neuper@37950:          p1':= sort (mv_geq LEX_) (the (expanded2poly p1 vars ));
neuper@37950:        	 p2':= sort (mv_geq LEX_) (the (expanded2poly p2 vars ));
neuper@37950: 	 p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
neuper@37950: 	 if (!p3)=[(1,mv_null2(vars))] then 
neuper@37950: 	     (
neuper@37950: 	      Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 
neuper@37950: 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
neuper@37950: 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
neuper@37950: 	      
neuper@37950: 	      if #1(hd(sort (mv_geq LEX_) (!p2')))(* mv_lc2(!p2',LEX_)*)>0 then
neuper@37950: 	      (
neuper@37950: 	       Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 	       $ 
neuper@37950: 	       (
neuper@37950: 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		poly2expanded(!p1',vars) $ 
neuper@37950: 		poly2expanded(!p3,vars) 
neuper@37950: 		) 
neuper@37950: 	       $ 
neuper@37950: 	       (
neuper@37950: 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		poly2expanded(!p2',vars) $ 
neuper@37950: 		poly2expanded(!p3,vars)
neuper@37950: 		) 	
neuper@37950: 	       )	
neuper@37950: 	      else
neuper@37950: 	      (
neuper@37950: 	       p1':=mv_skalar_mul(!p1',~1);
neuper@37950: 	       p2':=mv_skalar_mul(!p2',~1);
neuper@37950: 	       p3:=mv_skalar_mul(!p3,~1);
neuper@37950: 	       (
neuper@37950: 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2expanded(!p1',vars) $ 
neuper@37950: 		 poly2expanded(!p3,vars) 
neuper@37950: 		 ) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2expanded(!p2',vars) $ 
neuper@37950: 		 poly2expanded(!p3,vars)
neuper@37950: 		 ) 	
neuper@37950: 		)	
neuper@37950: 	       )	  
neuper@37950: 	      )
neuper@37950: 	     )
neuper@37950:     end
neuper@37950: | step_cancel_expanded _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction"); 
neuper@37950: 
neuper@37950: (*. calculates the greatest common divisor of numerator and denominator and divides each through it .*)
neuper@37950: fun direct_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) = 
neuper@37950:     let
neuper@37950: 	val p1' = ref [];
neuper@37950: 	val p2' = ref [];
neuper@37950: 	val p3  = ref []
neuper@37950: 	val vars = rev(get_vars(p1) union get_vars(p2));
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 p1':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p1 vars )),LEX_));
neuper@37950: 	 p2':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p2 vars )),LEX_));	 
neuper@37950: 	 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
neuper@37950: 
neuper@37950: 	 if (!p3)=[(1,mv_null2(vars))] then 
neuper@37950: 	     (
neuper@37950: 	      (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
neuper@37950: 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
neuper@37950: 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then	      
neuper@37950: 	      (
neuper@37950: 	       (
neuper@37950: 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 poly2term((!p1'),vars)
neuper@37950: 		 ) 
neuper@37950: 		$ 
neuper@37950: 		( 
neuper@37950: 		 poly2term((!p2'),vars)
neuper@37950: 		 ) 	
neuper@37950: 		)
neuper@37950: 	       ,
neuper@37950: 	       if mv_grad(!p3)>0 then 
neuper@37950: 		   [
neuper@37950: 		    (
neuper@37950: 		     Const ("Not",[bool]--->bool) $
neuper@37950: 		     (
neuper@37950: 		      Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
neuper@37950: 		      poly2term((!p3),vars) $
neuper@37950: 		      Free("0",HOLogic.realT)
neuper@37950: 		      )
neuper@37950: 		     )
neuper@37950: 		    ]
neuper@37950: 	       else
neuper@37950: 		   []
neuper@37950: 		   )
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   p1':=mv_skalar_mul(!p1',~1);
neuper@37950: 		   p2':=mv_skalar_mul(!p2',~1);
neuper@37950: 		   if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1); 
neuper@37950: 		       (
neuper@37950: 			Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 			$ 
neuper@37950: 			(
neuper@37950: 			 poly2term((!p1'),vars)
neuper@37950: 			 ) 
neuper@37950: 			$ 
neuper@37950: 			( 
neuper@37950: 			 poly2term((!p2'),vars)
neuper@37950: 			 ) 	
neuper@37950: 			,
neuper@37950: 			if mv_grad(!p3)>0 then 
neuper@37950: 			    [
neuper@37950: 			     (
neuper@37950: 			      Const ("Not",[bool]--->bool) $
neuper@37950: 			      (
neuper@37950: 			       Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
neuper@37950: 			       poly2term((!p3),vars) $
neuper@37950: 			       Free("0",HOLogic.realT)
neuper@37950: 			       )
neuper@37950: 			      )
neuper@37950: 			     ]
neuper@37950: 			else
neuper@37950: 			    []
neuper@37950: 			    )
neuper@37950: 		       )
neuper@37950: 		  )
neuper@37950: 	     )
neuper@37950:     end
neuper@37950:   | direct_cancel _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction"); 
neuper@37950: 
neuper@37950: (*. same es direct_cancel, this time for expanded forms (input+output).*) 
neuper@37950: fun direct_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) =  
neuper@37950:     let
neuper@37950: 	val p1' = ref [];
neuper@37950: 	val p2' = ref [];
neuper@37950: 	val p3  = ref []
neuper@37950: 	val vars = rev(get_vars(p1) union get_vars(p2));
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_));
neuper@37950: 	 p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_));	 
neuper@37950: 	 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));
neuper@37950: 
neuper@37950: 	 if (!p3)=[(1,mv_null2(vars))] then 
neuper@37950: 	     (
neuper@37950: 	      (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])
neuper@37950: 	      )
neuper@37950: 	 else
neuper@37950: 	     (
neuper@37950: 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));
neuper@37950: 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));
neuper@37950: 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then	      
neuper@37950: 	      (
neuper@37950: 	       (
neuper@37950: 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		$ 
neuper@37950: 		(
neuper@37950: 		 poly2expanded((!p1'),vars)
neuper@37950: 		 ) 
neuper@37950: 		$ 
neuper@37950: 		( 
neuper@37950: 		 poly2expanded((!p2'),vars)
neuper@37950: 		 ) 	
neuper@37950: 		)
neuper@37950: 	       ,
neuper@37950: 	       if mv_grad(!p3)>0 then 
neuper@37950: 		   [
neuper@37950: 		    (
neuper@37950: 		     Const ("Not",[bool]--->bool) $
neuper@37950: 		     (
neuper@37950: 		      Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
neuper@37950: 		      poly2expanded((!p3),vars) $
neuper@37950: 		      Free("0",HOLogic.realT)
neuper@37950: 		      )
neuper@37950: 		     )
neuper@37950: 		    ]
neuper@37950: 	       else
neuper@37950: 		   []
neuper@37950: 		   )
neuper@37950: 	      else
neuper@37950: 		  (
neuper@37950: 		   p1':=mv_skalar_mul(!p1',~1);
neuper@37950: 		   p2':=mv_skalar_mul(!p2',~1);
neuper@37950: 		   if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1); 
neuper@37950: 		       (
neuper@37950: 			Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 			$ 
neuper@37950: 			(
neuper@37950: 			 poly2expanded((!p1'),vars)
neuper@37950: 			 ) 
neuper@37950: 			$ 
neuper@37950: 			( 
neuper@37950: 			 poly2expanded((!p2'),vars)
neuper@37950: 			 ) 	
neuper@37950: 			,
neuper@37950: 			if mv_grad(!p3)>0 then 
neuper@37950: 			    [
neuper@37950: 			     (
neuper@37950: 			      Const ("Not",[bool]--->bool) $
neuper@37950: 			      (
neuper@37950: 			       Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $
neuper@37950: 			       poly2expanded((!p3),vars) $
neuper@37950: 			       Free("0",HOLogic.realT)
neuper@37950: 			       )
neuper@37950: 			      )
neuper@37950: 			     ]
neuper@37950: 			else
neuper@37950: 			    []
neuper@37950: 			    )
neuper@37950: 		       )
neuper@37950: 		  )
neuper@37950: 	     )
neuper@37950:     end
neuper@37950:   | direct_cancel_expanded _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction"); 
neuper@37950: 
neuper@37950: 
neuper@37950: (*. adds two fractions .*)
neuper@37950: fun add_fract ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =
neuper@37950:     let
neuper@37950: 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
neuper@37950: 	val t11'=ref (the(term2poly t11 vars));
neuper@37950: val _= writeln"### add_fract: done t11"
neuper@37950: 	val t12'=ref (the(term2poly t12 vars));
neuper@37950: val _= writeln"### add_fract: done t12"
neuper@37950: 	val t21'=ref (the(term2poly t21 vars));
neuper@37950: val _= writeln"### add_fract: done t21"
neuper@37950: 	val t22'=ref (the(term2poly t22 vars));
neuper@37950: val _= writeln"### add_fract: done t22"
neuper@37950: 	val den=ref [];
neuper@37950: 	val nom=ref [];
neuper@37950: 	val m1=ref [];
neuper@37950: 	val m2=ref [];
neuper@37950:     in
neuper@37950: 	
neuper@37950: 	(
neuper@37950: 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
neuper@37950: writeln"### add_fract: done sort mv_lcm";
neuper@37950: 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
neuper@37950: writeln"### add_fract: done sort mv_division t12";
neuper@37950: 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
neuper@37950: writeln"### add_fract: done sort mv_division t22";
neuper@37950: 	 nom :=sort (mv_geq LEX_) 
neuper@37950: 		    (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),
neuper@37950: 				       mv_mul(!t21',!m2,LEX_),
neuper@37950: 				       LEX_),
neuper@37950: 				LEX_));
neuper@37950: writeln"### add_fract: done sort mv_add";
neuper@37950: 	 (
neuper@37950: 	  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 	  $ 
neuper@37950: 	  (
neuper@37950: 	   poly2term((!nom),vars)
neuper@37950: 	   ) 
neuper@37950: 	  $ 
neuper@37950: 	  ( 
neuper@37950: 	   poly2term((!den),vars)
neuper@37950: 	   )	      
neuper@37950: 	  )
neuper@37950: 	 )	     
neuper@37950:     end 
neuper@37950:   | add_fract (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call");
neuper@37950: 
neuper@37950: (*. adds two expanded fractions .*)
neuper@37950: fun add_fract_exp ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =
neuper@37950:     let
neuper@37950: 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);
neuper@37950: 	val t11'=ref (the(expanded2poly t11 vars));
neuper@37950: 	val t12'=ref (the(expanded2poly t12 vars));
neuper@37950: 	val t21'=ref (the(expanded2poly t21 vars));
neuper@37950: 	val t22'=ref (the(expanded2poly t22 vars));
neuper@37950: 	val den=ref [];
neuper@37950: 	val nom=ref [];
neuper@37950: 	val m1=ref [];
neuper@37950: 	val m2=ref [];
neuper@37950:     in
neuper@37950: 	
neuper@37950: 	(
neuper@37950: 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));
neuper@37950: 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));
neuper@37950: 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));
neuper@37950: 	 nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_));
neuper@37950: 	 (
neuper@37950: 	  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 	  $ 
neuper@37950: 	  (
neuper@37950: 	   poly2expanded((!nom),vars)
neuper@37950: 	   ) 
neuper@37950: 	  $ 
neuper@37950: 	  ( 
neuper@37950: 	   poly2expanded((!den),vars)
neuper@37950: 	   )	      
neuper@37950: 	  )
neuper@37950: 	 )	     
neuper@37950:     end 
neuper@37950:   | add_fract_exp (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call");
neuper@37950: 
neuper@37950: (*. adds a list of terms .*)
neuper@37950: fun add_list_of_fractions []= (Free("0",HOLogic.realT),[])
neuper@37950:   | add_list_of_fractions [x]= direct_cancel x
neuper@37950:   | add_list_of_fractions (x::y::xs) = 
neuper@37950:     let
neuper@37950: 	val (t1a,rest1)=direct_cancel(x);
neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(x)";
neuper@37950: 	val (t2a,rest2)=direct_cancel(y);
neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(y)";
neuper@37950: 	val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs));
neuper@37950: val _= writeln"### add_list_of_fractions xs: has done add_list_of_fraction xs";
neuper@37950: 	val (t4a,rest4)=direct_cancel(t3a);
neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(t3a)";
neuper@37950: 	val rest=rest1 union rest2 union rest3 union rest4;
neuper@37950:     in
neuper@37950: 	(writeln"### add_list_of_fractions in";
neuper@37950: 	 (
neuper@37950: 	 (t4a,rest) 
neuper@37950: 	 )
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. adds a list of expanded terms .*)
neuper@37950: fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])
neuper@37950:   | add_list_of_fractions_exp [x]= direct_cancel_expanded x
neuper@37950:   | add_list_of_fractions_exp (x::y::xs) = 
neuper@37950:     let
neuper@37950: 	val (t1a,rest1)=direct_cancel_expanded(x);
neuper@37950: 	val (t2a,rest2)=direct_cancel_expanded(y);
neuper@37950: 	val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));
neuper@37950: 	val (t4a,rest4)=direct_cancel_expanded(t3a);
neuper@37950: 	val rest=rest1 union rest2 union rest3 union rest4;
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 (t4a,rest) 
neuper@37950: 	 )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the lcm of a list of mv_poly .*)
neuper@37950: fun calc_lcm ([x],var)= (x,var) 
neuper@37950:   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
neuper@37950: 
neuper@37950: (*. converts a list of terms to a list of mv_poly .*)
neuper@37950: fun t2d([],_)=[] 
neuper@37950:   | t2d((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars); 
neuper@37950: 
neuper@37950: (*. same as t2d, this time for expanded forms .*)
neuper@37950: fun t2d_exp([],_)=[]  
neuper@37950:   | t2d_exp((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
neuper@37950: 
neuper@37950: (*. converts a list of fract terms to a list of their denominators .*)
neuper@37950: fun termlist2denominators [] = ([],[])
neuper@37950:   | termlist2denominators xs = 
neuper@37950:     let	
neuper@37950: 	val xxs=ref xs;
neuper@37950: 	val var=ref [];
neuper@37950:     in
neuper@37950: 	var:=[];
neuper@37950: 	while length(!xxs)>0 do
neuper@37950: 	    (
neuper@37950: 	     let 
neuper@37950: 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
neuper@37950: 	     in
neuper@37950: 		 (
neuper@37950: 		  xxs:=tl(!xxs);
neuper@37950: 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
neuper@37950: 		  )
neuper@37950: 	     end
neuper@37950: 	     );
neuper@37950: 	    (t2d(xs,!var),!var)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the lcm of a list of mv_poly .*)
neuper@37950: fun calc_lcm ([x],var)= (x,var) 
neuper@37950:   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);
neuper@37950: 
neuper@37950: (*. converts a list of terms to a list of mv_poly .*)
neuper@37950: fun t2d([],_)=[] 
neuper@37950:   | t2d((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars); 
neuper@37950: 
neuper@37950: (*. same as t2d, this time for expanded forms .*)
neuper@37950: fun t2d_exp([],_)=[]  
neuper@37950:   | t2d_exp((t as (Const("HOL.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);
neuper@37950: 
neuper@37950: (*. converts a list of fract terms to a list of their denominators .*)
neuper@37950: fun termlist2denominators [] = ([],[])
neuper@37950:   | termlist2denominators xs = 
neuper@37950:     let	
neuper@37950: 	val xxs=ref xs;
neuper@37950: 	val var=ref [];
neuper@37950:     in
neuper@37950: 	var:=[];
neuper@37950: 	while length(!xxs)>0 do
neuper@37950: 	    (
neuper@37950: 	     let 
neuper@37950: 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
neuper@37950: 	     in
neuper@37950: 		 (
neuper@37950: 		  xxs:=tl(!xxs);
neuper@37950: 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
neuper@37950: 		  )
neuper@37950: 	     end
neuper@37950: 	     );
neuper@37950: 	    (t2d(xs,!var),!var)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. same as termlist2denminators, this time for expanded forms .*)
neuper@37950: fun termlist2denominators_exp [] = ([],[])
neuper@37950:   | termlist2denominators_exp xs = 
neuper@37950:     let	
neuper@37950: 	val xxs=ref xs;
neuper@37950: 	val var=ref [];
neuper@37950:     in
neuper@37950: 	var:=[];
neuper@37950: 	while length(!xxs)>0 do
neuper@37950: 	    (
neuper@37950: 	     let 
neuper@37950: 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);
neuper@37950: 	     in
neuper@37950: 		 (
neuper@37950: 		  xxs:=tl(!xxs);
neuper@37950: 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))
neuper@37950: 		  )
neuper@37950: 	     end
neuper@37950: 	     );
neuper@37950: 	    (t2d_exp(xs,!var),!var)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. reduces all fractions to the least common denominator .*)
neuper@37950: fun com_den(x::xs,denom,den,var)=
neuper@37950:     let 
neuper@37950: 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
neuper@37950: 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
neuper@37950: 	val p3= #1(mv_division(denom,p2,LEX_));
neuper@37950: 	val p1var=get_vars(p1');
neuper@37950:     in     
neuper@37950: 	if length(xs)>0 then 
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
neuper@37950: 		 $ 
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  poly2term(the (term2poly p1' p1var),p1var)
neuper@37950: 		  $ 
neuper@37950: 		  den	
neuper@37950: 		  )    
neuper@37950: 		 $ 
neuper@37950: 		 #1(com_den(xs,denom,den,var))
neuper@37950: 		,
neuper@37950: 		[]
neuper@37950: 		)
neuper@37950: 	    else
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		 $ 
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   poly2term(the (term2poly p1' p1var),p1var) $ 
neuper@37950: 		   poly2term(p3,var)
neuper@37950: 		   ) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   den
neuper@37950: 		   ) 	
neuper@37950: 		  )
neuper@37950: 		 $ 
neuper@37950: 		 #1(com_den(xs,denom,den,var))
neuper@37950: 		,
neuper@37950: 		[]
neuper@37950: 		)
neuper@37950: 	else
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  poly2term(the (term2poly p1' p1var),p1var)
neuper@37950: 		  $ 
neuper@37950: 		  den	
neuper@37950: 		  )
neuper@37950: 		 ,
neuper@37950: 		 []
neuper@37950: 		 )
neuper@37950: 	     else
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   poly2term(the (term2poly p1' p1var),p1var) $ 
neuper@37950: 		   poly2term(p3,var)
neuper@37950: 		   ) 
neuper@37950: 		  $ 
neuper@37950: 		  den 	
neuper@37950: 		  ,
neuper@37950: 		  []
neuper@37950: 		  )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. same as com_den, this time for expanded forms .*)
neuper@37950: fun com_den_exp(x::xs,denom,den,var)=
neuper@37950:     let 
neuper@37950: 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
neuper@37950: 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
neuper@37950: 	val p3= #1(mv_division(denom,p2,LEX_));
neuper@37950: 	val p1var=get_vars(p1');
neuper@37950:     in     
neuper@37950: 	if length(xs)>0 then 
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)
neuper@37950: 		 $ 
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  poly2expanded(the(expanded2poly p1' p1var),p1var)
neuper@37950: 		  $ 
neuper@37950: 		  den	
neuper@37950: 		  )    
neuper@37950: 		 $ 
neuper@37950: 		 #1(com_den_exp(xs,denom,den,var))
neuper@37950: 		,
neuper@37950: 		[]
neuper@37950: 		)
neuper@37950: 	    else
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		 $ 
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $ 
neuper@37950: 		   poly2expanded(p3,var)
neuper@37950: 		   ) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   den
neuper@37950: 		   ) 	
neuper@37950: 		  )
neuper@37950: 		 $ 
neuper@37950: 		 #1(com_den_exp(xs,denom,den,var))
neuper@37950: 		,
neuper@37950: 		[]
neuper@37950: 		)
neuper@37950: 	else
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  poly2expanded(the(expanded2poly p1' p1var),p1var)
neuper@37950: 		  $ 
neuper@37950: 		  den	
neuper@37950: 		  )
neuper@37950: 		 ,
neuper@37950: 		 []
neuper@37950: 		 )
neuper@37950: 	     else
neuper@37950: 		 (
neuper@37950: 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) 
neuper@37950: 		  $ 
neuper@37950: 		  (
neuper@37950: 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $ 
neuper@37950: 		   poly2expanded(p3,var)
neuper@37950: 		   ) 
neuper@37950: 		  $ 
neuper@37950: 		  den 	
neuper@37950: 		  ,
neuper@37950: 		  []
neuper@37950: 		  )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (* wird aktuell nicht mehr gebraucht, bei rückänderung schon 
neuper@37950: -------------------------------------------------------------
neuper@37950: (* WN0210???SK brauch ma des überhaupt *)
neuper@37950: fun com_den2(x::xs,denom,den,var)=
neuper@37950:     let 
neuper@37950: 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
neuper@37950: 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));
neuper@37950: 	val p3= #1(mv_division(denom,p2,LEX_));
neuper@37950: 	val p1var=get_vars(p1');
neuper@37950:     in     
neuper@37950: 	if length(xs)>0 then 
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2term(the(term2poly p1' p1var),p1var) $ 
neuper@37950: 		 com_den2(xs,denom,den,var)
neuper@37950: 		)
neuper@37950: 	    else
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 (
neuper@37950: 		   let 
neuper@37950: 		       val p3'=poly2term(p3,var);
neuper@37950: 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
neuper@37950: 		   in
neuper@37950: 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
neuper@37950: 		   end
neuper@37950: 		  ) $ 
neuper@37950: 		 com_den2(xs,denom,den,var)
neuper@37950: 		)
neuper@37950: 	else
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 poly2term(the(term2poly p1' p1var),p1var)
neuper@37950: 		 )
neuper@37950: 	     else
neuper@37950: 		 (
neuper@37950: 		   let 
neuper@37950: 		       val p3'=poly2term(p3,var);
neuper@37950: 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
neuper@37950: 		   in
neuper@37950: 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)
neuper@37950: 		   end
neuper@37950: 		  )
neuper@37950:     end;
neuper@37950: 
neuper@37950: (* WN0210???SK brauch ma des überhaupt *)
neuper@37950: fun com_den_exp2(x::xs,denom,den,var)=
neuper@37950:     let 
neuper@37950: 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;
neuper@37950: 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));
neuper@37950: 	val p3= #1(mv_division(denom,p2,LEX_));
neuper@37950: 	val p1var=get_vars p1';
neuper@37950:     in     
neuper@37950: 	if length(xs)>0 then 
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 poly2expanded(the (expanded2poly p1' p1var),p1var) $ 
neuper@37950: 		 com_den_exp2(xs,denom,den,var)
neuper@37950: 		)
neuper@37950: 	    else
neuper@37950: 		(
neuper@37950: 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 		 (
neuper@37950: 		   let 
neuper@37950: 		       val p3'=poly2expanded(p3,var);
neuper@37950: 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
neuper@37950: 		   in
neuper@37950: 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
neuper@37950: 		   end
neuper@37950: 		  ) $ 
neuper@37950: 		 com_den_exp2(xs,denom,den,var)
neuper@37950: 		)
neuper@37950: 	else
neuper@37950: 	    if p3=[(1,mv_null2(var))] then
neuper@37950: 		(
neuper@37950: 		 poly2expanded(the (expanded2poly p1' p1var),p1var)
neuper@37950: 		 )
neuper@37950: 	     else
neuper@37950: 		 (
neuper@37950: 		   let 
neuper@37950: 		       val p3'=poly2expanded(p3,var);
neuper@37950: 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');
neuper@37950: 		   in
neuper@37950: 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)
neuper@37950: 		   end
neuper@37950: 		  )
neuper@37950:     end;
neuper@37950: ---------------------------------------------------------*)
neuper@37950: 
neuper@37950: 
neuper@37950: (*. searches for an element y of a list ys, which has an gcd not 1 with x .*) 
neuper@37950: fun exists_gcd (x,[]) = false 
neuper@37950:   | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then  exists_gcd (x,ys)
neuper@37950: 			   else true;
neuper@37950: 
neuper@37950: (*. divides each element of the list xs with y .*)
neuper@37950: fun list_div ([],y) = [] 
neuper@37950:   | list_div (x::xs,y) = 
neuper@37950:     let
neuper@37950: 	val (d,r)=mv_division(x,y,LEX_);
neuper@37950:     in
neuper@37950: 	if r=[] then 
neuper@37950: 	    d::list_div(xs,y)
neuper@37950: 	else x::list_div(xs,y)
neuper@37950:     end;
neuper@37950:     
neuper@37950: (*. checks if x is in the list ys .*)
neuper@37950: fun in_list (x,[]) = false 
neuper@37950:   | in_list (x,y::ys) = if x=y then true
neuper@37950: 			else in_list(x,ys);
neuper@37950: 
neuper@37950: (*. deletes all equal elements of the list xs .*)
neuper@37950: fun kill_equal [] = [] 
neuper@37950:   | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)
neuper@37950: 			 else x::kill_equal(xs);
neuper@37950: 
neuper@37950: (*. searches for new factors .*)
neuper@37950: fun new_factors [] = []
neuper@37950:   | new_factors (list:mv_poly list):mv_poly list = 
neuper@37950:     let
neuper@37950: 	val l = kill_equal list;
neuper@37950: 	val len = length(l);
neuper@37950:     in
neuper@37950: 	if len>=2 then
neuper@37950: 	    (
neuper@37950: 	     let
neuper@37950: 		 val x::y::xs=l;
neuper@37950: 		 val gcd=mv_gcd x y;
neuper@37950: 	     in
neuper@37950: 		 if gcd=[(1,mv_null2(#2(hd(x))))] then 
neuper@37950: 		     ( 
neuper@37950: 		      if exists_gcd(x,xs) then new_factors (y::xs @ [x])
neuper@37950: 		      else x::new_factors(y::xs)
neuper@37950: 	             )
neuper@37950: 		 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))
neuper@37950: 	     end
neuper@37950: 	     )
neuper@37950: 	else
neuper@37950: 	    if len=1 then [hd(l)]
neuper@37950: 	    else []
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. gets the factors of a list .*)
neuper@37950: fun get_factors x = new_factors x; 
neuper@37950: 
neuper@37950: (*. multiplies the elements of the list .*)
neuper@37950: fun multi_list [] = []
neuper@37950:   | multi_list (x::xs) = if xs=[] then x
neuper@37950: 			 else mv_mul(x,multi_list xs,LEX_);
neuper@37950: 
neuper@37950: (*. makes a term out of the elements of the list (polynomial representation) .*)
neuper@37950: fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT) 
neuper@37950:   | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)
neuper@37950: 			       else
neuper@37950: 				   (
neuper@37950: 				    Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 				    poly2term(sort (mv_geq LEX_) (x),vars) $ 
neuper@37950: 				    make_term(xs,vars)
neuper@37950: 				    );
neuper@37950: 
neuper@37950: (*. factorizes the denominator (polynomial representation) .*)				
neuper@37950: fun factorize_den (l,den,vars) = 
neuper@37950:     let
neuper@37950: 	val factor_list=kill_equal( (get_factors l));
neuper@37950: 	val mlist=multi_list(factor_list);
neuper@37950: 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
neuper@37950:     in
neuper@37950: 	if rest=[] then
neuper@37950: 	    (
neuper@37950: 	     if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)
neuper@37950: 	     else make_term(last::factor_list,vars)
neuper@37950: 	     )
neuper@37950: 	else raise error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")
neuper@37950:     end; 
neuper@37950: 
neuper@37950: (*. makes a term out of the elements of the list (expanded polynomial representation) .*)
neuper@37950: fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT) 
neuper@37950:   | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)
neuper@37950: 			       else
neuper@37950: 				   (
neuper@37950: 				    Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 				    poly2expanded(sort (mv_geq LEX_) (x),vars) $ 
neuper@37950: 				    make_exp(xs,vars)
neuper@37950: 				    );
neuper@37950: 
neuper@37950: (*. factorizes the denominator (expanded polynomial representation) .*)	
neuper@37950: fun factorize_den_exp (l,den,vars) = 
neuper@37950:     let
neuper@37950: 	val factor_list=kill_equal( (get_factors l));
neuper@37950: 	val mlist=multi_list(factor_list);
neuper@37950: 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);
neuper@37950:     in
neuper@37950: 	if rest=[] then
neuper@37950: 	    (
neuper@37950: 	     if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)
neuper@37950: 	     else make_exp(last::factor_list,vars)
neuper@37950: 	     )
neuper@37950: 	else raise error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")
neuper@37950:     end; 
neuper@37950: 
neuper@37950: (*. calculates the common denominator of all elements of the list and multiplies .*)
neuper@37950: (*. the nominators and denominators with the correct factor .*)
neuper@37950: (*. (polynomial representation) .*)
neuper@37950: fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)
neuper@37950:   | step_add_list_of_fractions [x]= raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")
neuper@37950:   | step_add_list_of_fractions (xs) = 
neuper@37950:     let
neuper@37950:         val den_list=termlist2denominators (xs); (* list of denominators *)
neuper@37950: 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)
neuper@37950: 	val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
neuper@37950:     in
neuper@37950: 	com_den(xs,denom,den,var)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (*. calculates the common denominator of all elements of the list and multiplies .*)
neuper@37950: (*. the nominators and denominators with the correct factor .*)
neuper@37950: (*. (expanded polynomial representation) .*)
neuper@37950: fun step_add_list_of_fractions_exp []  = (Free("0",HOLogic.realT),[]:term list)
neuper@37950:   | step_add_list_of_fractions_exp [x] = raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")
neuper@37950:   | step_add_list_of_fractions_exp (xs)= 
neuper@37950:     let
neuper@37950:         val den_list=termlist2denominators_exp (xs); (* list of denominators *)
neuper@37950: 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)
neuper@37950: 	val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)
neuper@37950:     in
neuper@37950: 	com_den_exp(xs,denom,den,var)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (* wird aktuell nicht mehr gebraucht, bei rückänderung schon 
neuper@37950: -------------------------------------------------------------
neuper@37950: (* WN0210???SK brauch ma des überhaupt *)
neuper@37950: fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)
neuper@37950:   | step_add_list_of_fractions2 [x]=(x,[])
neuper@37950:   | step_add_list_of_fractions2 (xs) = 
neuper@37950:     let
neuper@37950:         val den_list=termlist2denominators (xs); (* list of denominators *)
neuper@37950: 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)
neuper@37950: 	val den=factorize_den(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 	 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $
neuper@37950: 	 poly2term(denom,var)
neuper@37950: 	,
neuper@37950: 	[]
neuper@37950: 	)
neuper@37950:     end;
neuper@37950: 
neuper@37950: (* WN0210???SK brauch ma des überhaupt *)
neuper@37950: fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)
neuper@37950:   | step_add_list_of_fractions2_exp [x]=(x,[])
neuper@37950:   | step_add_list_of_fractions2_exp (xs) = 
neuper@37950:     let
neuper@37950:         val den_list=termlist2denominators_exp (xs); (* list of denominators *)
neuper@37950: 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)
neuper@37950: 	val den=factorize_den_exp(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)
neuper@37950:     in
neuper@37950: 	(
neuper@37950: 	 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 	 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $
neuper@37950: 	 poly2expanded(denom,var)
neuper@37950: 	,
neuper@37950: 	[]
neuper@37950: 	)
neuper@37950:     end;
neuper@37950: ---------------------------------------------- *)
neuper@37950: 
neuper@37950: 
neuper@37950: (*. converts a term, which contains severel terms seperated by +, into a list of these terms .*)
neuper@37950: fun term2list (t as (Const("HOL.divide",_) $ _ $ _)) = [t]
neuper@37950:   | term2list (t as (Const("Atools.pow",_) $ _ $ _)) = 
neuper@37950:     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 	  t $ Free("1",HOLogic.realT)
neuper@37950:      ]
neuper@37950:   | term2list (t as (Free(_,_))) = 
neuper@37950:     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 	  t $  Free("1",HOLogic.realT)
neuper@37950:      ]
neuper@37950:   | term2list (t as (Const("op *",_) $ _ $ _)) = 
neuper@37950:     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ 
neuper@37950: 	  t $ Free("1",HOLogic.realT)
neuper@37950:      ]
neuper@37950:   | term2list (Const("op +",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)
neuper@37950:   | term2list (Const("op -",_) $ t1 $ t2) = 
neuper@37950:     raise error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")
neuper@37950:   | term2list _ = raise error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");
neuper@37950: 
neuper@37950: (*.factors out the gcd of nominator and denominator:
neuper@37950:    a/b = (a' * gcd)/(b' * gcd),  a,b,gcd  are poly[2].*)
neuper@37950: fun factout_p_  (thy:theory) t = SOME (step_cancel t,[]:term list); 
neuper@37950: fun factout_ (thy:theory) t = SOME (step_cancel_expanded t,[]:term list); 
neuper@37950: 
neuper@37950: (*.cancels a single fraction with normalform [2]
neuper@37950:    resulting in a canceled fraction [2], see factout_ .*)
neuper@37950: fun cancel_p_ (thy:theory) t = (*WN.2.6.03 no rewrite -> NONE !*)
neuper@37950:     (let val (t',asm) = direct_cancel(*_expanded ... corrected MG.21.8.03*) t
neuper@37950:      in if t = t' then NONE else SOME (t',asm) 
neuper@37950:      end) handle _ => NONE;
neuper@37950: (*.the same as above with normalform [3]
neuper@37950:   val cancel_ :
neuper@37950:       theory ->        (*10.02 unused                                    *)
neuper@37950:       term -> 	       (*fraction in normalform [3]                      *)
neuper@37950:       (term * 	       (*fraction in normalform [3]                      *)
neuper@37950:        term list)      (*casual asumptions in normalform [3]             *)
neuper@37950: 	  option       (*NONE: the function is not applicable            *).*)
neuper@37950: fun cancel_ (thy:theory) t = SOME (direct_cancel_expanded t) handle _ => NONE;
neuper@37950: 
neuper@37950: (*.transforms sums of at least 2 fractions [3] to
neuper@37950:    sums with the least common multiple as nominator.*)
neuper@37950: fun common_nominator_p_ (thy:theory) t =
neuper@37950: ((*writeln("### common_nominator_p_ called");*)
neuper@37950:     SOME (step_add_list_of_fractions(term2list(t))) handle _ => NONE
neuper@37950: );
neuper@37950: fun common_nominator_ (thy:theory) t =
neuper@37950:     SOME (step_add_list_of_fractions_exp(term2list(t))) handle _ => NONE;
neuper@37950: 
neuper@37950: (*.add 2 or more fractions
neuper@37950: val add_fraction_p_ :
neuper@37950:       theory ->        (*10.02 unused                                    *)
neuper@37950:       term -> 	       (*2 or more fractions with normalform [2]         *)
neuper@37950:       (term * 	       (*one fraction with normalform [2]                *)
neuper@37950:        term list)      (*casual assumptions in normalform [2] WN0210???SK  *)
neuper@37950: 	  option       (*NONE: the function is not applicable            *).*)
neuper@37950: fun add_fraction_p_ (thy:theory) t = 
neuper@37950: (writeln("### add_fraction_p_ called");
neuper@37950:     (let val ts = term2list t
neuper@37950:      in if 1 < length ts
neuper@37950: 	then SOME (add_list_of_fractions ts)
neuper@37950: 	else NONE (*raise error ("RATIONALS_ADD_EXCEPTION: nothing to add")*)
neuper@37950:      end) handle _ => NONE
neuper@37950: );
neuper@37950: (*.same as add_fraction_p_ but with normalform [3].*)
neuper@37950: (*SOME (step_add_list_of_fractions2(term2list(t))); *)
neuper@37950: fun add_fraction_ (thy:theory) t = 
neuper@37950:     if length(term2list(t))>1 
neuper@37950:     then SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE
neuper@37950:     else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)
neuper@37950: 	NONE;
neuper@37950: fun add_fraction_ (thy:theory) t = 
neuper@37950:     (if 1 < length (term2list t)
neuper@37950:      then SOME (add_list_of_fractions_exp (term2list t))
neuper@37950:      else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)
neuper@37950: 	 NONE) handle _ => NONE;
neuper@37950: 
neuper@37950: (*SOME (step_add_list_of_fractions2_exp(term2list(t))); *)
neuper@37950: 
neuper@37950: (*. brings the term into a normal form .*)
neuper@37950: fun norm_rational_ (thy:theory) t = 
neuper@37950:     SOME (add_list_of_fractions(term2list(t))) handle _ => NONE; 
neuper@37950: fun norm_expanded_rat_ (thy:theory) t = 
neuper@37950:     SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE; 
neuper@37950: 
neuper@37950: 
neuper@37950: (*.evaluates conditions in calculate_Rational.*)
neuper@37950: (*make local with FIXX@ME result:term *term list*)
neuper@37950: val calc_rat_erls = prep_rls(
neuper@37950:   Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950: 	 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [], *)
neuper@37950: 	 rules = 
neuper@37950: 	 [Calc ("op =",eval_equal "#equal_"),
neuper@37950: 	  Calc ("Atools.is'_const",eval_const "#is_const_"),
neuper@37978: 	  Thm ("not_true",num_str @{thm not_true}),
neuper@37978: 	  Thm ("not_false",num_str @{thm not_false})
neuper@37950: 	  ], 
neuper@37950: 	 scr = EmptyScr});
neuper@37950: 
neuper@37950: 
neuper@37950: (*.simplifies expressions with numerals;
neuper@37950:    does NOT rearrange the term by AC-rewriting; thus terms with variables 
neuper@37950:    need to have constants to be commuted together respectively.*)
neuper@37950: val calculate_Rational = prep_rls(
neuper@37950:     merge_rls "calculate_Rational"
neuper@37950: 	(Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950: 	      erls = calc_rat_erls, srls = Erls, (*asm_thm = [],*) 
neuper@37950: 	      calc = [], 
neuper@37950: 	      rules = 
neuper@37950: 	      [Calc ("HOL.divide"  ,eval_cancel "#divide_"),
neuper@37950: 	       
neuper@37965: 	       Thm ("minus_divide_left",
neuper@37978: 		    num_str (@{thm minus_divide_left} RS @{thm sym})),
neuper@37950: 	       (*SYM - ?x / ?y = - (?x / ?y)  may come from subst*)
neuper@37950: 	       
neuper@37969: 	       Thm ("rat_add",num_str @{thm rat_add}),
neuper@37950: 	       (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \
neuper@37950: 		 \"a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)
neuper@37969: 	       Thm ("rat_add1",num_str @{thm rat_add1}),
neuper@37950: 	       (*"[| a is_const; b is_const; c is_const |] ==> \
neuper@37950: 		 \"a / c + b / c = (a + b) / c"*)
neuper@37969: 	       Thm ("rat_add2",num_str @{thm rat_add2}),
neuper@37950: 	       (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> \
neuper@37950: 		 \?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)
neuper@37969: 	       Thm ("rat_add3",num_str @{thm rat_add3}),
neuper@37950: 	       (*"[| a is_const; b is_const; c is_const |] ==> \
neuper@37950: 		 \"a + b / c = (a * c) / c + b / c"\
neuper@37950: 		 \.... is_const to be omitted here FIXME*)
neuper@37950: 	       
neuper@37969: 	       Thm ("rat_mult",num_str @{thm rat_mult}),
neuper@37950: 	       (*a / b * (c / d) = a * c / (b * d)*)
neuper@37965: 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),
neuper@37950: 	       (*?x * (?y / ?z) = ?x * ?y / ?z*)
neuper@37965: 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),
neuper@37950: 	       (*?y / ?z * ?x = ?y * ?x / ?z*)
neuper@37950: 	       
neuper@37969: 	       Thm ("real_divide_divide1",num_str @{thm real_divide_divide1}),
neuper@37950: 	       (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)
neuper@37965: 	       Thm ("divide_divide_eq_left",num_str @{thm divide_divide_eq_left}),
neuper@37950: 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
neuper@37950: 	       
neuper@37969: 	       Thm ("rat_power", num_str @{thm rat_power}),
neuper@37950: 	       (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
neuper@37950: 	       
neuper@37969: 	       Thm ("mult_cross",num_str @{thm mult_cross}),
neuper@37950: 	       (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)
neuper@37969: 	       Thm ("mult_cross1",num_str @{thm mult_cross1}),
neuper@37950: 	       (*"   b ~= 0            ==> (a / b = c    ) = (a     = b * c)*)
neuper@37969: 	       Thm ("mult_cross2",num_str @{thm mult_cross2})
neuper@37950: 	       (*"           d ~= 0    ==> (a     = c / d) = (a * d =     c)*)
neuper@37950: 	       ], scr = EmptyScr})
neuper@37950: 	calculate_Poly);
neuper@37978: val ------------------------------------------------------+++ = "33333";
neuper@37950: 
neuper@37950: 
neuper@37950: (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)
neuper@37950: fun eval_is_expanded (thmid:string) _ 
neuper@37950: 		       (t as (Const("Rational.is'_expanded", _) $ arg)) thy = 
neuper@37950:     if is_expanded arg
neuper@37950:     then SOME (mk_thmid thmid "" 
neuper@37950: 			((Syntax.string_of_term (thy2ctxt thy)) arg) "", 
neuper@37950: 	       Trueprop $ (mk_equality (t, HOLogic.true_const)))
neuper@37950:     else SOME (mk_thmid thmid "" 
neuper@37950: 			((Syntax.string_of_term (thy2ctxt thy)) arg) "", 
neuper@37950: 	       Trueprop $ (mk_equality (t, HOLogic.false_const)))
neuper@37950:   | eval_is_expanded _ _ _ _ = NONE; 
neuper@37950: 
neuper@37950: val rational_erls = 
neuper@37950:     merge_rls "rational_erls" calculate_Rational 
neuper@37950: 	      (append_rls "is_expanded" Atools_erls 
neuper@37950: 			  [Calc ("Rational.is'_expanded", eval_is_expanded "")
neuper@37950: 			   ]);
neuper@37950: 
neuper@37950: 
neuper@37950: (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:
neuper@37950:  =================================================================
neuper@37950:  A[2] 'cancel_p': .
neuper@37950:  A[3] 'cancel': .
neuper@37950:  B[2] 'common_nominator_p': transforms summands in a term [2]
neuper@37950:          to fractions with the (least) common multiple as nominator.
neuper@37950:  B[3] 'norm_rational': normalizes arbitrary algebraic terms (without 
neuper@37950:          radicals and transzendental functions) to one canceled fraction,
neuper@37950: 	 nominator and denominator in polynomial form.
neuper@37950: 
neuper@37950: In order to meet isac's requirements for interactive and stepwise calculation,
neuper@37950: each 'reverse-rewerite-set' consists of an initialization for the interpreter 
neuper@37950: state and of 4 functions, each of which employs rewriting as much as possible.
neuper@37950: The signature of these functions are the same in each 'reverse-rewrite-set' 
neuper@37950: respectively.*)
neuper@37950: 
neuper@37950: (* ************************************************************************* *)
neuper@37950: 
neuper@37978: val ------------------------------------------------------++++ = "44444";
neuper@37950: 
neuper@37950: local(*. cancel_p
neuper@37950: ------------------------
neuper@37950: cancels a single fraction consisting of two (uni- or multivariate)
neuper@37950: polynomials WN0609???SK[2] into another such a fraction; examples:
neuper@37950: 
neuper@37950: 	   a^2 + -1*b^2         a + b
neuper@37950:         -------------------- = ---------
neuper@37950: 	a^2 + -2*a*b + b^2     a + -1*b
neuper@37950: 
neuper@37950:         a^2    a
neuper@37950:         --- = ---
neuper@37950:          a     1
neuper@37950: 
neuper@37950: Remark: the reverse ruleset does _NOT_ work properly with other input !.*)
neuper@37950: (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)
neuper@37950: 
neuper@37950: val {rules, rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "make_polynomial");
neuper@37950: (*WN060829 ... make_deriv does not terminate with 1st expl above,
neuper@37950:            see rational.sml --- investigate rulesets for cancel_p ---*)
neuper@37950: val {rules, rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "rev_rew_p");
neuper@37950: 
neuper@37950: (*.init_state = fn : term -> istate
neuper@37950: initialzies the state of the script interpreter. The state is:
neuper@37950: 
neuper@37950: type rrlsstate =      (*state for reverse rewriting*)
neuper@37950:      (term *          (*the current formula*)
neuper@37950:       term *          (*the final term*)
neuper@37950:       rule list       (*'reverse rule list' (#)*)
neuper@37950: 	    list *    (*may be serveral, eg. in norm_rational*)
neuper@37950:       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)
neuper@37950:        (term *        (*... rewrite with ...*)
neuper@37950: 	term list))   (*... assumptions*)
neuper@37950: 	  list);      (*derivation from given term to normalform
neuper@37950: 		       in reverse order with sym_thm;
neuper@37950:                        (#) could be extracted from here by (map #1)*).*)
neuper@37950: (* val {rules, rew_ord=(_,ro),...} =
neuper@37950:        rep_rls (assoc_rls "rev_rew_p")        (*USE ALWAYS, SEE val cancel_p*);
neuper@37972:    val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*);
neuper@37950:    val t = t;
neuper@37950:    *)
neuper@37950: fun init_state thy eval_rls ro t =
neuper@37950:     let val SOME (t',_) = factout_p_ thy t
neuper@37950:         val SOME (t'',asm) = cancel_p_ thy t
neuper@37950:         val der = reverse_deriv thy eval_rls rules ro NONE t'
neuper@37950:         val der = der @ [(Thm ("real_mult_div_cancel2",
neuper@37969: 			       num_str @{thm real_mult_div_cancel2}),
neuper@37950: 			  (t'',asm))]
neuper@37950:         val rs = (distinct_Thm o (map #1)) der
neuper@37950: 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
neuper@37950: 				      "sym_real_mult_0",
neuper@37950: 				      "sym_real_mult_1"
neuper@37950: 				      (*..insufficient,eg.make_Polynomial*)])rs
neuper@37950:     in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end;
neuper@37950: 
neuper@37950: (*.locate_rule = fn : rule list -> term -> rule
neuper@37950: 		      -> (rule * (term * term list) option) list.
neuper@37950:   checks a rule R for being a cancel-rule, and if it is,
neuper@37950:   then return the list of rules (+ the terms they are rewriting to)
neuper@37950:   which need to be applied before R should be applied.
neuper@37950:   precondition: the rule is applicable to the argument-term.
neuper@37950: arguments:
neuper@37950:   rule list: the reverse rule list
neuper@37950:   -> term  : ... to which the rule shall be applied
neuper@37950:   -> rule  : ... to be applied to term
neuper@37950: value:
neuper@37950:   -> (rule           : a rule rewriting to ...
neuper@37950:       * (term        : ... the resulting term ...
neuper@37950:          * term list): ... with the assumptions ( //#0).
neuper@37950:       ) list         : there may be several such rules;
neuper@37950: 		       the list is empty, if the rule has nothing to do
neuper@37950: 		       with cancelation.*)
neuper@37950: (* val () = ();
neuper@37950:    *)
neuper@37950: fun locate_rule thy eval_rls ro [rs] t r =
neuper@37950:     if (id_of_thm r) mem (map (id_of_thm)) rs
neuper@37950:     then let val ropt =
neuper@37950: 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
neuper@37950: 	 in case ropt of
neuper@37950: 		SOME ta => [(r, ta)]
neuper@37950: 	      | NONE => (writeln("### locate_rule:  rewrite "^
neuper@37950: 				 (id_of_thm r)^" "^(term2str t)^" = NONE");
neuper@37950: 			 []) end
neuper@37950:     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])
neuper@37950:   | locate_rule _ _ _ _ _ _ =
neuper@37950:     raise error ("locate_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.next_rule = fn : rule list -> term -> rule option
neuper@37950:   for a given term return the next rules to be done for cancelling.
neuper@37950: arguments:
neuper@37950:   rule list     : the reverse rule list
neuper@37950:   term          : the term for which ...
neuper@37950: value:
neuper@37950:   -> rule option: ... this rule is appropriate for cancellation;
neuper@37950: 		  there may be no such rule (if the term is canceled already.*)
neuper@37972: (* val thy = thy;
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel;
neuper@37950:    val ([rs],t) = (rss,f);
neuper@37950:    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
neuper@37950: 
neuper@37972:    val (thy, [rs]) = (thy, revsets);
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel;
neuper@37950:    nex [rs] t;
neuper@37950:    *)
neuper@37950: fun next_rule thy eval_rls ro [rs] t =
neuper@37950:     let val der = make_deriv thy eval_rls rs ro NONE t;
neuper@37950:     in case der of
neuper@37950: (* val (_,r,_)::_ = der;
neuper@37950:    *)
neuper@37950: 	   (_,r,_)::_ => SOME r
neuper@37950: 	 | _ => NONE
neuper@37950:     end
neuper@37950:   | next_rule _ _ _ _ _ =
neuper@37950:     raise error ("next_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.val attach_form = f : rule list -> term -> term
neuper@37950: 			 -> (rule * (term * term list)) list
neuper@37950:   checks an input term TI, if it may belong to a current cancellation, by
neuper@37950:   trying to derive it from the given term TG.
neuper@37950: arguments:
neuper@37950:   term   : TG, the last one in the cancellation agreed upon by user + math-eng
neuper@37950:   -> term: TI, the next one input by the user
neuper@37950: value:
neuper@37950:   -> (rule           : the rule to be applied in order to reach TI
neuper@37950:       * (term        : ... obtained by applying the rule ...
neuper@37950:          * term list): ... and the respective assumptions.
neuper@37950:       ) list         : there may be several such rules;
neuper@37950:                        the list is empty, if the users term does not belong
neuper@37950: 		       to a cancellation of the term last agreed upon.*)
neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
neuper@37950:     []:(rule * (term * term list)) list;
neuper@37950: 
neuper@37950: in
neuper@37950: 
neuper@37950: val cancel_p =
neuper@37950:     Rrls {id = "cancel_p", prepat=[],
neuper@37950: 	  rew_ord=("ord_make_polynomial",
neuper@37972: 		   ord_make_polynomial false thy),
neuper@37950: 	  erls = rational_erls,
neuper@37950: 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),
neuper@37950: 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),
neuper@37950: 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),
neuper@37950: 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],
neuper@37950: 	  (*asm_thm=[("real_mult_div_cancel2","")],*)
neuper@37950: 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,
neuper@37950: 		     normal_form = cancel_p_ thy,
neuper@37950: 		     locate_rule = locate_rule thy Atools_erls ro,
neuper@37950: 		     next_rule   = next_rule thy Atools_erls ro,
neuper@37950: 		     attach_form = attach_form}}
neuper@37950: end;(*local*)
neuper@37950: 
neuper@37978: val ------------------------------------------------------+++++ = "55555";
neuper@37950: 
neuper@37950: local(*.ad (1) 'cancel'
neuper@37950: ------------------------------
neuper@37950: cancels a single fraction consisting of two (uni- or multivariate)
neuper@37950: polynomials WN0609???SK[3] into another such a fraction; examples:
neuper@37950: 
neuper@37950: 	   a^2 - b^2           a + b
neuper@37950:         -------------------- = ---------
neuper@37950: 	a^2 - 2*a*b + b^2      a - *b
neuper@37950: 
neuper@37950: Remark: the reverse ruleset does _NOT_ work properly with other input !.*)
neuper@37950: (*WN 24.8.02: wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)
neuper@37950: 
neuper@37950: (*
neuper@37950: val SOME (Rls {rules=rules,rew_ord=(_,ro),...}) = 
neuper@37950:     assoc'(!ruleset',"expand_binoms");
neuper@37950: *)
neuper@37950: val {rules=rules,rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "expand_binoms");
neuper@37972: val thy = thy;
neuper@37950: 
neuper@37950: fun init_state thy eval_rls ro t =
neuper@37950:     let val SOME (t',_) = factout_ thy t;
neuper@37950:         val SOME (t'',asm) = cancel_ thy t;
neuper@37950:         val der = reverse_deriv thy eval_rls rules ro NONE t';
neuper@37950:         val der = der @ [(Thm ("real_mult_div_cancel2",
neuper@37969: 			       num_str @{thm real_mult_div_cancel2}),
neuper@37950: 			  (t'',asm))]
neuper@37950:         val rs = map #1 der;
neuper@37950:     in (t,t'',[rs],der) end;
neuper@37950: 
neuper@37950: fun locate_rule thy eval_rls ro [rs] t r =
neuper@37950:     if (id_of_thm r) mem (map (id_of_thm)) rs
neuper@37950:     then let val ropt = 
neuper@37950: 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
neuper@37950: 	 in case ropt of
neuper@37950: 		SOME ta => [(r, ta)]
neuper@37950: 	      | NONE => (writeln("### locate_rule:  rewrite "^
neuper@37950: 				 (id_of_thm r)^" "^(term2str t)^" = NONE");
neuper@37950: 			 []) end
neuper@37950:     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])
neuper@37950:   | locate_rule _ _ _ _ _ _ = 
neuper@37950:     raise error ("locate_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: fun next_rule thy eval_rls ro [rs] t =
neuper@37950:     let val der = make_deriv thy eval_rls rs ro NONE t;
neuper@37950:     in case der of 
neuper@37950: (* val (_,r,_)::_ = der;
neuper@37950:    *)
neuper@37950: 	   (_,r,_)::_ => SOME r
neuper@37950: 	 | _ => NONE
neuper@37950:     end
neuper@37950:   | next_rule _ _ _ _ _ = 
neuper@37950:     raise error ("next_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
neuper@37950:     []:(rule * (term * term list)) list;
neuper@37950: 
neuper@37950: val pat = (term_of o the o (parse thy)) "?r/?s";
neuper@37950: val pre1 = (term_of o the o (parse thy)) "?r is_expanded";
neuper@37950: val pre2 = (term_of o the o (parse thy)) "?s is_expanded";
neuper@37950: val prepat = [([pre1, pre2], pat)];
neuper@37950: 
neuper@37950: in
neuper@37950: 
neuper@37950: 
neuper@37950: val cancel = 
neuper@37950:     Rrls {id = "cancel", prepat=prepat,
neuper@37950: 	  rew_ord=("ord_make_polynomial",
neuper@37972: 		   ord_make_polynomial false thy),
neuper@37950: 	  erls = rational_erls, 
neuper@37950: 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),
neuper@37950: 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),
neuper@37950: 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),
neuper@37950: 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],
neuper@37950: 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,
neuper@37950: 		     normal_form = cancel_ thy, 
neuper@37950: 		     locate_rule = locate_rule thy Atools_erls ro,
neuper@37950: 		     next_rule   = next_rule thy Atools_erls ro,
neuper@37950: 		     attach_form = attach_form}}
neuper@37950: end;(*local*)
neuper@37950: 
neuper@37978: val ------------------------------------------------------+++++-+ = "66666";
neuper@37950: 
neuper@37950: 
neuper@37950: local(*.ad [2] 'common_nominator_p'
neuper@37950: ---------------------------------
neuper@37950: FIXME Beschreibung .*)
neuper@37950: 
neuper@37950: 
neuper@37950: val {rules=rules,rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "make_polynomial");
neuper@37950: (*WN060829 ... make_deriv does not terminate with 1st expl above,
neuper@37950:            see rational.sml --- investigate rulesets for cancel_p ---*)
neuper@37950: val {rules, rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "rev_rew_p");
neuper@37972: val thy = thy;
neuper@37950: 
neuper@37950: 
neuper@37950: (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option
neuper@37950:   as defined above*)
neuper@37950: 
neuper@37950: (*.init_state = fn : term -> istate
neuper@37950: initialzies the state of the interactive interpreter. The state is:
neuper@37950: 
neuper@37950: type rrlsstate =      (*state for reverse rewriting*)
neuper@37950:      (term *          (*the current formula*)
neuper@37950:       term *          (*the final term*)
neuper@37950:       rule list       (*'reverse rule list' (#)*)
neuper@37950: 	    list *    (*may be serveral, eg. in norm_rational*)
neuper@37950:       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)
neuper@37950:        (term *        (*... rewrite with ...*)
neuper@37950: 	term list))   (*... assumptions*)
neuper@37950: 	  list);      (*derivation from given term to normalform
neuper@37950: 		       in reverse order with sym_thm;
neuper@37950:                        (#) could be extracted from here by (map #1)*).*)
neuper@37950: fun init_state thy eval_rls ro t =
neuper@37950:     let val SOME (t',_) = common_nominator_p_ thy t;
neuper@37950:         val SOME (t'',asm) = add_fraction_p_ thy t;
neuper@37950:         val der = reverse_deriv thy eval_rls rules ro NONE t';
neuper@37950:         val der = der @ [(Thm ("real_mult_div_cancel2",
neuper@37969: 			       num_str @{thm real_mult_div_cancel2}),
neuper@37950: 			  (t'',asm))]
neuper@37950:         val rs = (distinct_Thm o (map #1)) der;
neuper@37950: 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
neuper@37950: 				      "sym_real_mult_0",
neuper@37950: 				      "sym_real_mult_1"]) rs;
neuper@37950:     in (t,t'',[rs(*here only _ONE_*)],der) end;
neuper@37950: 
neuper@37950: (* use"knowledge/Rational.ML";
neuper@37950:    *)
neuper@37950: 
neuper@37950: (*.locate_rule = fn : rule list -> term -> rule
neuper@37950: 		      -> (rule * (term * term list) option) list.
neuper@37950:   checks a rule R for being a cancel-rule, and if it is,
neuper@37950:   then return the list of rules (+ the terms they are rewriting to)
neuper@37950:   which need to be applied before R should be applied.
neuper@37950:   precondition: the rule is applicable to the argument-term.
neuper@37950: arguments:
neuper@37950:   rule list: the reverse rule list
neuper@37950:   -> term  : ... to which the rule shall be applied
neuper@37950:   -> rule  : ... to be applied to term
neuper@37950: value:
neuper@37950:   -> (rule           : a rule rewriting to ...
neuper@37950:       * (term        : ... the resulting term ...
neuper@37950:          * term list): ... with the assumptions ( //#0).
neuper@37950:       ) list         : there may be several such rules;
neuper@37950: 		       the list is empty, if the rule has nothing to do
neuper@37950: 		       with cancelation.*)
neuper@37950: (* val () = ();
neuper@37950:    *)
neuper@37950: fun locate_rule thy eval_rls ro [rs] t r =
neuper@37950:     if (id_of_thm r) mem (map (id_of_thm)) rs
neuper@37950:     then let val ropt =
neuper@37950: 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
neuper@37950: 	 in case ropt of
neuper@37950: 		SOME ta => [(r, ta)]
neuper@37950: 	      | NONE => (writeln("### locate_rule:  rewrite "^
neuper@37950: 				 (id_of_thm r)^" "^(term2str t)^" = NONE");
neuper@37950: 			 []) end
neuper@37950:     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])
neuper@37950:   | locate_rule _ _ _ _ _ _ =
neuper@37950:     raise error ("locate_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.next_rule = fn : rule list -> term -> rule option
neuper@37950:   for a given term return the next rules to be done for cancelling.
neuper@37950: arguments:
neuper@37950:   rule list     : the reverse rule list
neuper@37950:   term          : the term for which ...
neuper@37950: value:
neuper@37950:   -> rule option: ... this rule is appropriate for cancellation;
neuper@37950: 		  there may be no such rule (if the term is canceled already.*)
neuper@37972: (* val thy = thy;
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel;
neuper@37950:    val ([rs],t) = (rss,f);
neuper@37950:    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
neuper@37950: 
neuper@37972:    val (thy, [rs]) = (thy, revsets);
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel;
neuper@37950:    nex [rs] t;
neuper@37950:    *)
neuper@37950: fun next_rule thy eval_rls ro [rs] t =
neuper@37950:     let val der = make_deriv thy eval_rls rs ro NONE t;
neuper@37950:     in case der of
neuper@37950: (* val (_,r,_)::_ = der;
neuper@37950:    *)
neuper@37950: 	   (_,r,_)::_ => SOME r
neuper@37950: 	 | _ => NONE
neuper@37950:     end
neuper@37950:   | next_rule _ _ _ _ _ =
neuper@37950:     raise error ("next_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.val attach_form = f : rule list -> term -> term
neuper@37950: 			 -> (rule * (term * term list)) list
neuper@37950:   checks an input term TI, if it may belong to a current cancellation, by
neuper@37950:   trying to derive it from the given term TG.
neuper@37950: arguments:
neuper@37950:   term   : TG, the last one in the cancellation agreed upon by user + math-eng
neuper@37950:   -> term: TI, the next one input by the user
neuper@37950: value:
neuper@37950:   -> (rule           : the rule to be applied in order to reach TI
neuper@37950:       * (term        : ... obtained by applying the rule ...
neuper@37950:          * term list): ... and the respective assumptions.
neuper@37950:       ) list         : there may be several such rules;
neuper@37950:                        the list is empty, if the users term does not belong
neuper@37950: 		       to a cancellation of the term last agreed upon.*)
neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
neuper@37950:     []:(rule * (term * term list)) list;
neuper@37950: 
neuper@37950: val pat0 = (term_of o the o (parse thy)) "?r/?s+?u/?v";
neuper@37950: val pat1 = (term_of o the o (parse thy)) "?r/?s+?u   ";
neuper@37950: val pat2 = (term_of o the o (parse thy)) "?r   +?u/?v";
neuper@37950: val prepat = [([HOLogic.true_const], pat0),
neuper@37950: 	      ([HOLogic.true_const], pat1),
neuper@37950: 	      ([HOLogic.true_const], pat2)];
neuper@37950: 
neuper@37950: in
neuper@37950: 
neuper@37950: (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;
neuper@37950:   besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht
neuper@37950:   dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)
neuper@37950: val common_nominator_p =
neuper@37950:     Rrls {id = "common_nominator_p", prepat=prepat,
neuper@37950: 	  rew_ord=("ord_make_polynomial",
neuper@37972: 		   ord_make_polynomial false thy),
neuper@37950: 	  erls = rational_erls,
neuper@37950: 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),
neuper@37950: 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),
neuper@37950: 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),
neuper@37950: 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],
neuper@37950: 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,
neuper@37950: 		     normal_form = add_fraction_p_ thy,(*FIXME.WN0211*)
neuper@37950: 		     locate_rule = locate_rule thy Atools_erls ro,
neuper@37950: 		     next_rule   = next_rule thy Atools_erls ro,
neuper@37950: 		     attach_form = attach_form}}
neuper@37950: end;(*local*)
neuper@37950: 
neuper@37978: val ------------------------------------------------------+++++-++ = "77777";
neuper@37950: 
neuper@37950: local(*.ad [2] 'common_nominator'
neuper@37950: ---------------------------------
neuper@37950: FIXME Beschreibung .*)
neuper@37950: 
neuper@37950: 
neuper@37950: val {rules=rules,rew_ord=(_,ro),...} =
neuper@37950:     rep_rls (assoc_rls "make_polynomial");
neuper@37972: val thy = thy;
neuper@37950: 
neuper@37950: 
neuper@37950: (*.common_nominator_ = fn : theory -> term -> (term * term list) option
neuper@37950:   as defined above*)
neuper@37950: 
neuper@37950: (*.init_state = fn : term -> istate
neuper@37950: initialzies the state of the interactive interpreter. The state is:
neuper@37950: 
neuper@37978: val ------------------------------------------------------+++++-++++ = "99999";
neuper@37978: val ------------------------------------------------------+++++-+++++ = "00000";
neuper@37978: 
neuper@37950: type rrlsstate =      (*state for reverse rewriting*)
neuper@37950:      (term *          (*the current formula*)
neuper@37950:       term *          (*the final term*)
neuper@37950:       rule list       (*'reverse rule list' (#)*)
neuper@37950: 	    list *    (*may be serveral, eg. in norm_rational*)
neuper@37950:       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)
neuper@37950:        (term *        (*... rewrite with ...*)
neuper@37950: 	term list))   (*... assumptions*)
neuper@37950: 	  list);      (*derivation from given term to normalform
neuper@37950: 		       in reverse order with sym_thm;
neuper@37950:                        (#) could be extracted from here by (map #1)*).*)
neuper@37950: fun init_state thy eval_rls ro t =
neuper@37950:     let val SOME (t',_) = common_nominator_ thy t;
neuper@37950:         val SOME (t'',asm) = add_fraction_ thy t;
neuper@37950:         val der = reverse_deriv thy eval_rls rules ro NONE t';
neuper@37950:         val der = der @ [(Thm ("real_mult_div_cancel2",
neuper@37969: 			       num_str @{thm real_mult_div_cancel2}),
neuper@37950: 			  (t'',asm))]
neuper@37950:         val rs = (distinct_Thm o (map #1)) der;
neuper@37950: 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",
neuper@37950: 				      "sym_real_mult_0",
neuper@37950: 				      "sym_real_mult_1"]) rs;
neuper@37950:     in (t,t'',[rs(*here only _ONE_*)],der) end;
neuper@37950: 
neuper@37950: (*.locate_rule = fn : rule list -> term -> rule
neuper@37950: 		      -> (rule * (term * term list) option) list.
neuper@37950:   checks a rule R for being a cancel-rule, and if it is,
neuper@37950:   then return the list of rules (+ the terms they are rewriting to)
neuper@37950:   which need to be applied before R should be applied.
neuper@37950:   precondition: the rule is applicable to the argument-term.
neuper@37950: arguments:
neuper@37950:   rule list: the reverse rule list
neuper@37950:   -> term  : ... to which the rule shall be applied
neuper@37950:   -> rule  : ... to be applied to term
neuper@37950: value:
neuper@37950:   -> (rule           : a rule rewriting to ...
neuper@37950:       * (term        : ... the resulting term ...
neuper@37950:          * term list): ... with the assumptions ( //#0).
neuper@37950:       ) list         : there may be several such rules;
neuper@37950: 		       the list is empty, if the rule has nothing to do
neuper@37950: 		       with cancelation.*)
neuper@37950: (* val () = ();
neuper@37950:    *)
neuper@37950: fun locate_rule thy eval_rls ro [rs] t r =
neuper@37950:     if (id_of_thm r) mem (map (id_of_thm)) rs
neuper@37950:     then let val ropt =
neuper@37950: 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;
neuper@37950: 	 in case ropt of
neuper@37950: 		SOME ta => [(r, ta)]
neuper@37950: 	      | NONE => (writeln("### locate_rule:  rewrite "^
neuper@37950: 				 (id_of_thm r)^" "^(term2str t)^" = NONE");
neuper@37950: 			 []) end
neuper@37950:     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])
neuper@37950:   | locate_rule _ _ _ _ _ _ =
neuper@37950:     raise error ("locate_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.next_rule = fn : rule list -> term -> rule option
neuper@37950:   for a given term return the next rules to be done for cancelling.
neuper@37950: arguments:
neuper@37950:   rule list     : the reverse rule list
neuper@37950:   term          : the term for which ...
neuper@37950: value:
neuper@37950:   -> rule option: ... this rule is appropriate for cancellation;
neuper@37950: 		  there may be no such rule (if the term is canceled already.*)
neuper@37972: (* val thy = thy;
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel;
neuper@37950:    val ([rs],t) = (rss,f);
neuper@37950:    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)
neuper@37950: 
neuper@37972:    val (thy, [rs]) = (thy, revsets);
neuper@37950:    val Rrls {rew_ord=(_,ro),...} = cancel_p;
neuper@37950:    nex [rs] t;
neuper@37950:    *)
neuper@37950: fun next_rule thy eval_rls ro [rs] t =
neuper@37950:     let val der = make_deriv thy eval_rls rs ro NONE t;
neuper@37950:     in case der of
neuper@37950: (* val (_,r,_)::_ = der;
neuper@37950:    *)
neuper@37950: 	   (_,r,_)::_ => SOME r
neuper@37950: 	 | _ => NONE
neuper@37950:     end
neuper@37950:   | next_rule _ _ _ _ _ =
neuper@37950:     raise error ("next_rule: doesnt match rev-sets in istate");
neuper@37950: 
neuper@37950: (*.val attach_form = f : rule list -> term -> term
neuper@37950: 			 -> (rule * (term * term list)) list
neuper@37950:   checks an input term TI, if it may belong to a current cancellation, by
neuper@37950:   trying to derive it from the given term TG.
neuper@37950: arguments:
neuper@37950:   term   : TG, the last one in the cancellation agreed upon by user + math-eng
neuper@37950:   -> term: TI, the next one input by the user
neuper@37950: value:
neuper@37950:   -> (rule           : the rule to be applied in order to reach TI
neuper@37950:       * (term        : ... obtained by applying the rule ...
neuper@37950:          * term list): ... and the respective assumptions.
neuper@37950:       ) list         : there may be several such rules;
neuper@37950:                        the list is empty, if the users term does not belong
neuper@37950: 		       to a cancellation of the term last agreed upon.*)
neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)
neuper@37950:     []:(rule * (term * term list)) list;
neuper@37950: 
neuper@37950: val pat0 =  (term_of o the o (parse thy)) "?r/?s+?u/?v";
neuper@37950: val pat01 = (term_of o the o (parse thy)) "?r/?s-?u/?v";
neuper@37950: val pat1 =  (term_of o the o (parse thy)) "?r/?s+?u   ";
neuper@37950: val pat11 = (term_of o the o (parse thy)) "?r/?s-?u   ";
neuper@37950: val pat2 =  (term_of o the o (parse thy)) "?r   +?u/?v";
neuper@37950: val pat21 = (term_of o the o (parse thy)) "?r   -?u/?v";
neuper@37950: val prepat = [([HOLogic.true_const], pat0),
neuper@37950: 	      ([HOLogic.true_const], pat01),
neuper@37950: 	      ([HOLogic.true_const], pat1),
neuper@37950: 	      ([HOLogic.true_const], pat11),
neuper@37950: 	      ([HOLogic.true_const], pat2),
neuper@37950: 	      ([HOLogic.true_const], pat21)];
neuper@37950: 
neuper@37950: 
neuper@37950: in
neuper@37950: 
neuper@37950: val common_nominator =
neuper@37950:     Rrls {id = "common_nominator", prepat=prepat,
neuper@37950: 	  rew_ord=("ord_make_polynomial",
neuper@37972: 		   ord_make_polynomial false thy),
neuper@37950: 	  erls = rational_erls,
neuper@37950: 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),
neuper@37950: 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),
neuper@37950: 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),
neuper@37950: 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],
neuper@37950: 	  (*asm_thm=[("real_mult_div_cancel2","")],*)
neuper@37950: 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,
neuper@37950: 		     normal_form = add_fraction_ (*NOT common_nominator_*) thy,
neuper@37950: 		     locate_rule = locate_rule thy Atools_erls ro,
neuper@37950: 		     next_rule   = next_rule thy Atools_erls ro,
neuper@37950: 		     attach_form = attach_form}}
neuper@37950: 
neuper@37950: end;(*local*)
neuper@37950: end;(*struct*)
neuper@37950: 
neuper@37978: val ------------------------------------------------------+++++-+++ = "88888";
neuper@37978: 
neuper@37950: open RationalI;
neuper@37950: (*##*)
neuper@37950: 
neuper@37950: (*.the expression contains + - * ^ / only ?.*)
neuper@37950: fun is_ratpolyexp (Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("op +",_) $ Free _ $ Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("op -",_) $ Free _ $ Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("op *",_) $ Free _ $ Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("HOL.divide",_) $ Free _ $ Free _) = true
neuper@37950:   | is_ratpolyexp (Const ("op +",_) $ t1 $ t2) = 
neuper@37950:                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
neuper@37950:   | is_ratpolyexp (Const ("op -",_) $ t1 $ t2) = 
neuper@37950:                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
neuper@37950:   | is_ratpolyexp (Const ("op *",_) $ t1 $ t2) = 
neuper@37950:                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
neuper@37950:   | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) = 
neuper@37950:                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
neuper@37950:   | is_ratpolyexp (Const ("HOL.divide",_) $ t1 $ t2) = 
neuper@37950:                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))
neuper@37950:   | is_ratpolyexp _ = false;
neuper@37950: 
neuper@37950: (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)
neuper@37950: fun eval_is_ratpolyexp (thmid:string) _ 
neuper@37950: 		       (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =
neuper@37950:     if is_ratpolyexp arg
neuper@37950:     then SOME (mk_thmid thmid "" 
neuper@37950: 			((Syntax.string_of_term (thy2ctxt thy)) arg) "", 
neuper@37950: 	       Trueprop $ (mk_equality (t, HOLogic.true_const)))
neuper@37950:     else SOME (mk_thmid thmid "" 
neuper@37950: 			((Syntax.string_of_term (thy2ctxt thy)) arg) "", 
neuper@37950: 	       Trueprop $ (mk_equality (t, HOLogic.false_const)))
neuper@37950:   | eval_is_ratpolyexp _ _ _ _ = NONE; 
neuper@37950: 
neuper@37950: 
neuper@37950: 
neuper@37950: (*-------------------18.3.03 --> struct <-----------vvv--*)
neuper@37950: val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)
neuper@37950: 
neuper@37950: (*.discard binary minus, shift unary minus into -1*; 
neuper@37950:    unary minus before numerals are put into the numeral by parsing;
neuper@37950:    contains absolute minimum of thms for context in norm_Rational .*)
neuper@37950: val discard_minus = prep_rls(
neuper@37950:   Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37969:       rules = [Thm ("real_diff_minus", num_str @{thm real_diff_minus}),
neuper@37950: 	       (*"a - b = a + -1 * b"*)
neuper@37969: 	       Thm ("sym_real_mult_minus1",
neuper@37969:                      num_str (@{thm real_mult_minus1} RS @{thm sym}))
neuper@37950: 	       (*- ?z = "-1 * ?z"*)
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) 
neuper@37950:       "empty_script")
neuper@37950:       }):rls;
neuper@37950: (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)
neuper@37950: val powers_erls = prep_rls(
neuper@37950:   Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37950:       rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),
neuper@37950: 	       Calc ("Atools.is'_even",eval_is_even "#is_even_"),
neuper@37950: 	       Calc ("op <",eval_equ "#less_"),
neuper@37950: 	       Thm ("not_false", not_false),
neuper@37950: 	       Thm ("not_true", not_true),
neuper@37950: 	       Calc ("op +",eval_binop "#add_")
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) 
neuper@37950:       "empty_script")
neuper@37950:       }:rls);
neuper@37950: (*.all powers over + distributed; atoms over * collected, other distributed
neuper@37950:    contains absolute minimum of thms for context in norm_Rational .*)
neuper@37950: val powers = prep_rls(
neuper@37950:   Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = powers_erls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37969:       rules = [Thm ("realpow_multI", num_str @{thm realpow_multI}),
neuper@37950: 	       (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)
neuper@37969: 	       Thm ("realpow_pow",num_str @{thm realpow_pow}),
neuper@37950: 	       (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)
neuper@37969: 	       Thm ("realpow_oneI",num_str @{thm realpow_oneI}),
neuper@37950: 	       (*"r ^^^ 1 = r"*)
neuper@37969: 	       Thm ("realpow_minus_even",num_str @{thm realpow_minus_even}),
neuper@37950: 	       (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)
neuper@37969: 	       Thm ("realpow_minus_odd",num_str @{thm realpow_minus_odd}),
neuper@37950: 	       (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)
neuper@37950: 	       
neuper@37950: 	       (*----- collect atoms over * -----*)
neuper@37969: 	       Thm ("realpow_two_atom",num_str @{thm realpow_two_atom}),	
neuper@37950: 	       (*"r is_atom ==> r * r = r ^^^ 2"*)
neuper@37969: 	       Thm ("realpow_plus_1",num_str @{thm realpow_plus_1}),		
neuper@37950: 	       (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)
neuper@37969: 	       Thm ("realpow_addI_atom",num_str @{thm realpow_addI_atom}),
neuper@37950: 	       (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)
neuper@37950: 
neuper@37950: 	       (*----- distribute none-atoms -----*)
neuper@37969: 	       Thm ("realpow_def_atom",num_str @{thm realpow_def_atom}),
neuper@37950: 	       (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)
neuper@37969: 	       Thm ("realpow_eq_oneI",num_str @{thm realpow_eq_oneI}),
neuper@37950: 	       (*"1 ^^^ n = 1"*)
neuper@37950: 	       Calc ("op +",eval_binop "#add_")
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) 
neuper@37950:       "empty_script")
neuper@37950:       }:rls);
neuper@37950: (*.contains absolute minimum of thms for context in norm_Rational.*)
neuper@37950: val rat_mult_divide = prep_rls(
neuper@37950:   Rls {id = "rat_mult_divide", preconds = [], 
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37969:       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),
neuper@37950: 	       (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
neuper@37965: 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),
neuper@37950: 	       (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],
neuper@37950: 	       otherwise inv.to a / b / c = ...*)
neuper@37965: 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),
neuper@37950: 	       (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much
neuper@37950: 		     and does not commute a / b * c ^^^ 2 !*)
neuper@37950: 	       
neuper@37965: 	       Thm ("divide_divide_eq_right", real_divide_divide1_eq),
neuper@37950: 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
neuper@37965: 	       Thm ("divide_divide_eq_left", real_divide_divide2_eq),
neuper@37950: 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
neuper@37950: 	       Calc ("HOL.divide"  ,eval_cancel "#divide_")
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37950:       }:rls);
neuper@37950: (*.contains absolute minimum of thms for context in norm_Rational.*)
neuper@37950: val reduce_0_1_2 = prep_rls(
neuper@37950:   Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),
neuper@37950:       erls = e_rls,srls = Erls,calc = [],(*asm_thm = [],*)
neuper@37965:       rules = [(*Thm ("divide_1",num_str @{thm divide_1}),
neuper@37950: 		 "?x / 1 = ?x" unnecess.for normalform*)
neuper@37965: 	       Thm ("mult_1_left",num_str @{thm mult_1_left}),                 
neuper@37950: 	       (*"1 * z = z"*)
neuper@37969: 	       (*Thm ("real_mult_minus1",num_str @{thm real_mult_minus1}),
neuper@37950: 	       "-1 * z = - z"*)
neuper@37969: 	       (*Thm ("real_minus_mult_cancel",num_str @{thm real_minus_mult_cancel}),
neuper@37950: 	       "- ?x * - ?y = ?x * ?y"*)
neuper@37950: 
neuper@37965: 	       Thm ("mult_zero_left",num_str @{thm mult_zero_left}),        
neuper@37950: 	       (*"0 * z = 0"*)
neuper@37965: 	       Thm ("add_0_left",num_str @{thm add_0_left}),
neuper@37950: 	       (*"0 + z = z"*)
neuper@37965: 	       (*Thm ("right_minus",num_str @{thm right_minus}),
neuper@37950: 	       "?z + - ?z = 0"*)
neuper@37950: 
neuper@37969: 	       Thm ("sym_real_mult_2",
neuper@37969:                      num_str (@{thm real_mult_2} RS @{thm sym})),	
neuper@37950: 	       (*"z1 + z1 = 2 * z1"*)
neuper@37969: 	       Thm ("real_mult_2_assoc",num_str @{thm real_mult_2_assoc}),
neuper@37950: 	       (*"z1 + (z1 + k) = 2 * z1 + k"*)
neuper@37950: 
neuper@37965: 	       Thm ("divide_zero_left",num_str @{thm divide_zero_left})
neuper@37950: 	       (*"0 / ?x = 0"*)
neuper@37950: 	       ], scr = EmptyScr}:rls);
neuper@37950: 
neuper@37950: (*erls for calculate_Rational; 
neuper@37950:   make local with FIXX@ME result:term *term list WN0609???SKMG*)
neuper@37950: val norm_rat_erls = prep_rls(
neuper@37950:   Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37950:       rules = [Calc ("Atools.is'_const",eval_const "#is_const_")
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) 
neuper@37950:       "empty_script")
neuper@37950:       }:rls);
neuper@37950: (*.consists of rls containing the absolute minimum of thms.*)
neuper@37950: (*040209: this version has been used by RL for his equations,
neuper@37950: which is now replaced by MGs version below
neuper@37950: vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)
neuper@37950: val norm_Rational = prep_rls(
neuper@37950:   Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:       erls = norm_rat_erls, srls = Erls, calc = [], (*asm_thm = [],*)
neuper@37950:       rules = [(*sequence given by operator precedence*)
neuper@37950: 	       Rls_ discard_minus,
neuper@37950: 	       Rls_ powers,
neuper@37950: 	       Rls_ rat_mult_divide,
neuper@37950: 	       Rls_ expand,
neuper@37950: 	       Rls_ reduce_0_1_2,
neuper@37950: 	       (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)
neuper@37950: 	       Rls_ order_add_mult,
neuper@37950: 	       Rls_ collect_numerals,
neuper@37950: 	       Rls_ add_fractions_p,
neuper@37950: 	       Rls_ cancel_p
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) 
neuper@37950:       "empty_script")
neuper@37950:       }:rls);
neuper@37950: val norm_Rational_parenthesized = prep_rls(
neuper@37950:   Seq {id = "norm_Rational_parenthesized", preconds = []:term list, 
neuper@37950:        rew_ord = ("dummy_ord", dummy_ord),
neuper@37950:       erls = Atools_erls, srls = Erls,
neuper@37950:       calc = [], (*asm_thm = [],*)
neuper@37950:       rules = [Rls_  norm_Rational, (*from RL -- not the latest one*)
neuper@37950: 	       Rls_ discard_parentheses
neuper@37950: 	       ],
neuper@37950:       scr = EmptyScr
neuper@37950:       }:rls);      
neuper@37950: 
neuper@37950: 
neuper@37950: (*-------------------18.3.03 --> struct <-----------^^^--*)
neuper@37950: 
neuper@37950: 
neuper@37950: (*WN030318???SK: simplifies all but cancel and common_nominator*)
neuper@37950: val simplify_rational = 
neuper@37950:     merge_rls "simplify_rational" expand_binoms
neuper@37950:     (append_rls "divide" calculate_Rational
neuper@37965: 		[Thm ("divide_1",num_str @{thm divide_1}),
neuper@37950: 		 (*"?x / 1 = ?x"*)
neuper@37978: 		 Thm ("rat_mult",num_str @{thm rat_mult}),
neuper@37950: 		 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
neuper@37965: 		 Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),
neuper@37950: 		 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],
neuper@37950: 		 otherwise inv.to a / b / c = ...*)
neuper@37965: 		 Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),
neuper@37950: 		 (*"?a / ?b * ?c = ?a * ?c / ?b"*)
neuper@37969: 		 Thm ("add_minus",num_str @{thm add_minus}),
neuper@37950: 		 (*"?a + ?b - ?b = ?a"*)
neuper@37969: 		 Thm ("add_minus1",num_str @{thm add_minus1}),
neuper@37950: 		 (*"?a - ?b + ?b = ?a"*)
neuper@37978: 		 Thm ("divide_minus1",num_str @{thm divide_minus1})
neuper@37950: 		 (*"?x / -1 = - ?x"*)
neuper@37950: (*
neuper@37950: ,
neuper@37969: 		 Thm ("",num_str @{thm })
neuper@37950: *)
neuper@37950: 		 ]);
neuper@37950: 
neuper@37950: (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)
neuper@37950: 
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: (*                  Simplifier für beliebige Buchterme                *) 
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: (*----------------------- norm_Rational_mg ---------------------------*)
neuper@37950: (*. description of the simplifier see MG-DA.p.56ff .*)
neuper@37950: (* ------------------------------------------------------------------- *)
neuper@37950: val common_nominator_p_rls = prep_rls(
neuper@37950:   Rls {id = "common_nominator_p_rls", preconds = [],
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950: 	 erls = e_rls, srls = Erls, calc = [],
neuper@37950: 	 rules = 
neuper@37950: 	 [Rls_ common_nominator_p
neuper@37950: 	  (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?
neuper@37950: 	    FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)
neuper@37950: 	  ], 
neuper@37950: 	 scr = EmptyScr});
neuper@37950: (* ------------------------------------------------------------------- *)
neuper@37950: val cancel_p_rls = prep_rls(
neuper@37950:   Rls {id = "cancel_p_rls", preconds = [],
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950: 	 erls = e_rls, srls = Erls, calc = [],
neuper@37950: 	 rules = 
neuper@37950: 	 [Rls_ cancel_p
neuper@37950: 	  (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)
neuper@37950: 	  ], 
neuper@37950: 	 scr = EmptyScr});
neuper@37950: (* -------------------------------------------------------------------- *)
neuper@37950: (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
neuper@37950:     used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)
neuper@37950: val rat_mult_poly = prep_rls(
neuper@37950:   Rls {id = "rat_mult_poly", preconds = [],
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950: 	 erls =  append_rls "e_rls-is_polyexp" e_rls
neuper@37950: 	         [Calc ("Poly.is'_polyexp", eval_is_polyexp "")], 
neuper@37950: 	 srls = Erls, calc = [],
neuper@37950: 	 rules = 
neuper@37969: 	 [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),
neuper@37950: 	  (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
neuper@37969: 	  Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})
neuper@37950: 	  (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
neuper@37950: 	  ], 
neuper@37950: 	 scr = EmptyScr});
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;
neuper@37950:     used in looping part norm_Rational_rls, see example DA-M02-main.p.60 
neuper@37950:     .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls, 
neuper@37950:     I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028 
neuper@37950:     ... WN0609???MG.*)
neuper@37950: val rat_mult_div_pow = prep_rls(
neuper@37950:   Rls {id = "rat_mult_div_pow", preconds = [], 
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:        erls = e_rls,
neuper@37950:        (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls
neuper@37950: 			[Calc ("Poly.is'_polyexp", eval_is_polyexp "")],
neuper@37950:          with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get 
neuper@37950: 	 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.
neuper@37950:          thus we decided to go on with this flaw*)
neuper@37950: 		 srls = Erls, calc = [],
neuper@37969:       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),
neuper@37950: 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
neuper@37969: 	       Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),
neuper@37950: 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
neuper@37969: 	       Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}),
neuper@37950: 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
neuper@37950: 
neuper@37950: 	       Thm ("real_divide_divide1_mg", real_divide_divide1_mg),
neuper@37950: 	       (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)
neuper@37965: 	       Thm ("divide_divide_eq_right", real_divide_divide1_eq),
neuper@37950: 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
neuper@37965: 	       Thm ("divide_divide_eq_left", real_divide_divide2_eq),
neuper@37950: 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
neuper@37950: 	       Calc ("HOL.divide"  ,eval_cancel "#divide_"),
neuper@37950: 	      
neuper@37969: 	       Thm ("rat_power", num_str @{thm rat_power})
neuper@37950: 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
neuper@37950: 	       ],
neuper@37950:       scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37950:       }:rls);
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: val rat_reduce_1 = prep_rls(
neuper@37950:   Rls {id = "rat_reduce_1", preconds = [], 
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:        erls = e_rls, srls = Erls, calc = [], 
neuper@37965:        rules = [Thm ("divide_1",num_str @{thm divide_1}),
neuper@37950: 		(*"?x / 1 = ?x"*)
neuper@37965: 		Thm ("mult_1_left",num_str @{thm mult_1_left})           
neuper@37950: 		(*"1 * z = z"*)
neuper@37950: 		],
neuper@37950:        scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37950:        }:rls);
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: (*. looping part of norm_Rational(*_mg*) .*)
neuper@37950: val norm_Rational_rls = prep_rls(
neuper@37950:    Rls {id = "norm_Rational_rls", preconds = [], 
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:        erls = norm_rat_erls, srls = Erls, calc = [],
neuper@37950:        rules = [Rls_ common_nominator_p_rls,
neuper@37950: 		Rls_ rat_mult_div_pow,
neuper@37950: 		Rls_ make_rat_poly_with_parentheses,
neuper@37950: 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)
neuper@37950: 		Rls_ rat_reduce_1
neuper@37950: 		],
neuper@37950:        scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37950:        }:rls);
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)
neuper@37950:  just be renaming:*)
neuper@37950: val norm_Rational(*_mg*) = prep_rls(
neuper@37950:    Seq {id = "norm_Rational"(*_mg*), preconds = [], 
neuper@37950:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37950:        erls = norm_rat_erls, srls = Erls, calc = [],
neuper@37950:        rules = [Rls_ discard_minus_,
neuper@37950: 		Rls_ rat_mult_poly,(* removes double fractions like a/b/c    *)
neuper@37950: 		Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)
neuper@37950: 		Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)
neuper@37950: 		Rls_ norm_Rational_rls,   (* the main rls, looping (#)       *)
neuper@37950: 		Rls_ discard_parentheses_ (* mult only                       *)
neuper@37950: 		],
neuper@37950:        scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37950:        }:rls);
neuper@37950: (* ------------------------------------------------------------------ *)
neuper@37950: 
neuper@37950: 
neuper@37967: ruleset' := overwritelthy @{theory} (!ruleset',
neuper@37950:   [("calculate_Rational", calculate_Rational),
neuper@37950:    ("calc_rat_erls",calc_rat_erls),
neuper@37950:    ("rational_erls", rational_erls),
neuper@37950:    ("cancel_p", cancel_p),
neuper@37950:    ("cancel", cancel),
neuper@37950:    ("common_nominator_p", common_nominator_p),
neuper@37950:    ("common_nominator_p_rls", common_nominator_p_rls),
neuper@37950:    ("common_nominator"  , common_nominator),
neuper@37950:    ("discard_minus", discard_minus),
neuper@37950:    ("powers_erls", powers_erls),
neuper@37950:    ("powers", powers),
neuper@37950:    ("rat_mult_divide", rat_mult_divide),
neuper@37950:    ("reduce_0_1_2", reduce_0_1_2),
neuper@37950:    ("rat_reduce_1", rat_reduce_1),
neuper@37950:    ("norm_rat_erls", norm_rat_erls),
neuper@37950:    ("norm_Rational", norm_Rational),
neuper@37950:    ("norm_Rational_rls", norm_Rational_rls),
neuper@37950:    ("norm_Rational_parenthesized", norm_Rational_parenthesized),
neuper@37950:    ("rat_mult_poly", rat_mult_poly),
neuper@37950:    ("rat_mult_div_pow", rat_mult_div_pow),
neuper@37950:    ("cancel_p_rls", cancel_p_rls)
neuper@37950:    ]);
neuper@37950: 
neuper@37950: calclist':= overwritel (!calclist', 
neuper@37950:    [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))
neuper@37950:     ]);
neuper@37950: 
neuper@37950: (** problems **)
neuper@37950: 
neuper@37950: store_pbt
neuper@37972:  (prep_pbt thy "pbl_simp_rat" [] e_pblID
neuper@37950:  (["rational","simplification"],
neuper@37967:   [("#Given" ,["TERM t_"]),
neuper@37950:    ("#Where" ,["t_ is_ratpolyexp"]),
neuper@37950:    ("#Find"  ,["normalform n_"])
neuper@37950:   ],
neuper@37950:   append_rls "e_rls" e_rls [(*for preds in where_*)], 
neuper@37950:   SOME "Simplify t_", 
neuper@37950:   [["simplification","of_rationals"]]));
neuper@37950: 
neuper@37950: (** methods **)
neuper@37950: 
neuper@37950: (*WN061025 this methods script is copied from (auto-generated) script
neuper@37950:   of norm_Rational in order to ease repair on inform*)
neuper@37950: store_met
neuper@37972:     (prep_met thy "met_simp_rat" [] e_metID
neuper@37950: 	      (["simplification","of_rationals"],
neuper@37967: 	       [("#Given" ,["TERM t_"]),
neuper@37950: 		("#Where" ,["t_ is_ratpolyexp"]),
neuper@37950: 		("#Find"  ,["normalform n_"])
neuper@37950: 		],
neuper@37950: 	       {rew_ord'="tless_true",
neuper@37950: 		rls' = e_rls,
neuper@37950: 		calc = [], srls = e_rls, 
neuper@37950: 		prls = append_rls "simplification_of_rationals_prls" e_rls 
neuper@37950: 				[(*for preds in where_*)
neuper@37950: 				 Calc ("Rational.is'_ratpolyexp", 
neuper@37950: 				       eval_is_ratpolyexp "")],
neuper@37950: 		crls = e_rls, nrls = norm_Rational_rls},
neuper@37950: "Script SimplifyScript (t_::real) =                              " ^
neuper@37950: "  ((Try (Rewrite_Set discard_minus_ False) @@                   " ^
neuper@37950: "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^
neuper@37950: "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^
neuper@37950: "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^
neuper@37950: "    (Repeat                                                     " ^
neuper@37950: "     ((Try (Rewrite_Set common_nominator_p_rls False) @@        " ^
neuper@37950: "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^
neuper@37950: "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^
neuper@37950: "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^
neuper@37950: "       Try (Rewrite_Set rat_reduce_1 False)))) @@               " ^
neuper@37950: "    Try (Rewrite_Set discard_parentheses_ False))               " ^
neuper@37950: "    t_)"
neuper@37950: 	       ));
neuper@37950: *}
neuper@37950: 
neuper@37950: end




