(* rationals, fractions of multivariate polynomials over the real field

    author: isac team

    Copyright (c) isac team 2002, 2013

    Use is subject to license terms.

    depends on Poly (and not on Atools), because 

    fractions with _normalized_ polynomials are canceled, added, etc.

 *)

 theory Rational 

 imports Poly "~~/src/Tools/isac/Knowledge/GCD_Poly_ML"

 begin

 section {* Constants for evaluation by "Calc" *}

 consts

   is'_expanded    :: "real => bool" ("_ is'_expanded")     (*RL->Poly.thy*)

   is'_ratpolyexp  :: "real => bool" ("_ is'_ratpolyexp") 

   get_denominator :: "real => real"

   get_numerator   :: "real => real"

 ML {*

 (*.the expression contains + - * ^ / only ?.*)

 fun is_ratpolyexp (Free _) = true

   | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ Free _ $ Free _) = true

   | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ Free _ $ Free _) = true

   | is_ratpolyexp (Const ("Groups.times_class.times",_) $ Free _ $ Free _) = true

   | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true

   | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ Free _ $ Free _) = true

   | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ t1 $ t2) = 

                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

   | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) = 

                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

   | is_ratpolyexp (Const ("Groups.times_class.times",_) $ t1 $ t2) = 

                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

   | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) = 

                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

   | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ t1 $ t2) = 

                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

   | is_ratpolyexp _ = false;

 (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)

 fun eval_is_ratpolyexp (thmid:string) _ 

 		       (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =

     if is_ratpolyexp arg

     then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "", 

 	         Trueprop $ (mk_equality (t, @{term True})))

     else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "", 

 	         Trueprop $ (mk_equality (t, @{term False})))

   | eval_is_ratpolyexp _ _ _ _ = NONE; 

 (*("get_denominator", ("Rational.get_denominator", eval_get_denominator ""))*)

 fun eval_get_denominator (thmid:string) _ 

 		      (t as Const ("Rational.get_denominator", _) $

               (Const ("Fields.inverse_class.divide", _) $ num $

                 denom)) thy = 

       SOME (mk_thmid thmid "" (term_to_string''' thy denom) "", 

 	            Trueprop $ (mk_equality (t, denom)))

   | eval_get_denominator _ _ _ _ = NONE; 

 (*("get_numerator", ("Rational.get_numerator", eval_get_numerator ""))*)

 fun eval_get_numerator (thmid:string) _ 

       (t as Const ("Rational.get_numerator", _) $

           (Const ("Fields.inverse_class.divide", _) $num

             $denom )) thy = 

     SOME (mk_thmid thmid "" (term_to_string''' thy num) "", 

 	    Trueprop $ (mk_equality (t, num)))

   | eval_get_numerator _ _ _ _ = NONE; 

 *}

 section {* Theorems for rewriting *}

 axiomatization (* naming due to Isabelle2002, but not contained in Isabelle2002; 

                   many thms are due to RL and can be removed with updating the equation solver;

                   TODO: replace by equivalent thms in recent Isabelle201x *) 

 where

   mult_cross:      "[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)" and

   mult_cross1:     "   b ~= 0            ==> (a / b = c    ) = (a     = b * c)" and

   mult_cross2:     "           d ~= 0    ==> (a     = c / d) = (a * d =     c)" and

   add_minus:       "a + b - b = a"(*RL->Poly.thy*) and

   add_minus1:      "a - b + b = a"(*RL->Poly.thy*) and

   rat_mult:        "a / b * (c / d) = a * c / (b * d)"(*?Isa02*)  and

   rat_mult2:       "a / b *  c      = a * c /  b     "(*?Isa02*) and

   rat_mult_poly_l: "c is_polyexp ==> c * (a / b) = c * a /  b" and

   rat_mult_poly_r: "c is_polyexp ==> (a / b) * c = a * c /  b" and

 (*real_times_divide1_eq .. Isa02*) 

   real_times_divide_1_eq:  "-1 * (c / d) = -1 * c / d " and

   real_times_divide_num:   "a is_const ==> a * (c / d) = a * c / d " and

   real_mult_div_cancel2:   "k ~= 0 ==> m * k / (n * k) = m / n" and

 (*real_mult_div_cancel1:   "k ~= 0 ==> k * m / (k * n) = m / n"..Isa02*)

   real_divide_divide1:     "y ~= 0 ==> (u / v) / (y / z) = (u / v) * (z / y)" and

   real_divide_divide1_mg:  "y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)" and

 (*real_divide_divide2_eq:  "x / y / z = x / (y * z)"..Isa02*)

   rat_power:               "(a / b)^^^n = (a^^^n) / (b^^^n)" and

   rat_add:         "[| a is_const; b is_const; c is_const; d is_const |] ==> 

 	           a / c + b / d = (a * d + b * c) / (c * d)" and

   rat_add_assoc:   "[| a is_const; b is_const; c is_const; d is_const |] ==> 

 	           a / c +(b / d + e) = (a * d + b * c)/(d * c) + e" and

   rat_add1:        "[| a is_const; b is_const; c is_const |] ==> 

 	           a / c + b / c = (a + b) / c" and

   rat_add1_assoc:   "[| a is_const; b is_const; c is_const |] ==> 

 	           a / c + (b / c + e) = (a + b) / c + e" and

   rat_add2:        "[| a is_const; b is_const; c is_const |] ==> 

 	           a / c + b = (a + b * c) / c" and

   rat_add2_assoc:  "[| a is_const; b is_const; c is_const |] ==> 

 	           a / c + (b + e) = (a + b * c) / c + e" and

   rat_add3:        "[| a is_const; b is_const; c is_const |] ==> 

 	           a + b / c = (a * c + b) / c" and

   rat_add3_assoc:   "[| a is_const; b is_const; c is_const |] ==> 

 	           a + (b / c + e) = (a * c + b) / c + e"

 section {* Cancellation and addition of fractions *}

 subsection {* Conversion term <--> poly *}

 subsubsection {* Convert a term to the internal representation of a multivariate polynomial *}

 ML {*

 fun monom_of_term  vs (c, es) (Free (id, _)) =

     if is_numeral id 

     then (id |> int_of_str |> the |> curry op * c, es) (*several numerals in one monom*)

     else (c, list_update es (find_index (curry op = id) vs) 1)

   | monom_of_term  vs (c, es) (Const ("Atools.pow", _) $ Free (id, _) $ Free (e, _)) =

     (c, list_update es (find_index (curry op = id) vs) (the (int_of_str e)))

   | monom_of_term vs (c, es) (Const ("Groups.times_class.times", _) $ m1 $ m2) =

     let val (c', es') = monom_of_term vs (c, es) m1

     in monom_of_term vs (c', es') m2 end

   | monom_of_term _ _ t = raise ERROR ("poly malformed with " ^ term2str t)

 fun monoms_of_term vs (t as Free _) =

     [monom_of_term  vs (1, replicate (length vs) 0) t]

   | monoms_of_term vs (t as Const ("Atools.pow", _) $ _ $  _) =

     [monom_of_term  vs (1, replicate (length vs) 0) t]

   | monoms_of_term vs (t as Const ("Groups.times_class.times", _) $ _ $  _) =

     [monom_of_term  vs (1, replicate (length vs) 0) t]

   | monoms_of_term vs (Const ("Groups.plus_class.plus", _) $ ms1 $ ms2) =

     (monoms_of_term vs ms1) @ (monoms_of_term vs ms2)

   | monoms_of_term _ t = raise ERROR ("poly malformed with " ^ term2str t)

 (* convert a term to the internal representation of a multivariate polynomial;

   the conversion is quite liberal, see test --- fun poly_of_term ---:

 * the order of variables and the parentheses within a monomial are arbitrary

 * the coefficient may be somewhere

 * he order and the parentheses within monomials are arbitrary

 But the term must be completely expand + over * (laws of distributivity are not applicable).

 The function requires the free variables as strings already given, 

 because the gcd involves 2 polynomials (with the same length for their list of exponents).

 *)

 fun poly_of_term vs (t as Const ("Groups.plus_class.plus", _) $ _ $ _) =

     (SOME (t |> monoms_of_term vs |> order)

       handle ERROR _ => NONE)

   | poly_of_term vs t =

     (SOME [monom_of_term vs (1, replicate (length vs) 0) t]

       handle ERROR _ => NONE)

 fun is_poly t =

   let 

     val vs = t |> vars |> map str_of_free_opt (* tolerate Var in simplification *)

       |> filter is_some |> map the |> sort string_ord

   in 

     case poly_of_term vs t of SOME _ => true | NONE => false

   end

 val is_expanded = is_poly   (* TODO: check names *)

 val is_polynomial = is_poly (* TODO: check names *)

 *}

 subsubsection {* Convert internal representation of a multivariate polynomial to a term *}

 ML {*

 fun term_of_es _ _ _ [] = [] (*assumes same length for vs and es*)

   | term_of_es baseT expT (_ :: vs) (0 :: es) =

     [] @ term_of_es baseT expT vs es

   | term_of_es baseT expT (v :: vs) (1 :: es) =

     [(Free (v, baseT))] @ term_of_es baseT expT vs es

   | term_of_es baseT expT (v :: vs) (e :: es) =

     [Const ("Atools.pow", [baseT, expT] ---> baseT) $ 

       (Free (v, baseT)) $  (Free (isastr_of_int e, expT))]

     @ term_of_es baseT expT vs es

 fun term_of_monom baseT expT vs ((c, es): monom) =

     let val es' = term_of_es baseT expT vs es

     in 

       if c = 1 

       then 

         if es' = [] (*if es = [0,0,0,...]*)

         then Free (isastr_of_int c, baseT)

         else foldl (HOLogic.mk_binop "Groups.times_class.times") (hd es', tl es')

       else foldl (HOLogic.mk_binop "Groups.times_class.times") (Free (isastr_of_int c, baseT), es') 

     end

 fun term_of_poly baseT expT vs p =

   let val monos = map (term_of_monom baseT expT vs) p

   in foldl (HOLogic.mk_binop "Groups.plus_class.plus") (hd monos, tl monos) end

 *}

 subsection {* Apply gcd_poly for cancelling and adding fractions as terms *}

 ML {*

 fun mk_noteq_0 baseT t = 

   Const ("HOL.Not", HOLogic.boolT --> HOLogic.boolT) $ 

     (Const ("HOL.eq", [baseT, baseT] ---> HOLogic.boolT) $ t $ Free ("0", HOLogic.realT))

 fun mk_asms baseT ts =

   let val as' = filter_out is_num ts (* asm like "2 ~= 0" is needless *)

   in map (mk_noteq_0 baseT) as' end

 *}

 subsubsection {* Factor out gcd for cancellation *}

 ML {*

 fun check_fraction t =

   let val Const ("Fields.inverse_class.divide", _) $ numerator $ denominator = t

   in SOME (numerator, denominator) end

   handle Bind => NONE

 (* prepare a term for cancellation by factoring out the gcd

   assumes: is a fraction with outmost "/"*)

 fun factout_p_ (thy: theory) t =

   let val opt = check_fraction t

   in

     case opt of 

       NONE => NONE

     | SOME (numerator, denominator) =>

       let 

         val vs = t |> vars |> map str_of_free_opt (* tolerate Var in simplification *)

           |> filter is_some |> map the |> sort string_ord

         val baseT = type_of numerator

         val expT = HOLogic.realT

       in

         case (poly_of_term vs numerator, poly_of_term vs denominator) of

           (SOME a, SOME b) =>

             let

               val ((a', b'), c) = gcd_poly a b

               val es = replicate (length vs) 0 

             in

               if c = [(1, es)] orelse c = [(~1, es)]

               then NONE

               else 

                 let

                   val b't = term_of_poly baseT expT vs b'

                   val ct = term_of_poly baseT expT vs c

                   val t' = 

                     HOLogic.mk_binop "Fields.inverse_class.divide" 

                       (HOLogic.mk_binop "Groups.times_class.times"

                         (term_of_poly baseT expT vs a', ct),

                        HOLogic.mk_binop "Groups.times_class.times" (b't, ct))

                 in SOME (t', mk_asms baseT [b't, ct]) end

             end

         | _ => NONE : (term * term list) option

       end

   end

 *}

 subsubsection {* Cancel a fraction *}

 ML {*

 (* cancel a term by the gcd ("" denote terms with internal algebraic structure)

   cancel_p_ :: theory \<Rightarrow> term  \<Rightarrow> (term \<times> term list) option

   cancel_p_ thy "a / b" = SOME ("a' / b'", ["b' \<noteq> 0"])

   assumes: a is_polynomial  \<and>  b is_polynomial  \<and>  b \<noteq> 0

   yields

     SOME ("a' / b'", ["b' \<noteq> 0"]). gcd_poly a b \<noteq> 1  \<and>  gcd_poly a b \<noteq> -1  \<and>  

       a' * gcd_poly a b = a  \<and>  b' * gcd_poly a b = b

     \<or> NONE *)

 fun cancel_p_ (_: theory) t =

   let val opt = check_fraction t

   in

     case opt of 

       NONE => NONE

     | SOME (numerator, denominator) =>

       let 

         val vs = t |> vars |> map str_of_free_opt (* tolerate Var in simplification *)

           |> filter is_some |> map the |> sort string_ord

         val baseT = type_of numerator

         val expT = HOLogic.realT

       in

         case (poly_of_term vs numerator, poly_of_term vs denominator) of

           (SOME a, SOME b) =>

             let

               val ((a', b'), c) = gcd_poly a b

               val es = replicate (length vs) 0 

             in

               if c = [(1, es)] orelse c = [(~1, es)]

               then NONE

               else 

                 let

                   val bt' = term_of_poly baseT expT vs b'

                   val ct = term_of_poly baseT expT vs c

                   val t' = 

                     HOLogic.mk_binop "Fields.inverse_class.divide" 

                       (term_of_poly baseT expT vs a', bt')

                   val asm = mk_asms baseT [bt']

                 in SOME (t', asm) end

             end

         | _ => NONE : (term * term list) option

       end

   end

 *}

 subsubsection {* Factor out to a common denominator for addition *}

 ML {*

 (* addition of fractions allows (at most) one non-fraction (a monomial) *)

 fun check_frac_sum 

     (Const ("Groups.plus_class.plus", _) $ 

       (Const ("Fields.inverse_class.divide", _) $ n1 $ d1) $

       (Const ("Fields.inverse_class.divide", _) $ n2 $ d2))

     = SOME ((n1, d1), (n2, d2))

   | check_frac_sum 

     (Const ("Groups.plus_class.plus", _) $ 

       nofrac $ 

       (Const ("Fields.inverse_class.divide", _) $ n2 $ d2))

     = SOME ((nofrac, Free ("1", HOLogic.realT)), (n2, d2))

   | check_frac_sum 

     (Const ("Groups.plus_class.plus", _) $ 

       (Const ("Fields.inverse_class.divide", _) $ n1 $ d1) $ 

       nofrac)

     = SOME ((n1, d1), (nofrac, Free ("1", HOLogic.realT)))

   | check_frac_sum _ = NONE  

 (* prepare a term for addition by providing the least common denominator as a product

   assumes: is a term with outmost "+" and at least one outmost "/" in respective summands*)

 fun common_nominator_p_ (_: theory) t =

   let val opt = check_frac_sum t

   in

     case opt of 

       NONE => NONE

     | SOME ((n1, d1), (n2, d2)) =>

       let 

         val vs = t |> vars |> map str_of_free_opt (* tolerate Var in simplification *)

           |> filter is_some |> map the |> sort string_ord

       in

         case (poly_of_term vs d1, poly_of_term vs d2) of

           (SOME a, SOME b) =>

             let

               val ((a', b'), c) = gcd_poly a b

               val (baseT, expT) = (type_of n1, HOLogic.realT)

               val [d1', d2', c'] = map (term_of_poly baseT expT vs) [a', b', c]

               (*----- minimum of parentheses & nice result, but breaks tests: -------------

               val denom = HOLogic.mk_binop "Groups.times_class.times" 

                 (HOLogic.mk_binop "Groups.times_class.times" (d1', d2'), c') -------------*)

               val denom =

                 if c = [(1, replicate (length vs) 0)]

                 then HOLogic.mk_binop "Groups.times_class.times" (d1', d2')

                 else

                   HOLogic.mk_binop "Groups.times_class.times" (c',

                   HOLogic.mk_binop "Groups.times_class.times" (d1', d2')) (*--------------*)

               val t' =

                 HOLogic.mk_binop "Groups.plus_class.plus"

                   (HOLogic.mk_binop "Fields.inverse_class.divide"

                     (HOLogic.mk_binop "Groups.times_class.times" (n1, d2'), denom),

                   HOLogic.mk_binop "Fields.inverse_class.divide" 

                     (HOLogic.mk_binop "Groups.times_class.times" (n2, d1'), denom))

               val asm = mk_asms baseT [d1', d2', c']

             in SOME (t', asm) end

         | _ => NONE : (term * term list) option

       end

   end

 *}

 subsubsection {* Addition of at least one fraction within a sum *}

 ML {*

 (* add fractions

   assumes: is a term with outmost "+" and at least one outmost "/" in respective summands

   NOTE: the case "(_ + _) + _" need not be considered due to iterated addition.*)

 fun add_fraction_p_ (_: theory) t =

   case check_frac_sum t of 

     NONE => NONE

   | SOME ((n1, d1), (n2, d2)) =>

     let 

       val vs = t |> vars |> map str_of_free_opt (* tolerate Var in simplification *)

         |> filter is_some |> map the |> sort string_ord

     in

       case (poly_of_term vs n1, poly_of_term vs d1, poly_of_term vs n2, poly_of_term vs d2) of

         (SOME _, SOME a, SOME _, SOME b) =>

           let

             val ((a', b'), c) = gcd_poly a b

             val (baseT, expT) = (type_of n1, HOLogic.realT)

             val nomin = term_of_poly baseT expT vs 

               (((the (poly_of_term vs n1)) %%*%% b') %%+%% ((the (poly_of_term vs n2)) %%*%% a')) 

             val denom = term_of_poly baseT expT vs ((c %%*%% a') %%*%% b')

             val t' = HOLogic.mk_binop "Fields.inverse_class.divide" (nomin, denom)

           in SOME (t', mk_asms baseT [denom]) end

       | _ => NONE : (term * term list) option

     end

 *}

 section {* Embed cancellation and addition into rewriting *}

 ML {* val thy = @{theory} *}

 subsection {* Rulesets and predicate for embedding *}

 ML {*

 (* evaluates conditions in calculate_Rational *)

 val calc_rat_erls =

   prep_rls

     (Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

       erls = e_rls, srls = Erls, calc = [], errpatts = [],

       rules = 

         [Calc ("HOL.eq", eval_equal "#equal_"),

         Calc ("Atools.is'_const", eval_const "#is_const_"),

         Thm ("not_true", num_str @{thm not_true}),

         Thm ("not_false", num_str @{thm not_false})], 

       scr = EmptyScr});

 (* simplifies expressions with numerals;

    does NOT rearrange the term by AC-rewriting; thus terms with variables 

    need to have constants to be commuted together respectively           *)

 val calculate_Rational =

   prep_rls (merge_rls "calculate_Rational"

     (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

       erls = calc_rat_erls, srls = Erls,

       calc = [], errpatts = [],

       rules = 

         [Calc ("Fields.inverse_class.divide", eval_cancel "#divide_e"),

         Thm ("minus_divide_left", num_str (@{thm minus_divide_left} RS @{thm sym})),

           (*SYM - ?x / ?y = - (?x / ?y)  may come from subst*)

         Thm ("rat_add", num_str @{thm rat_add}),

           (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \

           \a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)

         Thm ("rat_add1", num_str @{thm rat_add1}),

           (*"[| a is_const; b is_const; c is_const |] ==> a / c + b / c = (a + b) / c"*)

         Thm ("rat_add2", num_str @{thm rat_add2}),

           (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> ?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)

         Thm ("rat_add3", num_str @{thm rat_add3}),

           (*"[| a is_const; b is_const; c is_const |] ==> a + b / c = (a * c) / c + b / c"\

           .... is_const to be omitted here FIXME*)

         Thm ("rat_mult", num_str @{thm rat_mult}), 

           (*a / b * (c / d) = a * c / (b * d)*)

         Thm ("times_divide_eq_right", num_str @{thm times_divide_eq_right}),

           (*?x * (?y / ?z) = ?x * ?y / ?z*)

         Thm ("times_divide_eq_left", num_str @{thm times_divide_eq_left}),

           (*?y / ?z * ?x = ?y * ?x / ?z*)

         Thm ("real_divide_divide1", num_str @{thm real_divide_divide1}),

           (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)

         Thm ("divide_divide_eq_left", num_str @{thm divide_divide_eq_left}),

           (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

         Thm ("rat_power", num_str @{thm rat_power}),

           (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

         Thm ("mult_cross", num_str @{thm mult_cross}),

           (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)

         Thm ("mult_cross1", num_str @{thm mult_cross1}),

           (*"   b ~= 0            ==> (a / b = c    ) = (a     = b * c)*)

         Thm ("mult_cross2", num_str @{thm mult_cross2})

           (*"           d ~= 0    ==> (a     = c / d) = (a * d =     c)*)], 

       scr = EmptyScr})

     calculate_Poly);

 (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)

 fun eval_is_expanded (thmid:string) _ 

 		       (t as (Const("Rational.is'_expanded", _) $ arg)) thy = 

     if is_expanded arg

     then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "", 

 	         Trueprop $ (mk_equality (t, @{term True})))

     else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "", 

 	         Trueprop $ (mk_equality (t, @{term False})))

   | eval_is_expanded _ _ _ _ = NONE; 

 calclist':= overwritel (!calclist', 

    [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))]);

 val rational_erls = 

   merge_rls "rational_erls" calculate_Rational 

     (append_rls "is_expanded" Atools_erls 

       [Calc ("Rational.is'_expanded", eval_is_expanded "")]);

 *}

 subsection {* Embed cancellation into rewriting *}

 ML {*

 local (* cancel_p *)

 val {rules = rules, rew_ord = (_, ro), ...} = rep_rls (assoc_rls "rev_rew_p");

 fun init_state thy eval_rls ro t =

   let

     val SOME (t', _) = factout_p_ thy t;

     val SOME (t'', asm) = cancel_p_ thy t;

     val der = reverse_deriv thy eval_rls rules ro NONE t';

     val der = der @ 

       [(Thm ("real_mult_div_cancel2", num_str @{thm real_mult_div_cancel2}), (t'', asm))]

     val rs = (distinct_Thm o (map #1)) der

   	val rs = filter_out (eq_Thms 

   	  ["sym_real_add_zero_left", "sym_real_mult_0", "sym_real_mult_1"]) rs

   in (t, t'', [rs(*one in order to ease locate_rule*)], der) end;

 fun locate_rule thy eval_rls ro [rs] t r =

     if member op = ((map (id_of_thm)) rs) (id_of_thm r)

     then 

       let val ropt = rewrite_ thy ro eval_rls true (thm_of_thm r) t;

       in

         case ropt of SOME ta => [(r, ta)]

 	      | NONE => (tracing 

 	          ("### locate_rule:  rewrite " ^ id_of_thm r ^ " " ^ term2str t ^ " = NONE"); []) 

 			end

     else (tracing ("### locate_rule:  " ^ id_of_thm r ^ " not mem rrls"); [])

   | locate_rule _ _ _ _ _ _ = error "locate_rule: doesnt match rev-sets in istate";

 fun next_rule thy eval_rls ro [rs] t =

     let

       val der = make_deriv thy eval_rls rs ro NONE t;

     in case der of (_, r, _) :: _ => SOME r | _ => NONE end

   | next_rule _ _ _ _ _ = error ("next_rule: doesnt match rev-sets in istate");

 fun attach_form (_:rule list list) (_:term) (_:term) = 

   [(*TODO*)]: (rule * (term * term list)) list;

 in

 val cancel_p = 

   Rrls {id = "cancel_p", prepat = [],

 	rew_ord=("ord_make_polynomial", ord_make_polynomial false thy),

 	erls = rational_erls, 

 	calc = 

 	  [("PLUS", ("Groups.plus_class.plus", eval_binop "#add_")),

 	  ("TIMES" , ("Groups.times_class.times", eval_binop "#mult_")),

 	  ("DIVIDE", ("Fields.inverse_class.divide", eval_cancel "#divide_e")),

 	  ("POWER", ("Atools.pow", eval_binop "#power_"))],

     errpatts = [],

 	scr =

 	  Rfuns {init_state  = init_state thy Atools_erls ro,

 		normal_form = cancel_p_ thy, 

 		locate_rule = locate_rule thy Atools_erls ro,

 		next_rule   = next_rule thy Atools_erls ro,

 		attach_form = attach_form}}

 end; (* local cancel_p *)

 *}

 subsection {* Embed addition into rewriting *}

 ML {*

 local (* add_fractions_p *)

 val {rules = rules, rew_ord = (_, ro), ...} = rep_rls (assoc_rls "make_polynomial");

 val {rules, rew_ord=(_,ro),...} = rep_rls (assoc_rls "rev_rew_p");

 fun init_state thy eval_rls ro t =

   let 

     val SOME (t',_) = common_nominator_p_ thy t;

     val SOME (t'', asm) = add_fraction_p_ thy t;

     val der = reverse_deriv thy eval_rls rules ro NONE t';

     val der = der @ 

       [(Thm ("real_mult_div_cancel2", num_str @{thm real_mult_div_cancel2}), (t'',asm))]

     val rs = (distinct_Thm o (map #1)) der;

     val rs = filter_out (eq_Thms 

       ["sym_real_add_zero_left", "sym_real_mult_0", "sym_real_mult_1"]) rs;

   in (t, t'', [rs(*here only _ONE_*)], der) end;

 fun locate_rule thy eval_rls ro [rs] t r =

     if member op = ((map (id_of_thm)) rs) (id_of_thm r)

     then 

       let val ropt = rewrite_ thy ro eval_rls true (thm_of_thm r) t;

       in 

         case ropt of

           SOME ta => [(r, ta)]

 	      | NONE => 

 	        (tracing ("### locate_rule:  rewrite " ^ id_of_thm r ^ " " ^ term2str t ^ " = NONE");

 	        []) end

     else (tracing ("### locate_rule:  " ^ id_of_thm r ^ " not mem rrls"); [])

   | locate_rule _ _ _ _ _ _ = error "locate_rule: doesnt match rev-sets in istate";

 fun next_rule thy eval_rls ro [rs] t =

     let val der = make_deriv thy eval_rls rs ro NONE t;

     in 

       case der of

 	      (_,r,_)::_ => SOME r

 	    | _ => NONE

     end

   | next_rule _ _ _ _ _ = error ("next_rule: doesnt match rev-sets in istate");

 val pat0 = parse_patt thy "?r/?s+?u/?v :: real";

 val pat1 = parse_patt thy "?r/?s+?u    :: real";

 val pat2 = parse_patt thy "?r   +?u/?v :: real";

 val prepat = [([@{term True}], pat0),

 	      ([@{term True}], pat1),

 	      ([@{term True}], pat2)];

 in

 val add_fractions_p =

   Rrls {id = "add_fractions_p", prepat=prepat,

     rew_ord = ("ord_make_polynomial", ord_make_polynomial false thy),

     erls = rational_erls,

     calc = [("PLUS", ("Groups.plus_class.plus", eval_binop "#add_")),

       ("TIMES", ("Groups.times_class.times", eval_binop "#mult_")),

       ("DIVIDE", ("Fields.inverse_class.divide", eval_cancel "#divide_e")),

       ("POWER", ("Atools.pow", eval_binop "#power_"))],

     errpatts = [],

     scr = Rfuns {init_state  = init_state thy Atools_erls ro,

       normal_form = add_fraction_p_ thy,

       locate_rule = locate_rule thy Atools_erls ro,

       next_rule   = next_rule thy Atools_erls ro,

       attach_form = attach_form}}

 end; (*local add_fractions_p *)

 *}

 subsection {* Cancelling and adding all occurrences in a term /////////////////////////////*}

 ML {*

 (*copying cancel_p_rls + add her caused error in interface.sml*)

 *}

 section {* Rulesets for general simplification *}

 ML {* 

 (*-------------------18.3.03 --> struct <-----------vvv--

 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)

  -------------------18.3.03 --> struct <-----------vvv--*)

 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)

 val powers_erls = prep_rls(

   Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

       erls = e_rls, srls = Erls, calc = [], errpatts = [],

       rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),

 	       Calc ("Atools.is'_even",eval_is_even "#is_even_"),

 	       Calc ("Orderings.ord_class.less",eval_equ "#less_"),

 	       Thm ("not_false", num_str @{thm not_false}),

 	       Thm ("not_true", num_str @{thm not_true}),

 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

 	       ],

       scr = EmptyScr

       }:rls);

 (*.all powers over + distributed; atoms over * collected, other distributed

    contains absolute minimum of thms for context in norm_Rational .*)

 val powers = prep_rls(

   Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

       erls = powers_erls, srls = Erls, calc = [], errpatts = [],

       rules = [Thm ("realpow_multI", num_str @{thm realpow_multI}),

 	       (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)

 	       Thm ("realpow_pow",num_str @{thm realpow_pow}),

 	       (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)

 	       Thm ("realpow_oneI",num_str @{thm realpow_oneI}),

 	       (*"r ^^^ 1 = r"*)

 	       Thm ("realpow_minus_even",num_str @{thm realpow_minus_even}),

 	       (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)

 	       Thm ("realpow_minus_odd",num_str @{thm realpow_minus_odd}),

 	       (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)

 	       (*----- collect atoms over * -----*)

 	       Thm ("realpow_two_atom",num_str @{thm realpow_two_atom}),	

 	       (*"r is_atom ==> r * r = r ^^^ 2"*)

 	       Thm ("realpow_plus_1",num_str @{thm realpow_plus_1}),		

 	       (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)

 	       Thm ("realpow_addI_atom",num_str @{thm realpow_addI_atom}),

 	       (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)

 	       (*----- distribute none-atoms -----*)

 	       Thm ("realpow_def_atom",num_str @{thm realpow_def_atom}),

 	       (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)

 	       Thm ("realpow_eq_oneI",num_str @{thm realpow_eq_oneI}),

 	       (*"1 ^^^ n = 1"*)

 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

 	       ],

       scr = EmptyScr

       }:rls);

 (*.contains absolute minimum of thms for context in norm_Rational.*)

 val rat_mult_divide = prep_rls(

   Rls {id = "rat_mult_divide", preconds = [], 

       rew_ord = ("dummy_ord", dummy_ord), 

       erls = e_rls, srls = Erls, calc = [], errpatts = [],

       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),

 	       (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

 	       (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

 	       otherwise inv.to a / b / c = ...*)

 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

 	       (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much

 		     and does not commute a / b * c ^^^ 2 !*)

 	       Thm ("divide_divide_eq_right", 

                      num_str @{thm divide_divide_eq_right}),

 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

 	       Thm ("divide_divide_eq_left",

                      num_str @{thm divide_divide_eq_left}),

 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")

 	       ],

       scr = EmptyScr

       }:rls);

 (*.contains absolute minimum of thms for context in norm_Rational.*)

 val reduce_0_1_2 = prep_rls(

   Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

       erls = e_rls, srls = Erls, calc = [], errpatts = [],

       rules = [(*Thm ("divide_1",num_str @{thm divide_1}),

 		 "?x / 1 = ?x" unnecess.for normalform*)

 	       Thm ("mult_1_left",num_str @{thm mult_1_left}),                 

 	       (*"1 * z = z"*)

 	       (*Thm ("real_mult_minus1",num_str @{thm real_mult_minus1}),

 	       "-1 * z = - z"*)

 	       (*Thm ("real_minus_mult_cancel",num_str @{thm real_minus_mult_cancel}),

 	       "- ?x * - ?y = ?x * ?y"*)

 	       Thm ("mult_zero_left",num_str @{thm mult_zero_left}),        

 	       (*"0 * z = 0"*)

 	       Thm ("add_0_left",num_str @{thm add_0_left}),

 	       (*"0 + z = z"*)

 	       (*Thm ("right_minus",num_str @{thm right_minus}),

 	       "?z + - ?z = 0"*)

 	       Thm ("sym_real_mult_2",

                      num_str (@{thm real_mult_2} RS @{thm sym})),	

 	       (*"z1 + z1 = 2 * z1"*)

 	       Thm ("real_mult_2_assoc",num_str @{thm real_mult_2_assoc}),

 	       (*"z1 + (z1 + k) = 2 * z1 + k"*)

 	       Thm ("divide_zero_left",num_str @{thm divide_zero_left})

 	       (*"0 / ?x = 0"*)

 	       ], scr = EmptyScr}:rls);

 (*erls for calculate_Rational; 

   make local with FIXX@ME result:term *term list WN0609???SKMG*)

 val norm_rat_erls = prep_rls(

   Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

       erls = e_rls, srls = Erls, calc = [], errpatts = [],

       rules = [Calc ("Atools.is'_const",eval_const "#is_const_")

 	       ], scr = EmptyScr}:rls);

 (* consists of rls containing the absolute minimum of thms *)

 (*040209: this version has been used by RL for his equations,

 which is now replaced by MGs version "norm_Rational" below *)

 val norm_Rational_min = prep_rls(

   Rls {id = "norm_Rational_min", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

       erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

       rules = [(*sequence given by operator precedence*)

 	       Rls_ discard_minus,

 	       Rls_ powers,

 	       Rls_ rat_mult_divide,

 	       Rls_ expand,

 	       Rls_ reduce_0_1_2,

 	       Rls_ order_add_mult,

 	       Rls_ collect_numerals,

 	       Rls_ add_fractions_p,

 	       Rls_ cancel_p

 	       ],

       scr = EmptyScr}:rls);

 val norm_Rational_parenthesized = prep_rls(

   Seq {id = "norm_Rational_parenthesized", preconds = []:term list, 

        rew_ord = ("dummy_ord", dummy_ord),

       erls = Atools_erls, srls = Erls,

       calc = [], errpatts = [],

       rules = [Rls_  norm_Rational_min,

 	       Rls_ discard_parentheses

 	       ],

       scr = EmptyScr}:rls);      

 (*WN030318???SK: simplifies all but cancel and common_nominator*)

 val simplify_rational = 

     merge_rls "simplify_rational" expand_binoms

     (append_rls "divide" calculate_Rational

 		[Thm ("divide_1",num_str @{thm divide_1}),

 		 (*"?x / 1 = ?x"*)

 		 Thm ("rat_mult",num_str @{thm rat_mult}),

 		 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

 		 Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

 		 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

 		 otherwise inv.to a / b / c = ...*)

 		 Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

 		 (*"?a / ?b * ?c = ?a * ?c / ?b"*)

 		 Thm ("add_minus",num_str @{thm add_minus}),

 		 (*"?a + ?b - ?b = ?a"*)

 		 Thm ("add_minus1",num_str @{thm add_minus1}),

 		 (*"?a - ?b + ?b = ?a"*)

 		 Thm ("divide_minus1",num_str @{thm divide_minus1})

 		 (*"?x / -1 = - ?x"*)

 		 ]);

 *}

 ML {* 

 val add_fractions_p_rls = prep_rls(

   Rls {id = "add_fractions_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

 	  erls = e_rls, srls = Erls, calc = [], errpatts = [],

 	  rules = [Rls_ add_fractions_p], 

 	  scr = EmptyScr});

 (* "Rls" causes repeated application of cancel_p to one and the same term *)

 val cancel_p_rls = prep_rls(

   Rls 

     {id = "cancel_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

     erls = e_rls, srls = Erls, calc = [], errpatts = [],

     rules = [Rls_ cancel_p], 

 	  scr = EmptyScr});

 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

     used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)

 val rat_mult_poly = prep_rls(

   Rls {id = "rat_mult_poly", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

 	  erls = append_rls "e_rls-is_polyexp" e_rls [Calc ("Poly.is'_polyexp", eval_is_polyexp "")], 

 	  srls = Erls, calc = [], errpatts = [],

 	  rules = 

 	    [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

 	    (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

 	    Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})

 	    (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) ], 

 	  scr = EmptyScr});

 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

     used in looping part norm_Rational_rls, see example DA-M02-main.p.60 

     .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls, 

     I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028 

     ... WN0609???MG.*)

 val rat_mult_div_pow = prep_rls(

   Rls {id = "rat_mult_div_pow", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

     erls = e_rls, srls = Erls, calc = [], errpatts = [],

     rules = [Thm ("rat_mult", num_str @{thm rat_mult}),

       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

       Thm ("rat_mult_poly_l", num_str @{thm rat_mult_poly_l}),

       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

       Thm ("rat_mult_poly_r", num_str @{thm rat_mult_poly_r}),

       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

       Thm ("real_divide_divide1_mg", num_str @{thm real_divide_divide1_mg}),

       (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)

       Thm ("divide_divide_eq_right", num_str @{thm divide_divide_eq_right}),

       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

       Thm ("divide_divide_eq_left", num_str @{thm divide_divide_eq_left}),

       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

       Calc ("Fields.inverse_class.divide", eval_cancel "#divide_e"),

       Thm ("rat_power", num_str @{thm rat_power})

       (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

       ],

     scr = EmptyScr}:rls);

 val rat_reduce_1 = prep_rls(

   Rls {id = "rat_reduce_1", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

     erls = e_rls, srls = Erls, calc = [], errpatts = [], 

     rules = 

       [Thm ("divide_1", num_str @{thm divide_1}),

       (*"?x / 1 = ?x"*)

       Thm ("mult_1_left", num_str @{thm mult_1_left})           

       (*"1 * z = z"*)

       ],

     scr = EmptyScr}:rls);

 (* looping part of norm_Rational *)

 val norm_Rational_rls = prep_rls (

   Rls {id = "norm_Rational_rls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 

     erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

     rules = [Rls_ add_fractions_p_rls,

       Rls_ rat_mult_div_pow,

       Rls_ make_rat_poly_with_parentheses,

       Rls_ cancel_p_rls,

       Rls_ rat_reduce_1

       ],

     scr = EmptyScr}:rls);

 val norm_Rational = prep_rls (

   Seq 

     {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord", dummy_ord), 

     erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

     rules = [Rls_ discard_minus,

       Rls_ rat_mult_poly,             (* removes double fractions like a/b/c *)

       Rls_ make_rat_poly_with_parentheses,

       Rls_ cancel_p_rls,

       Rls_ norm_Rational_rls,         (* the main rls, looping (#) *)

       Rls_ discard_parentheses1       (* mult only *)

       ],

     scr = EmptyScr}:rls);

 *}

 ML {*

 ruleset' := overwritelthy @{theory} (!ruleset',

   [("calculate_Rational", calculate_Rational),

    ("calc_rat_erls",calc_rat_erls),

    ("rational_erls", rational_erls),

    ("cancel_p", cancel_p),

    ("add_fractions_p", add_fractions_p),

    ("add_fractions_p_rls", add_fractions_p_rls),

    (*WN120410 ("discard_minus", discard_minus), is defined in Poly.thy*)

    ("powers_erls", powers_erls),

    ("powers", powers),

    ("rat_mult_divide", rat_mult_divide),

    ("reduce_0_1_2", reduce_0_1_2),

    ("rat_reduce_1", rat_reduce_1),

    ("norm_rat_erls", norm_rat_erls),

    ("norm_Rational", norm_Rational),

    ("norm_Rational_rls", norm_Rational_rls),

    ("norm_Rational_parenthesized", norm_Rational_parenthesized),

    ("rat_mult_poly", rat_mult_poly),

    ("rat_mult_div_pow", rat_mult_div_pow),

    ("cancel_p_rls", cancel_p_rls)

    ]);

 *}

 section {* A problem for simplification of rationals *}

 ML {*

 store_pbt

  (prep_pbt thy "pbl_simp_rat" [] e_pblID

  (["rational","simplification"],

   [("#Given" ,["Term t_t"]),

    ("#Where" ,["t_t is_ratpolyexp"]),

    ("#Find"  ,["normalform n_n"])

   ],

   append_rls "e_rls" e_rls [(*for preds in where_*)], 

   SOME "Simplify t_t", 

   [["simplification","of_rationals"]]));

 *}

 section {* A methods for simplification of rationals *}

 ML {*

 (*WN061025 this methods script is copied from (auto-generated) script

   of norm_Rational in order to ease repair on inform*)

 store_met (prep_met thy "met_simp_rat" [] e_metID (["simplification","of_rationals"],

 	  [("#Given" ,["Term t_t"]),

 	   ("#Where" ,["t_t is_ratpolyexp"]),

 	   ("#Find"  ,["normalform n_n"])],

 	  {rew_ord'="tless_true", rls' = e_rls, calc = [], srls = e_rls, 

 	   prls = append_rls "simplification_of_rationals_prls" e_rls 

 				    [(*for preds in where_*) Calc ("Rational.is'_ratpolyexp", eval_is_ratpolyexp "")],

 				   crls = e_rls, errpats = [],

    nrls = norm_Rational_rls},

    "Script SimplifyScript (t_t::real) =                             " ^

    "  ((Try (Rewrite_Set discard_minus False) @@                   " ^

    "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^

    "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^

    "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^

    "    (Repeat                                                     " ^

    "     ((Try (Rewrite_Set add_fractions_p_rls False) @@        " ^

    "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^

    "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^

    "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^

    "       Try (Rewrite_Set rat_reduce_1 False)))) @@               " ^

    "    Try (Rewrite_Set discard_parentheses1 False))               " ^

    "   t_t)"));

 *}

 end