wneuper/isa: src/Tools/isac/Knowledge/Rational.thy@fb5a2c25e3f0

     1 (* rationals, i.e. fractions of multivariate polynomials over the real field

     2    author: isac team

     3    Copyright (c) isac team 2002

     4    Use is subject to license terms.

     6    depends on Poly (and not on Atools), because

     7    fractions with _normalized_ polynomials are canceled, added, etc.

     9    ATTENTION WN130616: WITH Unsynchronized.ref Rational.thy TAKES ~1MIN FOR EVALUATION

    10 *)

    12 theory Rational

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

    14 begin

    16 consts

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

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

    20   get_denominator :: "real => real"

    21   get_numerator   :: "real => real"

    24 axioms(*axiomatization where*) (*.not contained in Isabelle2002,

    25           stated as axioms, TODO?: prove as theorems*)

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

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

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

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

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

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

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

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

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

    40 (*real_times_divide1_eq .. Isa02*)

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

    42   real_times_divide_num:   "a is_const ==>

    43 	          	   a     * (c / d) = a * c /      d " (*and*)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    72 section {* Cancellation and addition of fractions *}

    73 subsection {* Auxiliary functions and data *}

    74 subsubsection {* Conversion term <--> poly *}

    75 ML {*

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

    77     if is_numeral id

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

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

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

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

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

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

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

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

    87 fun monoms_of_term vs (t as Free _) =

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

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

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

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

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

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

    94     (monoms_of_term vs ms1) @ (monoms_of_term vs ms2)

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

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

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

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

   100 * the coefficient may be somewhere

   101 * he order and the parentheses within monomials are arbitrary

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

   104 The function requires the free variables as strings already given,

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

   106 *)

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

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

   109       handle ERROR _ => NONE)

   110   | poly_of_term vs t =

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

   112       handle ERROR _ => NONE)

   114 fun is_poly t =

   115   let val vs = t |> vars |> map free2str |> sort string_ord

   116   in

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

   118   end

   119 val is_expanded = is_poly

   121 (* convert internal representation of a multivariate polynomial to a term*)

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

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

   124     [] @ term_of_es baseT expT vs es

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

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

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

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

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

   130     @ term_of_es baseT expT vs es

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

   133     let val es' = term_of_es baseT expT vs es

   134     in

   135       if c = 1

   136       then

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

   138         then Free (isastr_of_int c, baseT)

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

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

   141     end

   143 fun term_of_poly baseT expT vs p =

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

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

   146 *}

   148 text {*calculate in rationals: gcd, lcm, etc.

   149       (c) Stefan Karnel 2002

   150       Institute for Mathematics D and Institute for Software Technology,

   151       TU-Graz SS 2002 *}

   153 text {* Remark on notions in the documentation below:

   154     referring to the remark on 'polynomials' in Poly.sml we use

   155     [2] 'polynomial' normalform (Polynom)

   156     [3] 'expanded_term' normalform (Ausmultiplizierter Term),

   157     where normalform [2] is a special case of [3], i.e. [3] implies [2].

   158     Instead of

   159       'fraction with numerator and nominator both in normalform [2]'

   160       'fraction with numerator and nominator both in normalform [3]'

   161     we say:

   162       'fraction in normalform [2]'

   163       'fraction in normalform [3]'

   164     or

   165       'fraction [2]'

   166       'fraction [3]'.

   167     a 'simple fraction' is a term with '/' as outmost operator and

   168     numerator and nominator in normalform [2] or [3].

   169 *}

   171 ML {*

   172 val thy = @{theory};

   174 signature RATIONALI =

   175 sig

   176   type mv_monom

   177   type mv_poly

   178   val add_fraction_p_ : theory -> term -> (term * term list) option

   179   val calculate_Rational : rls

   180   val calc_rat_erls:rls

   181   val cancel_p : rls

   182   val cancel_p_ : theory -> term -> (term * term list) option

   183   val common_nominator_p : rls

   184   val common_nominator_p_ : theory -> term -> (term * term list) option

   185   val eval_is_expanded : string -> 'a -> term -> theory ->

   186 			 (string * term) option

   187   val expanded2polynomial : term -> term option

   188   val factout_p_ : theory -> term -> (term * term list) option

   189   val is_expanded : term -> bool

   190   val is_polynomial : term -> bool

   192   val mv_gcd : (int * int list) list -> mv_poly -> mv_poly

   193   val mv_lcm : mv_poly -> mv_poly -> mv_poly

   195   val norm_expanded_rat_ : theory -> term -> (term * term list) option

   196 (*WN0602.2.6.pull into struct !!!

   197   val norm_Rational : rls(*.normalizes an arbitrary rational term without

   198                            roots into a simple and canceled fraction

   199                            with normalform [2].*)

   200 *)

   201 (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME

   202       rls               (*.normalizes an rational term [2] without

   203                            roots into a simple and canceled fraction

   204                            with normalform [2].*)

   205 *)

   206   val norm_rational_ : theory -> term -> (term * term list) option

   207   val polynomial2expanded : term -> term option

   208   val rational_erls :

   209       rls             (*.evaluates an arbitrary rational term with numerals.*)

   211 (*WN0210???SK: fehlen Funktionen, die exportiert werden sollen ? *)

   212 end (* sig *)

   214 (*.**************************************************************************

   215 survey on the functions

   216 ~~~~~~~~~~~~~~~~~~~~~~~

   217  [2] 'polynomial'   :rls               | [3]'expanded_term':rls

   218 --------------------:------------------+-------------------:-----------------

   219  factout_p_         :                  | factout_          :

   220  cancel_p_          :                  | cancel_           :

   221                     :cancel_p          |                   :cancel

   222 --------------------:------------------+-------------------:-----------------

   223  common_nominator_p_:                  | common_nominator_ :

   224                     :common_nominator_p|                   :common_nominator

   225  add_fraction_p_    :                  | add_fraction_     :

   226 --------------------:------------------+-------------------:-----------------

   227 ???SK                 :norm_rational_p   |                   :norm_rational

   229 This survey shows only the principal functions for reuse, and the identifiers

   230 of the rls exported. The list below shows some more useful functions.

   233 conversion from Isabelle-term to internal representation

   234 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   236 ... BITTE FORTSETZEN ...

   238 polynomial2expanded = ...

   239 expanded2polynomial = ...

   241 remark: polynomial2expanded o expanded2polynomial = I,

   242         where 'o' is function chaining, and 'I' is identity WN0210???SK

   244 functions for greatest common divisor and canceling

   245 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   246 ################################################################################

   247 ##   search Isabelle2009-2/src/HOL/Multivariate_Analysis

   248 ##   Amine Chaieb, Robert Himmelmann, John Harrison.

   249 ################################################################################

   250 mv_gcd

   251 factout_

   252 factout_p_

   253 cancel_

   254 cancel_p_

   256 functions for least common multiple and addition of fractions

   257 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   258 mv_lcm

   259 common_nominator_

   260 common_nominator_p_

   261 add_fraction_       (*.add 2 or more fractions.*)

   262 add_fraction_p_     (*.add 2 or more fractions.*)

   264 functions for normalform of rationals

   265 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   266 WN0210???SK interne Funktionen f"ur norm_rational:

   267           schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?

   269 norm_rational_

   270 norm_expanded_rat_

   272 **************************************************************************.*)

   275 (*##*)

   276 structure RationalI : RATIONALI =

   277 struct

   278 (*##*)

   280 infix mem ins union; (*WN100819 updating to Isabelle2009-2*)

   281 fun x mem [] = false

   282   | x mem (y :: ys) = x = y orelse x mem ys;

   286 val is_expanded = is_poly

   287 val is_polynomial = is_poly

   289 fun mk_noteq_0 baseT t =

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

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

   293 fun check_fraction t =

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

   295   in SOME (numerator, denominator) end

   296   handle Bind => NONE

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

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

   300 fun factout_p_ (thy: theory) t =

   301   let val opt = check_fraction t

   302   in

   303     case opt of

   304       NONE => NONE

   305     | SOME (numerator, denominator) =>

   306       let

   307         val vs = t |> vars |> map free2str |> sort string_ord

   308         val baseT = type_of numerator

   309         val expT = HOLogic.realT

   310       in

   311         case (poly_of_term vs numerator, poly_of_term vs denominator) of

   312           (SOME a, SOME b) =>

   313             let

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

   315               val b't = term_of_poly baseT expT vs b'

   316               val ct = term_of_poly baseT expT vs c

   317               val t' =

   318                 HOLogic.mk_binop "Fields.inverse_class.divide"

   319                   (HOLogic.mk_binop "Groups.times_class.times" (term_of_poly baseT expT vs a', ct),

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

   321               val asm = map (mk_noteq_0 baseT) [b't, ct]

   322             in SOME (t', asm) end

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

   324       end

   325   end

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

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

   329 fun cancel_p_ (thy: theory) t =

   330   let val opt = check_fraction t

   331   in

   332     case opt of

   333       NONE => NONE

   334     | SOME (numerator, denominator) =>

   335       let

   336         val vs = t |> vars |> map free2str |> sort string_ord

   337         val baseT = type_of numerator

   338         val expT = HOLogic.realT

   339       in

   340         case (poly_of_term vs numerator, poly_of_term vs denominator) of

   341           (SOME a, SOME b) =>

   342             let

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

   344               val b't = term_of_poly baseT expT vs b'

   345               val ct = term_of_poly baseT expT vs c

   346               val t' =

   347                 HOLogic.mk_binop "Fields.inverse_class.divide"

   348                   (term_of_poly baseT expT vs a', b't)

   349               val asm = map (mk_noteq_0 baseT) [b't, ct]

   350             in SOME (t', asm) end

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

   352       end

   353   end

   355 (* addition of fractions allows (at most) one non-fraction ---postponed after 1st integration*)

   356 fun norm_frac_sum

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

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

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

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

   361   | norm_frac_sum

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

   363       nofrac $

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

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

   366   | norm_frac_sum

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

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

   369       nofrac)

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

   371   | norm_frac_sum _ = NONE

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

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

   375 fun common_nominator_p_ (thy: theory) t =

   376   let val opt = norm_frac_sum t

   377   in

   378     case opt of

   379       NONE => NONE

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

   381       let

   382         val vs = t |> vars |> map free2str |> sort string_ord

   383         val baseT = type_of n1

   384         val expT = HOLogic.realT

   385       in

   386         case (poly_of_term vs d1, poly_of_term vs d2) of

   387           (SOME a, SOME b) =>

   388             let

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

   390               val d1' = term_of_poly baseT expT vs a'

   391               val d2' = term_of_poly baseT expT vs b'

   392               val c' = term_of_poly baseT expT vs c

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

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

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

   396                 --------------------------------------------------------------------------*)

   397               val denom = HOLogic.mk_binop "Groups.times_class.times" (c',

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

   399               (*--------------------------------------------------------------------------*)

   400               val t' =

   401                 HOLogic.mk_binop "Groups.plus_class.plus"

   402                   (HOLogic.mk_binop "Fields.inverse_class.divide"

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

   404                   HOLogic.mk_binop "Fields.inverse_class.divide"

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

   406               val asm = map (mk_noteq_0 baseT) [d1', d2', c']

   407             in SOME (t', asm) end

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

   409       end

   410   end

   412 (* add fractions

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

   414 fun add_fraction_p_ (thy: theory) t =

   415   let val opt = norm_frac_sum t

   416   in

   417     case opt of

   418       NONE => NONE

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

   420       let

   421         val vs = t |> vars |> map free2str |> sort string_ord

   422         val baseT = type_of n1

   423         val expT = HOLogic.realT

   424       in

   425         case (poly_of_term vs d1, poly_of_term vs d2) of

   426           (SOME a, SOME b) =>

   427             let

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

   429               val nomin = term_of_poly baseT expT vs

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

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

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

   433               val asm = [mk_noteq_0 baseT  denom]

   434             in SOME (t', asm) end

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

   436       end

   437   end

   439 fun (x ins xs) = if x mem xs then xs else x :: xs;

   440 fun xs union [] = xs

   441   | [] union ys = ys

   442   | (x :: xs) union ys = xs union (x ins ys);

   444 (*. gcd of integers .*)

   445 (* die gcd Funktion von Isabelle funktioniert nicht richtig !!! *)

   446 fun gcd_int a b = if b=0 then a

   447 		  else gcd_int b (a mod b);

   449 (*. univariate polynomials (uv) .*)

   450 (*. univariate polynomials are represented as a list

   451     of the coefficent in reverse maximum degree order .*)

   452 (*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*)

   453 type uv_poly = int list;

   455 (*. adds two uv polynomials .*)

   456 fun uv_mod_add_poly ([]:uv_poly,p2:uv_poly) = p2:uv_poly

   457   | uv_mod_add_poly (p1,[]) = p1

   458   | uv_mod_add_poly (x::p1,y::p2) = (x+y)::(uv_mod_add_poly(p1,p2));

   460 (*. multiplies a uv polynomial with a skalar s .*)

   461 fun uv_mod_smul_poly ([]:uv_poly,s:int) = []:uv_poly

   462   | uv_mod_smul_poly (x::p,s) = (x*s)::(uv_mod_smul_poly(p,s));

   464 (*. calculates the remainder of a polynomial divided by a skalar s .*)

   465 fun uv_mod_rem_poly ([]:uv_poly,s) = []:uv_poly

   466   | uv_mod_rem_poly (x::p,s) = (x mod s)::(uv_mod_smul_poly(p,s));

   468 (*. calculates the degree of a uv polynomial .*)

   469 fun uv_mod_deg ([]:uv_poly) = 0

   470   | uv_mod_deg p = length(p)-1;

   472 (*. calculates the remainder of x/p and represents it as

   473     value between -p/2 and p/2 .*)

   474 fun uv_mod_mod2(x,p)=

   475     let

   476 	val y=(x mod p);

   477     in

   478 	if (y)>(p div 2) then (y)-p else

   479 	    (

   480 	     if (y)<(~p div 2) then p+(y) else (y)

   481 	     )

   482     end;

   484 (*.calculates the remainder for each element of a integer list divided by p.*)

   485 fun uv_mod_list_modp [] p = []

   486   | uv_mod_list_modp (x::xs) p = (uv_mod_mod2(x,p))::(uv_mod_list_modp xs p);

   488 (*. appends an integer at the end of a integer list .*)

   489 fun uv_mod_null (p1:int list,0) = p1

   490   | uv_mod_null (p1:int list,n1:int) = uv_mod_null(p1,n1-1) @ [0];

   492 (*. uv polynomial division, result is (quotient, remainder) .*)

   493 (*. only for uv_mod_divides .*)

   494 (* FIXME: Division von x^9+x^5+1 durch x-1000 funktioniert nicht,

   495    integer zu klein  *)

   496 fun uv_mod_pdiv (p1:uv_poly) ([]:uv_poly) =

   497     error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

   498   | uv_mod_pdiv p1 [x] =

   499     let

   500 	val xs= Unsynchronized.ref  [];

   501     in

   502 	if x<>0 then

   503 	    (

   504 	     xs:=(uv_mod_rem_poly(p1,x));

   505 	     while length(!xs)>0 andalso hd(!xs)=0 do xs:=tl(!xs)

   506 	     )

   507 	else error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero");

   508 	([]:uv_poly,!xs:uv_poly)

   509     end

   510   | uv_mod_pdiv p1 p2 =

   511     let

   512 	val n= uv_mod_deg(p2);

   513 	val m= Unsynchronized.ref (uv_mod_deg(p1));

   514 	val p1'= Unsynchronized.ref  (rev(p1));

   515 	val p2'=(rev(p2));

   516 	val lc2=hd(p2');

   517 	val q= Unsynchronized.ref  [];

   518 	val c= Unsynchronized.ref  0;

   519 	val output= Unsynchronized.ref  ([],[]);

   520     in

   521 	(

   522 	 if (!m)=0 orelse p2=[0]

   523          then error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero")

   524 	 else

   525 	     (

   526 	      if (!m)<n then

   527 		  (

   528 		   output:=([0],p1)

   529 		   )

   530 	      else

   531 		  (

   532 		   while (!m)>=n do

   533 		       (

   534 			c:=hd(!p1') div hd(p2');

   535 			if !c<>0 then

   536 			    (

   537 			     p1':=uv_mod_add_poly(!p1',uv_mod_null(uv_mod_smul_poly(p2',~(!c)),!m-n));

   538 			     while length(!p1')>0 andalso hd(!p1')=0  do p1':= tl(!p1');

   539 			     m:=uv_mod_deg(!p1')

   540 			     )

   541 			else m:=0

   542 			);

   543     		   output:=(rev(!q),rev(!p1'))

   544 		   )

   545 	      );

   546 	     !output

   547 	 )

   548     end;

   550 (*. divides p1 by p2 in Zp .*)

   551 fun uv_mod_pdivp (p1:uv_poly) (p2:uv_poly) p =

   552     let

   553 	val n=uv_mod_deg(p2);

   554 	val m= Unsynchronized.ref  (uv_mod_deg(uv_mod_list_modp p1 p));

   555 	val p1'= Unsynchronized.ref  (rev(p1));

   556 	val p2'=(rev(uv_mod_list_modp p2 p));

   557 	val lc2=hd(p2');

   558 	val q= Unsynchronized.ref  [];

   559 	val c= Unsynchronized.ref  0;

   560 	val output= Unsynchronized.ref  ([],[]);

   561     in

   562 	(

   563 	 if (!m)=0 orelse p2=[0] then error ("RATIONALS_UV_MOD_PDIVP_EXCEPTION: Division by zero")

   564 	 else

   565 	     (

   566 	      if (!m)<n then

   567 		  (

   568 		   output:=([0],p1)

   569 		   )

   570 	      else

   571 		  (

   572 		   while (!m)>=n do

   573 		       (

   574 			c:=uv_mod_mod2(hd(!p1')*(power lc2 1), p);

   575 			q:=(!c)::(!q);

   576 			p1':=uv_mod_list_modp(tl(uv_mod_add_poly(uv_mod_smul_poly(!p1',lc2),

   577 								  uv_mod_smul_poly(uv_mod_smul_poly(p2',hd(!p1')),~1)))) p;

   578 			m:=(!m)-1

   579 			);

   581 		   while !p1'<>[] andalso hd(!p1')=0 do

   582 		       (

   583 			p1':=tl(!p1')

   584 			);

   586     		   output:=(rev(uv_mod_list_modp (!q) (p)),rev(!p1'))

   587 		   )

   588 	      );

   589 	     !output:uv_poly * uv_poly

   590 	 )

   591     end;

   593 (*. calculates the remainder of p1/p2 .*)

   594 fun uv_mod_prest (p1:uv_poly) ([]:uv_poly) = error("UV_MOD_PREST_EXCEPTION: Division by zero")

   595   | uv_mod_prest [] p2 = []:uv_poly

   596   | uv_mod_prest p1 p2 = (#2(uv_mod_pdiv p1 p2));

   598 (*. calculates the remainder of p1/p2 in Zp .*)

   599 fun uv_mod_prestp (p1:uv_poly) ([]:uv_poly) p= error("UV_MOD_PRESTP_EXCEPTION: Division by zero")

   600   | uv_mod_prestp [] p2 p= []:uv_poly

   601   | uv_mod_prestp p1 p2 p = #2(uv_mod_pdivp p1 p2 p);

   603 (*. calculates the content of a uv polynomial .*)

   604 fun uv_mod_cont ([]:uv_poly) = 0

   605   | uv_mod_cont (x::p)= gcd_int x (uv_mod_cont(p));

   607 (*. divides each coefficient of a uv polynomial by y .*)

   608 fun uv_mod_div_list (p:uv_poly,0) = error("UV_MOD_DIV_LIST_EXCEPTION: Division by zero")

   609   | uv_mod_div_list ([],y)   = []:uv_poly

   610   | uv_mod_div_list (x::p,y) = (x div y)::uv_mod_div_list(p,y);

   612 (*. calculates the primitiv part of a uv polynomial .*)

   613 fun uv_mod_pp ([]:uv_poly) = []:uv_poly

   614   | uv_mod_pp p =

   615     let

   616 	val c= Unsynchronized.ref  0;

   617     in

   618 	(

   619 	 c:=uv_mod_cont(p);

   621 	 if !c=0 then error ("RATIONALS_UV_MOD_PP_EXCEPTION: content is 0")

   622 	 else uv_mod_div_list(p,!c)

   623 	)

   624     end;

   626 (*. gets the leading coefficient of a uv polynomial .*)

   627 fun uv_mod_lc ([]:uv_poly) = 0

   628   | uv_mod_lc p  = hd(rev(p));

   630 (*. calculates the euklidean polynomial remainder sequence in Zp .*)

   631 fun uv_mod_prs_euklid_p(p1:uv_poly,p2:uv_poly,p)=

   632     let

   633 	val f = Unsynchronized.ref  [];

   634 	val f'= Unsynchronized.ref  p2;

   635 	val fi= Unsynchronized.ref  [];

   636     in

   637 	(

   638 	 f:=p2::p1::[];

   639  	 while uv_mod_deg(!f')>0 do

   640 	     (

   641 	      f':=uv_mod_prestp (hd(tl(!f))) (hd(!f)) p;

   642 	      if (!f')<>[] then

   643 		  (

   644 		   fi:=(!f');

   645 		   f:=(!fi)::(!f)

   646 		   )

   647 	      else ()

   648 	      );

   649 	      (!f)

   651 	 )

   652     end;

   654 (*. calculates the gcd of p1 and p2 in Zp .*)

   655 fun uv_mod_gcd_modp ([]:uv_poly) (p2:uv_poly) p = p2:uv_poly

   656   | uv_mod_gcd_modp p1 [] p= p1

   657   | uv_mod_gcd_modp p1 p2 p=

   658     let

   659 	val p1'= Unsynchronized.ref [];

   660 	val p2'= Unsynchronized.ref [];

   661 	val pc= Unsynchronized.ref [];

   662 	val g= Unsynchronized.ref  [];

   663 	val d= Unsynchronized.ref  0;

   664 	val prs= Unsynchronized.ref  [];

   665     in

   666 	(

   667 	 if uv_mod_deg(p1)>=uv_mod_deg(p2) then

   668 	     (

   669 	      p1':=uv_mod_list_modp (uv_mod_pp(p1)) p;

   670 	      p2':=uv_mod_list_modp (uv_mod_pp(p2)) p

   671 	      )

   672 	 else

   673 	     (

   674 	      p1':=uv_mod_list_modp (uv_mod_pp(p2)) p;

   675 	      p2':=uv_mod_list_modp (uv_mod_pp(p1)) p

   676 	      );

   677 	 d:=uv_mod_mod2((gcd_int (uv_mod_cont(p1))) (uv_mod_cont(p2)), p) ;

   678 	 if !d>(p div 2) then d:=(!d)-p else ();

   680 	 prs:=uv_mod_prs_euklid_p(!p1',!p2',p);

   682 	 if hd(!prs)=[] then pc:=hd(tl(!prs))

   683 	 else pc:=hd(!prs);

   685 	 g:=uv_mod_smul_poly(uv_mod_pp(!pc),!d);

   686 	 !g

   687 	 )

   688     end;

   690 (*. calculates the minimum of two real values x and y .*)

   691 fun uv_mod_r_min(x,y):Real.real = if x>y then y else x;

   693 (*. calculates the minimum of two integer values x and y .*)

   694 fun uv_mod_min(x,y) = if x>y then y else x;

   696 (*. adds the squared values of a integer list .*)

   697 fun uv_mod_add_qu [] = 0.0

   698   | uv_mod_add_qu (x::p) =  Real.fromInt(x)*Real.fromInt(x) + uv_mod_add_qu p;

   700 (*. calculates the euklidean norm .*)

   701 fun uv_mod_norm ([]:uv_poly) = 0.0

   702   | uv_mod_norm p = Math.sqrt(uv_mod_add_qu(p));

   704 (*. multipies two values a and b .*)

   705 fun uv_mod_multi a b = a * b;

   707 (*. decides if x is a prim, the list contains all primes which are lower then x .*)

   708 fun uv_mod_prim(x,[])= false

   709   | uv_mod_prim(x,[y])=if ((x mod y) <> 0) then true

   710 		else false

   711   | uv_mod_prim(x,y::ys) = if uv_mod_prim(x,[y])

   712 			then

   713 			    if uv_mod_prim(x,ys) then true

   714 			    else false

   715 		    else false;

   717 (*. gets the first prime, which is greater than p and does not divide g .*)

   718 fun uv_mod_nextprime(g,p)=

   719     let

   720 	val list= Unsynchronized.ref  [2];

   721 	val exit= Unsynchronized.ref  0;

   722 	val i = Unsynchronized.ref 2

   723     in

   724 	while (!i<p) do (* calculates the primes lower then p *)

   725 	    (

   726 	     if uv_mod_prim(!i,!list) then

   727 		 (

   728 		  if (p mod !i <> 0)

   729 		      then

   730 			  (

   731 			   list:= (!i)::(!list);

   732 			   i:= (!i)+1

   733 			   )

   734 		  else i:=(!i)+1

   735 		  )

   736 	     else i:= (!i)+1

   737 		 );

   738 	    i:=(p+1);

   739 	    while (!exit=0) do   (* calculate next prime which does not divide g *)

   740 	    (

   741 	     if uv_mod_prim(!i,!list) then

   742 		 (

   743 		  if (g mod !i <> 0)

   744 		      then

   745 			  (

   746 			   list:= (!i)::(!list);

   747 			   exit:= (!i)

   748 			   )

   749 		  else i:=(!i)+1

   750 		      )

   751 	     else i:= (!i)+1

   752 		 );

   753 	    !exit

   754     end;

   756 (*. decides if p1 is a factor of p2 in Zp .*)

   757 fun uv_mod_dividesp ([]:uv_poly) (p2:uv_poly) p= error("UV_MOD_DIVIDESP: Division by zero")

   758   | uv_mod_dividesp p1 p2 p= if uv_mod_prestp p2 p1 p = [] then true else false;

   760 (*. decides if p1 is a factor of p2 .*)

   761 fun uv_mod_divides ([]:uv_poly) (p2:uv_poly) = error("UV_MOD_DIVIDES: Division by zero")

   762   | uv_mod_divides p1 p2 = if uv_mod_prest p2 p1  = [] then true else false;

   764 (*. chinese remainder algorithm .*)

   765 fun uv_mod_cra2(r1,r2,m1,m2)=

   766     let

   767 	val c= Unsynchronized.ref  0;

   768 	val r1'= Unsynchronized.ref  0;

   769 	val d= Unsynchronized.ref  0;

   770 	val a= Unsynchronized.ref  0;

   771     in

   772 	(

   773 	 while (uv_mod_mod2((!c)*m1,m2))<>1 do

   774 	     (

   775 	      c:=(!c)+1

   776 	      );

   777 	 r1':= uv_mod_mod2(r1,m1);

   778 	 d:=uv_mod_mod2(((r2-(!r1'))*(!c)),m2);

   779 	 !r1'+(!d)*m1

   780 	 )

   781     end;

   783 (*. applies the chinese remainder algorithmen to the coefficients of x1 and x2 .*)

   784 fun uv_mod_cra_2 ([],[],m1,m2) = []

   785   | uv_mod_cra_2 ([],x2,m1,m2) = error("UV_MOD_CRA_2_EXCEPTION: invalid call x1")

   786   | uv_mod_cra_2 (x1,[],m1,m2) = error("UV_MOD_CRA_2_EXCEPTION: invalid call x2")

   787   | 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));

   789 (*. calculates the gcd of two uv polynomials p1' and p2' with the modular algorithm .*)

   790 fun uv_mod_gcd (p1':uv_poly) (p2':uv_poly) =

   791     let

   792 	val p1= Unsynchronized.ref  (uv_mod_pp(p1'));

   793 	val p2= Unsynchronized.ref  (uv_mod_pp(p2'));

   794 	val c=gcd_int (uv_mod_cont(p1')) (uv_mod_cont(p2'));

   795 	val temp= Unsynchronized.ref  [];

   796 	val cp= Unsynchronized.ref  [];

   797 	val qp= Unsynchronized.ref  [];

   798 	val q= Unsynchronized.ref [];

   799 	val pn= Unsynchronized.ref  0;

   800 	val d= Unsynchronized.ref  0;

   801 	val g1= Unsynchronized.ref  0;

   802 	val p= Unsynchronized.ref  0;

   803 	val m= Unsynchronized.ref  0;

   804 	val exit= Unsynchronized.ref  0;

   805 	val i= Unsynchronized.ref  1;

   806     in

   807 	if length(!p1)>length(!p2) then ()

   808 	else

   809 	    (

   810 	     temp:= !p1;

   811 	     p1:= !p2;

   812 	     p2:= !temp

   813 	     );

   816 	d:=gcd_int (uv_mod_lc(!p1)) (uv_mod_lc(!p2));

   817 	g1:=uv_mod_lc(!p1)*uv_mod_lc(!p2);

   818 	p:=4;

   820 	m:=Real.ceil(2.0 * Real.fromInt(!d) *

   821 	  Real.fromInt(power 2 (uv_mod_min(uv_mod_deg(!p2),uv_mod_deg(!p1)))) *

   822 	  Real.fromInt(!d) *

   823 	  uv_mod_r_min(uv_mod_norm(!p1) / Real.fromInt(abs(uv_mod_lc(!p1))),

   824 	  uv_mod_norm(!p2) / Real.fromInt(abs(uv_mod_lc(!p2)))));

   826 	while (!exit=0) do

   827 	    (

   828 	     p:=uv_mod_nextprime(!d,!p);

   829 	     cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)) ;

   830 	     if abs(uv_mod_lc(!cp))<>1 then  (* leading coefficient = 1 ? *)

   831 		 (

   832 		  i:=1;

   833 		  while (!i)<(!p) andalso (abs(uv_mod_mod2((uv_mod_lc(!cp)*(!i)),(!p)))<>1) do

   834 		      (

   835 		       i:=(!i)+1

   836 		       );

   837 		      cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)

   838 		  )

   839 	     else ();

   841 	     qp:= ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp));

   843 	     if uv_mod_deg(!qp)=0 then (q:=[1]; exit:=1) else ();

   845 	     pn:=(!p);

   846 	     q:=(!qp);

   848 	     while !pn<= !m andalso !m>(!p) andalso !exit=0 do

   849 		 (

   850 		  p:=uv_mod_nextprime(!d,!p);

   851  		  cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p));

   852  		  if uv_mod_lc(!cp)<>1 then  (* leading coefficient = 1 ? *)

   853  		      (

   854  		       i:=1;

   855  		       while (!i)<(!p) andalso ((uv_mod_mod2((uv_mod_lc(!q)*(!i)),(!p)))<>1) do

   856  			   (

   857  			    i:=(!i)+1

   858 		           );

   859 		       cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p)

   860  		      )

   861  		  else ();

   863 		  qp:=uv_mod_list_modp ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp)  ) (!p);

   864  		  if uv_mod_deg(!qp)>uv_mod_deg(!q) then

   865  		      (

   866  		       (*print("degree to high!!!\n")*)

   867  		       )

   868  		  else

   869  		      (

   870  		       if uv_mod_deg(!qp)=uv_mod_deg(!q) then

   871  			   (

   872  			    q:=uv_mod_cra_2(!q,!qp,!pn,!p);

   873 			    pn:=(!pn) * !p;

   874 			    q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn)); (* found already gcd ? *)

   875 			    if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then (exit:=1) else ()

   876 		 	    )

   877 		       else

   878 			   (

   879 			    if  uv_mod_deg(!qp)<uv_mod_deg(!q) then

   880 				(

   881 				 pn:= !p;

   882 				 q:= !qp

   883 				 )

   884 			    else ()

   885 			    )

   886 		       )

   887 		  );

   888  	     q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn));

   889 	     if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then exit:=1 else ()

   890 	     );

   891 	    uv_mod_smul_poly(!q,c):uv_poly

   892     end;

   894 (*. multivariate polynomials .*)

   895 (*. multivariate polynomials are represented as a list of the pairs,

   896  first is the coefficent and the second is a list of the exponents .*)

   897 (*. 5 * x^5 * y^3 + 4 * x^3 * z^2 + 2 * x^2 * y * z^3 - z - 19

   898  => [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*)

   900 (*. global variables .*)

   901 (*. order indicators .*)

   902 val LEX_=0; (* lexicographical term order *)

   903 val GGO_=1; (* greatest degree order *)

   905 (*. datatypes for internal representation.*)

   906 type mv_monom = (int *      (*.coefficient or the monom.*)

   907 		 int list); (*.list of exponents)      .*)

   908 fun mv_monom2str (i, is) = "("^ int2str i^"," ^ ints2str' is ^ ")";

   910 type mv_poly = mv_monom list;

   911 fun mv_poly2str p = (strs2str' o (map mv_monom2str)) p;

   913 (*. help function for monom_greater and geq .*)

   914 fun mv_mg_hlp([]) = EQUAL

   915   | mv_mg_hlp(x::list)=if x<0 then LESS

   916 		    else if x>0 then GREATER

   917 			 else mv_mg_hlp(list);

   919 (*. adds a list of values .*)

   920 fun mv_addlist([]) = 0

   921   | mv_addlist(p1) = hd(p1)+mv_addlist(tl(p1));

   923 (*. tests if the monomial M1 is greater as the monomial M2 and returns a boolean value .*)

   924 (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)

   925 fun mv_monom_greater((M1x,M1l):mv_monom,(M2x,M2l):mv_monom,order)=

   926     if order=LEX_ then

   927 	(

   928 	 if length(M1l)<>length(M2l) then error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")

   929 	 else if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true

   930 	     )

   931     else

   932 	if order=GGO_ then

   933 	    (

   934 	     if length(M1l)<>length(M2l) then error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error")

   935 	     else

   936 		 if mv_addlist(M1l)=mv_addlist(M2l)  then if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true

   937 		 else if mv_addlist(M1l)>mv_addlist(M2l) then true else false

   938 	     )

   939 	else error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");

   941 (*. tests if the monomial X is greater as the monomial Y and returns a order value (GREATER,EQUAL,LESS) .*)

   942 (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*)

   943 fun mv_geq order ((x1,x):mv_monom,(x2,y):mv_monom) =

   944 let

   945     val temp= Unsynchronized.ref  EQUAL;

   946 in

   947     if order=LEX_ then

   948 	(

   949 	 if length(x)<>length(y) then

   950 	     error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   951 	 else

   952 	     (

   953 	      temp:=mv_mg_hlp((map op- (x~~y)));

   954 	      if !temp=EQUAL then

   955 		  ( if x1=x2 then EQUAL

   956 		    else if x1>x2 then GREATER

   957 			 else LESS

   958 			     )

   959 	      else (!temp)

   960 	      )

   961 	     )

   962     else

   963 	if order=GGO_ then

   964 	    (

   965 	     if length(x)<>length(y) then

   966 	      error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   967 	     else

   968 		 if mv_addlist(x)=mv_addlist(y) then

   969 		     (mv_mg_hlp((map op- (x~~y))))

   970 		 else if mv_addlist(x)>mv_addlist(y) then GREATER else LESS

   971 		     )

   972 	else error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")

   973 end;

   975 (*. cuts the first variable from a polynomial .*)

   976 fun mv_cut([]:mv_poly)=[]:mv_poly

   977   | mv_cut((x,[])::list) = error ("RATIONALS_MV_CUT_EXCEPTION: Invalid list ")

   978   | mv_cut((x,y::ys)::list)=(x,ys)::mv_cut(list);

   980 (*. leading power product .*)

   981 fun mv_lpp([]:mv_poly,order)  = []

   982   | mv_lpp([(x,y)],order) = y

   983   | mv_lpp(p1,order)  = #2(hd(rev(sort (mv_geq order) p1)));

   985 (*. leading monomial .*)

   986 fun mv_lm([]:mv_poly,order)  = (0,[]):mv_monom

   987   | mv_lm([x],order) = x

   988   | mv_lm(p1,order)  = hd(rev(sort (mv_geq order) p1));

   990 (*. leading coefficient in term order .*)

   991 fun mv_lc2([]:mv_poly,order)  = 0

   992   | mv_lc2([(x,y)],order) = x

   993   | mv_lc2(p1,order)  = #1(hd(rev(sort (mv_geq order) p1)));

   996 (*. reverse the coefficients in mv polynomial .*)

   997 fun mv_rev_to([]:mv_poly) = []:mv_poly

   998   | mv_rev_to((c,e)::xs) = (c,rev(e))::mv_rev_to(xs);

  1000 (*. leading coefficient in reverse term order .*)

  1001 fun mv_lc([]:mv_poly,order)  = []:mv_poly

  1002   | mv_lc([(x,y)],order) = mv_rev_to(mv_cut(mv_rev_to([(x,y)])))

  1003   | mv_lc(p1,order)  =

  1004     let

  1005 	val p1o= Unsynchronized.ref  (rev(sort (mv_geq order) (mv_rev_to(p1))));

  1006 	val lp=hd(#2(hd(!p1o)));

  1007 	val lc= Unsynchronized.ref  [];

  1008     in

  1009 	(

  1010 	 while (length(!p1o)>0 andalso hd(#2(hd(!p1o)))=lp) do

  1011 	     (

  1012 	      lc:=hd(mv_cut([hd(!p1o)]))::(!lc);

  1013 	      p1o:=tl(!p1o)

  1014 	      );

  1015 	 if !lc=[] then error ("RATIONALS_MV_LC_EXCEPTION: lc is empty") else ();

  1016 	 mv_rev_to(!lc)

  1017 	 )

  1018     end;

  1020 (*. compares two powerproducts .*)

  1021 fun mv_monom_equal((_,xlist):mv_monom,(_,ylist):mv_monom) = (foldr and_) (((map op=) (xlist~~ylist)),true);

  1023 (*. help function for mv_add .*)

  1024 fun mv_madd([]:mv_poly,[]:mv_poly,order) = []:mv_poly

  1025   | mv_madd([(0,_)],p2,order) = p2

  1026   | mv_madd(p1,[(0,_)],order) = p1

  1027   | mv_madd([],p2,order) = p2

  1028   | mv_madd(p1,[],order) = p1

  1029   | mv_madd(p1,p2,order) =

  1030     (

  1031      if mv_monom_greater(hd(p1),hd(p2),order)

  1032 	 then hd(p1)::mv_madd(tl(p1),p2,order)

  1033      else if mv_monom_equal(hd(p1),hd(p2))

  1034 	      then if mv_lc2(p1,order)+mv_lc2(p2,order)<>0

  1035 		       then (mv_lc2(p1,order)+mv_lc2(p2,order),mv_lpp(p1,order))::mv_madd(tl(p1),tl(p2),order)

  1036 		   else mv_madd(tl(p1),tl(p2),order)

  1037 	  else hd(p2)::mv_madd(p1,tl(p2),order)

  1038 	      )

  1040 (*. adds two multivariate polynomials .*)

  1041 fun mv_add([]:mv_poly,p2:mv_poly,order) = p2

  1042   | mv_add(p1,[],order) = p1

  1043   | mv_add(p1,p2,order) = mv_madd(rev(sort (mv_geq order) p1),rev(sort (mv_geq order) p2), order);

  1045 (*. monom multiplication .*)

  1046 fun mv_mmul((x1,y1):mv_monom,(x2,y2):mv_monom)=(x1*x2,(map op+) (y1~~y2)):mv_monom;

  1048 (*. deletes all monomials with coefficient 0 .*)

  1049 fun mv_shorten([]:mv_poly,order) = []:mv_poly

  1050   | mv_shorten(x::xs,order)=mv_madd([x],mv_shorten(xs,order),order);

  1052 (*. zeros a list .*)

  1053 fun mv_null2([])=[]

  1054   | mv_null2(x::l)=0::mv_null2(l);

  1056 (*. multiplies two multivariate polynomials .*)

  1057 fun mv_mul([]:mv_poly,[]:mv_poly,_) = []:mv_poly

  1058   | mv_mul([],y::p2,_) = [(0,mv_null2(#2(y)))]

  1059   | mv_mul(x::p1,[],_) = [(0,mv_null2(#2(x)))]

  1060   | mv_mul(x::p1,y::p2,order) = mv_shorten(rev(sort (mv_geq order) (mv_mmul(x,y) :: (mv_mul(p1,y::p2,order) @

  1061 									    mv_mul([x],p2,order)))),order);

  1063 (*. gets the maximum value of a list .*)

  1064 fun mv_getmax([])=0

  1065   | mv_getmax(x::p1)= let

  1066 		       val m=mv_getmax(p1);

  1067 		   in

  1068 		       if m>x then m

  1069 		       else x

  1070 		   end;

  1071 (*. calculates the maximum degree of an multivariate polynomial .*)

  1072 fun mv_grad([]:mv_poly) = 0

  1073   | mv_grad(p1:mv_poly)= mv_getmax((map mv_addlist) ((map #2) p1));

  1075 (*. converts the sign of a value .*)

  1076 fun mv_minus(x)=(~1) * x;

  1078 (*. converts the sign of all coefficients of a polynomial .*)

  1079 fun mv_minus2([]:mv_poly)=[]:mv_poly

  1080   | mv_minus2(p1)=(mv_minus(#1(hd(p1))),#2(hd(p1)))::(mv_minus2(tl(p1)));

  1082 (*. searches for a negativ value in a list .*)

  1083 fun mv_is_negativ([])=false

  1084   | mv_is_negativ(x::xs)=if x<0 then true else mv_is_negativ(xs);

  1086 (*. division of monomials .*)

  1087 fun mv_mdiv((0,[]):mv_monom,_:mv_monom)=(0,[]):mv_monom

  1088   | mv_mdiv(_,(0,[]))= error ("RATIONALS_MV_MDIV_EXCEPTION Division by 0 ")

  1089   | mv_mdiv(p1:mv_monom,p2:mv_monom)=

  1090     let

  1091 	val c= Unsynchronized.ref  (#1(p2));

  1092 	val pp= Unsynchronized.ref  [];

  1093     in

  1094 	(

  1095 	 if !c=0 then error("MV_MDIV_EXCEPTION Dividing by zero")

  1096 	 else c:=(#1(p1) div #1(p2));

  1097 	     if #1(p2)<>0 then

  1098 		 (

  1099 		  pp:=(#2(mv_mmul((1,#2(p1)),(1,(map mv_minus) (#2(p2))))));

  1100 		  if mv_is_negativ(!pp) then (0,!pp)

  1101 		  else (!c,!pp)

  1102 		      )

  1103 	     else error("MV_MDIV_EXCEPTION Dividing by empty Polynom")

  1104 		 )

  1105     end;

  1107 (*. prints a polynom for (internal use only) .*)

  1108 fun mv_print_poly([]:mv_poly)=tracing("[]\n")

  1109   | mv_print_poly((x,y)::[])= tracing("("^Int.toString(x)^","^ints2str(y)^")\n")

  1110   | mv_print_poly((x,y)::p1) = (tracing("("^Int.toString(x)^","^ints2str(y)^"),");mv_print_poly(p1));

  1113 (*. division of two multivariate polynomials .*)

  1114 fun mv_division([]:mv_poly,g:mv_poly,order)=([]:mv_poly,[]:mv_poly)

  1115   | mv_division(f,[],order)= error ("RATIONALS_MV_DIVISION_EXCEPTION Division by zero")

  1116   | mv_division(f,g,order)=

  1117     let

  1118 	val r= Unsynchronized.ref  [];

  1119 	val q= Unsynchronized.ref  [];

  1120 	val g'= Unsynchronized.ref  ([] : mv_monom list);

  1121 	val k= Unsynchronized.ref  0;

  1122 	val m= Unsynchronized.ref  (0,[0]);

  1123 	val exit= Unsynchronized.ref  0;

  1124     in

  1125 	r := rev(sort (mv_geq order) (mv_shorten(f,order)));

  1126 	g':= rev(sort (mv_geq order) (mv_shorten(g,order)));

  1127 	if #1(hd(!g'))=0 then error("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero") else ();

  1128 	if  (mv_monom_greater (hd(!g'),hd(!r),order)) then ([(0,mv_null2(#2(hd(f))))],(!r))

  1129 	else

  1130 	    (

  1131 	     exit:=0;

  1132 	     while (if (!exit)=0 then not(mv_monom_greater (hd(!g'),hd(!r),order)) else false) do

  1133 		 (

  1134 		  if (#1(mv_lm(!g',order)))<>0 then m:=mv_mdiv(mv_lm(!r,order),mv_lm(!g',order))

  1135 		  else error ("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero");

  1136 		  if #1(!m)<>0 then

  1137 		      (

  1138 		       q:=(!m)::(!q);

  1139 		       r:=mv_add((!r),mv_minus2(mv_mul(!g',[!m],order)),order)

  1140 		       )

  1141 		  else exit:=1;

  1142 		  if (if length(!r)<>0 then length(!g')<>0 else false) then ()

  1143 		  else (exit:=1)

  1144 		  );

  1145 		 (rev(!q),!r)

  1146 		 )

  1147     end;

  1149 (*. multiplies a polynomial with an integer .*)

  1150 fun mv_skalar_mul([]:mv_poly,c) = []:mv_poly

  1151   | mv_skalar_mul((x,y)::p1,c) = ((x * c),y)::mv_skalar_mul(p1,c);

  1153 (*. inserts the a first variable into an polynomial with exponent v .*)

  1154 fun mv_correct([]:mv_poly,v:int)=[]:mv_poly

  1155   | mv_correct((x,y)::list,v:int)=(x,v::y)::mv_correct(list,v);

  1157 (*. multivariate case .*)

  1159 (*. decides if x is a factor of y .*)

  1160 fun mv_divides([]:mv_poly,[]:mv_poly)=  error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")

  1161   | mv_divides(x,[]) =  error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero")

  1162   | mv_divides(x:mv_poly,y:mv_poly) = #2(mv_division(y,x,LEX_))=[];

  1164 (*. gets the maximum of a and b .*)

  1165 fun mv_max(a,b) = if a>b then a else b;

  1167 (*. gets the maximum exponent of a mv polynomial in the lexicographic term order .*)

  1168 fun mv_deg([]:mv_poly) = 0

  1169   | mv_deg(p1)=

  1170     let

  1171 	val p1'=mv_shorten(p1,LEX_);

  1172     in

  1173 	if length(p1')=0 then 0

  1174 	else mv_max(hd(#2(hd(p1'))),mv_deg(tl(p1')))

  1175     end;

  1177 (*. gets the maximum exponent of a mv polynomial in the reverse lexicographic term order .*)

  1178 fun mv_deg2([]:mv_poly) = 0

  1179   | mv_deg2(p1)=

  1180     let

  1181 	val p1'=mv_shorten(p1,LEX_);

  1182     in

  1183 	if length(p1')=0 then 0

  1184 	else mv_max(hd(rev(#2(hd(p1')))),mv_deg2(tl(p1')))

  1185     end;

  1187 (*. evaluates the mv polynomial at the value v of the main variable .*)

  1188 fun mv_subs([]:mv_poly,v) = []:mv_poly

  1189   | mv_subs((c,e)::p1:mv_poly,v) = mv_skalar_mul(mv_cut([(c,e)]),power v (hd(e))) @ mv_subs(p1,v);

  1191 (*. calculates the content of a uv-polynomial in mv-representation .*)

  1192 fun uv_content2([]:mv_poly) = 0

  1193   | uv_content2((c,e)::p1) = (gcd_int c (uv_content2(p1)));

  1195 (*. converts a uv-polynomial from mv-representation to  uv-representation .*)

  1196 fun uv_to_list ([]:mv_poly)=[]:uv_poly

  1197   | uv_to_list ((c1,e1)::others) =

  1198     let

  1199 	val count= Unsynchronized.ref  0;

  1200 	val max=mv_grad((c1,e1)::others);

  1201 	val help= Unsynchronized.ref  ((c1,e1)::others);

  1202 	val list= Unsynchronized.ref  [];

  1203     in

  1204 	if length(e1)>1 then error ("RATIONALS_TO_LIST_EXCEPTION: not univariate")

  1205 	else if length(e1)=0 then [c1]

  1206 	     else

  1207 		 (

  1208 		  count:=0;

  1209 		  while (!count)<=max do

  1210 		      (

  1211 		       if length(!help)>0 andalso hd(#2(hd(!help)))=max-(!count) then

  1212 			   (

  1213 			    list:=(#1(hd(!help)))::(!list);

  1214 			    help:=tl(!help)

  1215 			    )

  1216 		       else

  1217 			   (

  1218 			    list:= 0::(!list)

  1219 			    );

  1220 		       count := (!count) + 1

  1221 		       );

  1222 		      (!list)

  1223 		      )

  1224     end;

  1226 (*. converts a uv-polynomial from uv-representation to mv-representation .*)

  1227 fun uv_to_poly ([]:uv_poly) = []:mv_poly

  1228   | uv_to_poly p1 =

  1229     let

  1230 	val count= Unsynchronized.ref  0;

  1231 	val help= Unsynchronized.ref  p1;

  1232 	val list= Unsynchronized.ref  [];

  1233     in

  1234 	while length(!help)>0 do

  1235 	    (

  1236 	     if hd(!help)=0 then ()

  1237 	     else list:=(hd(!help),[!count])::(!list);

  1238 	     count:=(!count)+1;

  1239 	     help:=tl(!help)

  1240 	     );

  1241 	    (!list)

  1242     end;

  1244 (*. univariate gcd calculation from polynomials in multivariate representation .*)

  1245 fun uv_gcd ([]:mv_poly) p2 = p2

  1246   | uv_gcd p1 ([]:mv_poly) = p1

  1247   | uv_gcd p1 [(c,[e])] =

  1248     let

  1249 	val list= Unsynchronized.ref  (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));

  1250 	val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));

  1251     in

  1252 	[(gcd_int (uv_content2(p1)) c,[min])]

  1253     end

  1254   | uv_gcd [(c,[e])] p2 =

  1255     let

  1256 	val list= Unsynchronized.ref  (rev(sort (mv_geq LEX_) (mv_shorten(p2,LEX_))));

  1257 	val min=uv_mod_min(e,(hd(#2(hd(rev(!list))))));

  1258     in

  1259 	[(gcd_int (uv_content2(p2)) c,[min])]

  1260     end

  1261   | uv_gcd p11 p22 = uv_to_poly(uv_mod_gcd (uv_to_list(mv_shorten(p11,LEX_))) (uv_to_list(mv_shorten(p22,LEX_))));

  1263 (*. help function for the newton interpolation .*)

  1264 fun mv_newton_help ([]:mv_poly list,k:int) = []:mv_poly list

  1265   | mv_newton_help (pl:mv_poly list,k) =

  1266     let

  1267 	val x= Unsynchronized.ref  (rev(pl));

  1268 	val t= Unsynchronized.ref  [];

  1269 	val y= Unsynchronized.ref  [];

  1270 	val n= Unsynchronized.ref  1;

  1271 	val n1= Unsynchronized.ref [];

  1272     in

  1273 	(

  1274 	 while length(!x)>1 do

  1275 	     (

  1276 	      if length(hd(!x))>0 then n1:=mv_null2(#2(hd(hd(!x))))

  1277 	      else if length(hd(tl(!x)))>0 then n1:=mv_null2(#2(hd(hd(tl(!x)))))

  1278 		   else n1:=[];

  1279 	      t:= #1(mv_division(mv_add(hd(!x),mv_skalar_mul(hd(tl(!x)),~1),LEX_),[(k,!n1)],LEX_));

  1280 	      y:=(!t)::(!y);

  1281 	      x:=tl(!x)

  1282 	      );

  1283 	 (!y)

  1284 	 )

  1285     end;

  1287 (*. help function for the newton interpolation .*)

  1288 fun mv_newton_add ([]:mv_poly list) t = []:mv_poly

  1289   | mv_newton_add [x:mv_poly] t = x

  1290   | mv_newton_add (pl:mv_poly list) t =

  1291     let

  1292 	val expos= Unsynchronized.ref  [];

  1293 	val pll= Unsynchronized.ref  pl;

  1294     in

  1295 	(

  1297 	 while length(!pll)>0 andalso hd(!pll)=[]  do

  1298 	     (

  1299 	      pll:=tl(!pll)

  1300 	      );

  1301 	 if length(!pll)>0 then expos:= #2(hd(hd(!pll))) else expos:=[];

  1302 	 mv_add(hd(pl),

  1303 		mv_mul(

  1304 		       mv_add(mv_correct(mv_cut([(1,mv_null2(!expos))]),1),[(~t,mv_null2(!expos))],LEX_),

  1305 		       mv_newton_add (tl(pl)) (t+1),

  1306 		       LEX_

  1307 		       ),

  1308 		LEX_)

  1309 	 )

  1310     end;

  1312 (*. calculates the newton interpolation with polynomial coefficients .*)

  1313 (*. step-depth is 1 and if the result is not an integerpolynomial .*)

  1314 (*. this function returns [] .*)

  1315 fun mv_newton ([]:(mv_poly) list) = []:mv_poly

  1316   | mv_newton ([mp]:(mv_poly) list) = mp:mv_poly

  1317   | mv_newton pl =

  1318     let

  1319 	val c= Unsynchronized.ref  pl;

  1320 	val c1= Unsynchronized.ref  [];

  1321 	val n=length(pl);

  1322 	val k= Unsynchronized.ref  1;

  1323 	val i= Unsynchronized.ref  n;

  1324 	val ppl= Unsynchronized.ref  [];

  1325     in

  1326 	c1:=hd(pl)::[];

  1327 	c:=mv_newton_help(!c,!k);

  1328 	c1:=(hd(!c))::(!c1);

  1329 	while(length(!c)>1 andalso !k<n) do

  1330 	    (

  1331 	     k:=(!k)+1;

  1332 	     while  length(!c)>0 andalso hd(!c)=[] do c:=tl(!c);

  1333 	     if !c=[] then () else c:=mv_newton_help(!c,!k);

  1334 	     ppl:= !c;

  1335 	     if !c=[] then () else  c1:=(hd(!c))::(!c1)

  1336 	     );

  1337 	while  hd(!c1)=[] do c1:=tl(!c1);

  1338 	c1:=rev(!c1);

  1339 	ppl:= !c1;

  1340 	mv_newton_add (!c1) 1

  1341     end;

  1343 (*. sets the exponents of the first variable to zero .*)

  1344 fun mv_null3([]:mv_poly)    = []:mv_poly

  1345   | mv_null3((x,y)::xs) = (x,0::tl(y))::mv_null3(xs);

  1347 (*. calculates the minimum exponents of a multivariate polynomial .*)

  1348 fun mv_min_pp([]:mv_poly)=[]

  1349   | mv_min_pp((c,e)::xs)=

  1350     let

  1351 	val y= Unsynchronized.ref  xs;

  1352 	val x= Unsynchronized.ref  [];

  1353     in

  1354 	(

  1355 	 x:=e;

  1356 	 while length(!y)>0 do

  1357 	     (

  1358 	      x:=(map uv_mod_min) ((!x) ~~ (#2(hd(!y))));

  1359 	      y:=tl(!y)

  1360 	      );

  1361 	 !x

  1362 	 )

  1363     end;

  1365 (*. checks if all elements of the list have value zero .*)

  1366 fun list_is_null [] = true

  1367   | list_is_null (x::xs) = (x=0 andalso list_is_null(xs));

  1369 (* check if main variable is zero*)

  1370 fun main_zero (ms : mv_poly) = (list_is_null o (map (hd o #2))) ms;

  1372 (*. calculates the content of an polynomial .*)

  1373 fun mv_content([]:mv_poly) = []:mv_poly

  1374   | mv_content(p1) =

  1375     let

  1376 	val list= Unsynchronized.ref  (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_))));

  1377 	val test= Unsynchronized.ref  (hd(#2(hd(!list))));

  1378 	val result= Unsynchronized.ref  [];

  1379 	val min=(hd(#2(hd(rev(!list)))));

  1380     in

  1381 	(

  1382 	 if length(!list)>1 then

  1383 	     (

  1384 	      while (if length(!list)>0 then (hd(#2(hd(!list)))=(!test)) else false) do

  1385 		  (

  1386 		   result:=(#1(hd(!list)),tl(#2(hd(!list))))::(!result);

  1388 		   if length(!list)<1 then list:=[]

  1389 		   else list:=tl(!list)

  1391 		       );

  1392 		  if length(!list)>0 then

  1393 		   (

  1394 		    list:=mv_gcd (!result) (mv_cut(mv_content(!list)))

  1395 		    )

  1396 		  else list:=(!result);

  1397 		  list:=mv_correct(!list,0);

  1398 		  (!list)

  1399 		  )

  1400 	 else

  1401 	     (

  1402 	      mv_null3(!list)

  1403 	      )

  1404 	     )

  1405     end

  1407 (*. calculates the primitiv part of a polynomial .*)

  1408 and mv_pp([]:mv_poly) = []:mv_poly

  1409   | mv_pp(p1) = let

  1410 		    val cont= Unsynchronized.ref  [];

  1411 		    val pp= Unsynchronized.ref [];

  1412 		in

  1413 		    cont:=mv_content(p1);

  1414 		    pp:=(#1(mv_division(p1,!cont,LEX_)));

  1415 		    if !pp=[]

  1416 			then error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")

  1417 		    else (!pp)

  1418 		end

  1420 (*. calculates the gcd of two multivariate polynomials with a modular approach .*)

  1421 and mv_gcd ([]:mv_poly) ([]:mv_poly) :mv_poly= []:mv_poly

  1422   | mv_gcd ([]:mv_poly) (p2) :mv_poly= p2:mv_poly

  1423   | mv_gcd (p1:mv_poly) ([]) :mv_poly= p1:mv_poly

  1424   | mv_gcd ([(x,xs)]:mv_poly) ([(y,ys)]):mv_poly =

  1425      let

  1426       val xpoly:mv_poly = [(x,xs)];

  1427       val ypoly:mv_poly = [(y,ys)];

  1428      in

  1429 	(

  1430 	 if xs=ys then [((gcd_int x y),xs)]

  1431 	 else [((gcd_int x y),(map uv_mod_min)(xs~~ys))]:mv_poly

  1432         )

  1433     end

  1434   | mv_gcd (p1:mv_poly) ([(y,ys)]) :mv_poly=

  1435 	(

  1436 	 [(gcd_int (uv_content2(p1)) (y),(map  uv_mod_min)(mv_min_pp(p1)~~ys))]:mv_poly

  1437 	)

  1438   | mv_gcd ([(y,ys)]:mv_poly) (p2):mv_poly =

  1439 	(

  1440          [(gcd_int (uv_content2(p2)) (y),(map  uv_mod_min)(mv_min_pp(p2)~~ys))]:mv_poly

  1441         )

  1442   | mv_gcd (p1':mv_poly) (p2':mv_poly):mv_poly=

  1443     let

  1444 	val vc=length(#2(hd(p1')));

  1445 	val cont =

  1446 		  (

  1447                    if main_zero(mv_content(p1')) andalso

  1448                      (main_zero(mv_content(p2'))) then

  1449                      mv_correct((mv_gcd (mv_cut(mv_content(p1'))) (mv_cut(mv_content(p2')))),0)

  1450                    else

  1451                      mv_gcd (mv_content(p1')) (mv_content(p2'))

  1452                   );

  1453 	val p1= #1(mv_division(p1',mv_content(p1'),LEX_));

  1454 	val p2= #1(mv_division(p2',mv_content(p2'),LEX_));

  1455 	val gcd= Unsynchronized.ref  [];

  1456 	val candidate= Unsynchronized.ref  [];

  1457 	val interpolation_list= Unsynchronized.ref  [];

  1458 	val delta= Unsynchronized.ref  [];

  1459         val p1r = Unsynchronized.ref [];

  1460         val p2r = Unsynchronized.ref [];

  1461         val p1r' = Unsynchronized.ref [];

  1462         val p2r' = Unsynchronized.ref [];

  1463 	val factor= Unsynchronized.ref  [];

  1464 	val r= Unsynchronized.ref  0;

  1465 	val gcd_r= Unsynchronized.ref  [];

  1466 	val d= Unsynchronized.ref  0;

  1467 	val exit= Unsynchronized.ref  0;

  1468 	val current_degree= Unsynchronized.ref  99999; (*. FIXME: unlimited ! .*)

  1469     in

  1470 	(

  1471 	 if vc<2 then (* areUnivariate(p1',p2') *)

  1472 	     (

  1473 	      gcd:=uv_gcd (mv_shorten(p1',LEX_)) (mv_shorten(p2',LEX_))

  1474 	      )

  1475 	 else

  1476 	     (

  1477 	      while !exit=0 do

  1478 		  (

  1479 		   r:=(!r)+1;

  1480                    p1r := mv_lc(p1,LEX_);

  1481 		   p2r := mv_lc(p2,LEX_);

  1482                    if main_zero(!p1r) andalso

  1483                       main_zero(!p2r)

  1484                    then

  1485                        (

  1486                         delta := mv_correct((mv_gcd (mv_cut (!p1r)) (mv_cut (!p2r))),0)

  1487                        )

  1488                    else

  1489                        (

  1490 		        delta := mv_gcd (!p1r) (!p2r)

  1491                        );

  1492 		   (*if mv_shorten(mv_subs(!p1r,!r),LEX_)=[] andalso

  1493 		      mv_shorten(mv_subs(!p2r,!r),LEX_)=[] *)

  1494 		   if mv_lc2(mv_shorten(mv_subs(!p1r,!r),LEX_),LEX_)=0 andalso

  1495 		      mv_lc2(mv_shorten(mv_subs(!p2r,!r),LEX_),LEX_)=0

  1496                    then

  1497                        (

  1498 		       )

  1499 		   else

  1500 		       (

  1501 			gcd_r:=mv_shorten(mv_gcd (mv_shorten(mv_subs(p1,!r),LEX_))

  1502 					         (mv_shorten(mv_subs(p2,!r),LEX_)) ,LEX_);

  1503 			gcd_r:= #1(mv_division(mv_mul(mv_correct(mv_subs(!delta,!r),0),!gcd_r,LEX_),

  1504 					       mv_correct(mv_lc(!gcd_r,LEX_),0),LEX_));

  1505 			d:=mv_deg2(!gcd_r); (* deg(gcd_r,z) *)

  1506 			if (!d < !current_degree) then

  1507 			    (

  1508 			     current_degree:= !d;

  1509 			     interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)

  1510 			     )

  1511 			else

  1512 			    (

  1513 			     if (!d = !current_degree) then

  1514 				 (

  1515 				  interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list)

  1516 				  )

  1517 			     else ()

  1518 				 )

  1519 			    );

  1520 		      if length(!interpolation_list)> uv_mod_min(mv_deg(p1),mv_deg(p2)) then

  1521 			  (

  1522 			   candidate := mv_newton(rev(!interpolation_list));

  1523 			   if !candidate=[] then ()

  1524 			   else

  1525 			       (

  1526 				candidate:=mv_pp(!candidate);

  1527 				if mv_divides(!candidate,p1) andalso mv_divides(!candidate,p2) then

  1528 				    (

  1529 				     gcd:= mv_mul(!candidate,cont,LEX_);

  1530 				     exit:=1

  1531 				     )

  1532 				else ()

  1533 				    );

  1534 			       interpolation_list:=[mv_correct(!gcd_r,0)]

  1535 			       )

  1536 		      else ()

  1537 			  )

  1538 	     );

  1539 	     (!gcd):mv_poly

  1540 	     )

  1541     end;

  1544 (*. calculates the least common divisor of two polynomials .*)

  1545 fun mv_lcm (p1:mv_poly) (p2:mv_poly) :mv_poly =

  1546     (

  1547      #1(mv_division(mv_mul(p1,p2,LEX_),mv_gcd p1 p2,LEX_))

  1548      );

  1550 (* gets the variables (strings) of a term *)

  1552 fun get_vars(term1) = (map free2str) (vars term1); (*["a","b","c"]; *)

  1554 (*. counts the negative coefficents in a polynomial .*)

  1555 fun count_neg ([]:mv_poly) = 0

  1556   | count_neg ((c,e)::xs) = if c<0 then 1+count_neg xs

  1557 			  else count_neg xs;

  1559 (*. help function for is_polynomial

  1560     checks the order of the operators .*)

  1561 fun test_polynomial (Const ("uminus",_) $ Free (str,_)) _ = true (*WN.13.3.03*)

  1562   | test_polynomial (t as Free(str,_)) v = true

  1563   | test_polynomial (t as Const ("Groups.times_class.times",_) $ t1 $ t2) v = if v="^" then false

  1564 						     else (test_polynomial t1 "*") andalso (test_polynomial t2 "*")

  1565   | test_polynomial (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

  1566 							  else (test_polynomial t1 " ") andalso (test_polynomial t2 " ")

  1567   | test_polynomial (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_polynomial t1 "^") andalso (test_polynomial t2 "^")

  1568   | test_polynomial _ v = false;

  1570 (*. tests if a term is a polynomial .*)

  1571 fun is_polynomial t = test_polynomial t " ";

  1573 (*. help function for is_expanded

  1574     checks the order of the operators .*)

  1575 fun test_exp (t as Free(str,_)) v = true

  1576   | test_exp (t as Const ("Groups.times_class.times",_) $ t1 $ t2) v = if v="^" then false

  1577 						     else (test_exp t1 "*") andalso (test_exp t2 "*")

  1578   | test_exp (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

  1579 							  else (test_exp t1 " ") andalso (test_exp t2 " ")

  1580   | test_exp (t as Const ("Groups.minus_class.minus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

  1581 							  else (test_exp t1 " ") andalso (test_exp t2 " ")

  1582   | test_exp (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_exp t1 "^") andalso (test_exp t2 "^")

  1583   | test_exp  _ v = false;

  1586 (*. help function for check_coeff:

  1587     converts the term to a list of coefficients .*)

  1588 fun term2coef' (t as Free(str,_(*typ*))) v :mv_poly option =

  1589     let

  1590 	val x= Unsynchronized.ref  NONE;

  1591 	val len= Unsynchronized.ref  0;

  1592 	val vl= Unsynchronized.ref  [];

  1593 	val vh= Unsynchronized.ref  [];

  1594 	val i= Unsynchronized.ref  0;

  1595     in

  1596 	if is_numeral str then

  1597 	    (

  1598 	     SOME [(((the o int_of_str) str),mv_null2(v))] handle _ => NONE

  1599 		 )

  1600 	else (* variable *)

  1601 	    (

  1602 	     len:=length(v);

  1603 	     vh:=v;

  1604 	     while ((!len)>(!i)) do

  1605 		 (

  1606 		  if str=hd((!vh)) then

  1607 		      (

  1608 		       vl:=1::(!vl)

  1609 		       )

  1610 		  else

  1611 		      (

  1612 		       vl:=0::(!vl)

  1613 		       );

  1614 		      vh:=tl(!vh);

  1615 		      i:=(!i)+1

  1616 		      );

  1617 		 SOME [(1,rev(!vl))] handle _ => NONE

  1618 	    )

  1619     end

  1620   | term2coef' (Const ("Groups.times_class.times",_) $ t1 $ t2) v :mv_poly option=

  1621     let

  1622 	val t1pp= Unsynchronized.ref  [];

  1623 	val t2pp= Unsynchronized.ref  [];

  1624 	val t1c= Unsynchronized.ref  0;

  1625 	val t2c= Unsynchronized.ref  0;

  1626     in

  1627 	(

  1628 	 t1pp:=(#2(hd(the(term2coef' t1 v))));

  1629 	 t2pp:=(#2(hd(the(term2coef' t2 v))));

  1630 	 t1c:=(#1(hd(the(term2coef' t1 v))));

  1631 	 t2c:=(#1(hd(the(term2coef' t2 v))));

  1633 	 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] handle _ => NONE

  1635 	 )

  1636     end

  1637   | term2coef' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ (t2 as Free (str2,_))) v :mv_poly option=

  1638     let

  1639 	val x= Unsynchronized.ref  NONE;

  1640 	val len= Unsynchronized.ref  0;

  1641 	val vl= Unsynchronized.ref  [];

  1642 	val vh= Unsynchronized.ref  [];

  1643 	val vtemp= Unsynchronized.ref  [];

  1644 	val i= Unsynchronized.ref  0;

  1645     in

  1646     (

  1647      if (not o is_numeral) str1 andalso is_numeral str2 then

  1648 	 (

  1649 	  len:=length(v);

  1650 	  vh:=v;

  1652 	  while ((!len)>(!i)) do

  1653 	      (

  1654 	       if str1=hd((!vh)) then

  1655 		   (

  1656 		    vl:=((the o int_of_str) str2)::(!vl)

  1657 		    )

  1658 	       else

  1659 		   (

  1660 		    vl:=0::(!vl)

  1661 		    );

  1662 		   vh:=tl(!vh);

  1663 		   i:=(!i)+1

  1664 		   );

  1665 	      SOME [(1,rev(!vl))] handle _ => NONE

  1666 	      )

  1667      else error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term")

  1668 	 )

  1669     end

  1670   | term2coef' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option=

  1671     (

  1672      SOME ((the(term2coef' t1 v)) @ (the(term2coef' t2 v))) handle _ => NONE

  1673 	 )

  1674   | term2coef' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option=

  1675     (

  1676      SOME ((the(term2coef' t1 v)) @ mv_skalar_mul((the(term2coef' t2 v)),1)) handle _ => NONE

  1677 	 )

  1678   | term2coef' (term) v = error ("RATIONALS_TERM2COEF_EXCEPTION 2: Invalid term");

  1680 (*. checks if all coefficients of a polynomial are positiv (except the first) .*)

  1681 fun check_coeff t = (* erste Koeffizient kann <0 sein !!! *)

  1682     if count_neg(tl(the(term2coef' t (get_vars(t)))))=0 then true

  1683     else false;

  1685 (*. checks for expanded term [3] .*)

  1686 fun is_expanded t = test_exp t " " andalso check_coeff(t);

  1688 (*WN.7.3.03 Hilfsfunktion f"ur term2poly'*)

  1689 fun mk_monom v' p vs =

  1690     let fun conv p (v: string) = if v'= v then p else 0

  1691     in map (conv p) vs end;

  1692 (* mk_monom "y" 5 ["a","b","x","y","z"];

  1693 val it = [0,0,0,5,0] : int list*)

  1695 (*. this function converts the term representation into the internal representation mv_poly .*)

  1696 fun term2poly' (Const ("uminus",_) $ Free (str,_)) v = (*WN.7.3.03*)

  1697     if is_numeral str

  1698     then SOME [((the o int_of_str) ("-"^str), mk_monom "#" 0 v)]

  1699     else SOME [(~1, mk_monom str 1 v)]

  1701   | term2poly' (Free(str,_)) v :mv_poly option =

  1702     let

  1703 	val x= Unsynchronized.ref  NONE;

  1704 	val len= Unsynchronized.ref  0;

  1705 	val vl= Unsynchronized.ref  [];

  1706 	val vh= Unsynchronized.ref  [];

  1707 	val i= Unsynchronized.ref  0;

  1708     in

  1709 	if is_numeral str then

  1710 	    (

  1711 	     SOME [(((the o int_of_str) str),mv_null2 v)] handle _ => NONE

  1712 		 )

  1713 	else (* variable *)

  1714 	    (

  1715 	     len:=length v;

  1716 	     vh:= v;

  1717 	     while ((!len)>(!i)) do

  1718 		 (

  1719 		  if str=hd((!vh)) then

  1720 		      (

  1721 		       vl:=1::(!vl)

  1722 		       )

  1723 		  else

  1724 		      (

  1725 		       vl:=0::(!vl)

  1726 		       );

  1727 		      vh:=tl(!vh);

  1728 		      i:=(!i)+1

  1729 		      );

  1730 		 SOME [(1,rev(!vl))] handle _ => NONE

  1731 	    )

  1732     end

  1733   | term2poly' (Const ("Groups.times_class.times",_) $ t1 $ t2) v :mv_poly option=

  1734     let

  1735 	val t1pp= Unsynchronized.ref  [];

  1736 	val t2pp= Unsynchronized.ref  [];

  1737 	val t1c= Unsynchronized.ref  0;

  1738 	val t2c= Unsynchronized.ref  0;

  1739     in

  1740 	(

  1741 	 t1pp:=(#2(hd(the(term2poly' t1 v))));

  1742 	 t2pp:=(#2(hd(the(term2poly' t2 v))));

  1743 	 t1c:=(#1(hd(the(term2poly' t1 v))));

  1744 	 t2c:=(#1(hd(the(term2poly' t2 v))));

  1746 	 SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )]

  1747 	 handle _ => NONE

  1749 	 )

  1750     end

  1751   | term2poly' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $

  1752 		      (t2 as Free (str2,_))) v :mv_poly option=

  1753     let

  1754 	val x= Unsynchronized.ref  NONE;

  1755 	val len= Unsynchronized.ref  0;

  1756 	val vl= Unsynchronized.ref  [];

  1757 	val vh= Unsynchronized.ref  [];

  1758 	val vtemp= Unsynchronized.ref  [];

  1759 	val i= Unsynchronized.ref  0;

  1760     in

  1761     (

  1762      if (not o is_numeral) str1 andalso is_numeral str2 then

  1763 	 (

  1764 	  len:=length(v);

  1765 	  vh:=v;

  1767 	  while ((!len)>(!i)) do

  1768 	      (

  1769 	       if str1=hd((!vh)) then

  1770 		   (

  1771 		    vl:=((the o int_of_str) str2)::(!vl)

  1772 		    )

  1773 	       else

  1774 		   (

  1775 		    vl:=0::(!vl)

  1776 		    );

  1777 		   vh:=tl(!vh);

  1778 		   i:=(!i)+1

  1779 		   );

  1780 	      SOME [(1,rev(!vl))] handle _ => NONE

  1781 	      )

  1782      else error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term")

  1783 	 )

  1784     end

  1785   | term2poly' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option =

  1786     (

  1787      SOME ((the(term2poly' t1 v)) @ (the(term2poly' t2 v))) handle _ => NONE

  1788 	 )

  1789   | term2poly' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option =

  1790     (

  1791      SOME ((the(term2poly' t1 v)) @ mv_skalar_mul((the(term2poly' t2 v)),~1)) handle _ => NONE

  1792 	 )

  1793   | term2poly' (term) v = error ("RATIONALS_TERM2POLY_EXCEPTION 2: Invalid term");

  1795 (*. translates an Isabelle term into internal representation.

  1796     term2poly

  1797     fn : term ->              (*normalform [2]                    *)

  1798     	 string list ->       (*for ...!!! BITTE DIE ERKLÄRUNG,

  1799     			       DIE DU MIR LETZTES MAL GEGEBEN HAST*)

  1800     	 mv_monom list        (*internal representation           *)

  1801     		  option      (*the translation may fail with NONE*)

  1802 .*)

  1803 fun term2poly (t:term) v =

  1804      if is_polynomial t then term2poly' t v

  1805      else error ("term2poly: invalid = "^(term2str t));

  1807 (*. same as term2poly with automatic detection of the variables .*)

  1808 fun term2polyx t = term2poly t (((map free2str) o vars) t);

  1810 (*. checks if the term is in expanded polynomial form and converts it into the internal representation .*)

  1811 fun expanded2poly (t:term) v =

  1812     (*if is_expanded t then*) term2poly' t v

  1813     (*else error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*);

  1815 (*. same as expanded2poly with automatic detection of the variables .*)

  1816 fun expanded2polyx t = expanded2poly t (((map free2str) o vars) t);

  1818 (*. converts a powerproduct into term representation .*)

  1819 fun powerproduct2term(xs,v) =

  1820     let

  1821 	val xss= Unsynchronized.ref  xs;

  1822 	val vv= Unsynchronized.ref  v;

  1823     in

  1824 	(

  1825 	 while hd(!xss)=0 do

  1826 	     (

  1827 	      xss:=tl(!xss);

  1828 	      vv:=tl(!vv)

  1829 	      );

  1831 	 if list_is_null(tl(!xss)) then

  1832 	     (

  1833 	      if hd(!xss)=1 then Free(hd(!vv), HOLogic.realT)

  1834 	      else

  1835 		  (

  1836 		   Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1837 		   Free(hd(!vv), HOLogic.realT) $  Free(str_of_int (hd(!xss)),HOLogic.realT)

  1838 		   )

  1839 	      )

  1840 	 else

  1841 	     (

  1842 	      if hd(!xss)=1 then

  1843 		  (

  1844 		   Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1845 		   Free(hd(!vv), HOLogic.realT) $

  1846 		   powerproduct2term(tl(!xss),tl(!vv))

  1847 		   )

  1848 	      else

  1849 		  (

  1850 		   Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1851 		   (

  1852 		    Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1853 		    Free(hd(!vv), HOLogic.realT) $  Free(str_of_int (hd(!xss)),HOLogic.realT)

  1854 		    ) $

  1855 		    powerproduct2term(tl(!xss),tl(!vv))

  1856 		   )

  1857 	      )

  1858 	 )

  1859     end;

  1861 (*. converts a monom into term representation .*)

  1862 (*fun monom2term ((c,e):mv_monom, v:string list) =

  1863     if c=0 then Free(str_of_int 0,HOLogic.realT)

  1864     else

  1865 	(

  1866 	 if list_is_null(e) then

  1867 	     (

  1868 	      Free(str_of_int c,HOLogic.realT)

  1869 	      )

  1870 	 else

  1871 	     (

  1872 	      if c=1 then

  1873 		  (

  1874 		   powerproduct2term(e,v)

  1875 		   )

  1876 	      else

  1877 		  (

  1878 		   Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1879 		   Free(str_of_int c,HOLogic.realT)  $

  1880 		   powerproduct2term(e,v)

  1881 		   )

  1882 		  )

  1883 	     );*)

  1886 (*fun monom2term ((i, is):mv_monom, v) =

  1887     if list_is_null is

  1888     then

  1889 	if i >= 0

  1890 	then Free (str_of_int i, HOLogic.realT)

  1891 	else Const ("uminus", HOLogic.realT --> HOLogic.realT) $

  1892 		   Free ((str_of_int o abs) i, HOLogic.realT)

  1893     else

  1894 	if i > 0

  1895 	then Const ("Groups.times_class.times", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $

  1896 		   (Free (str_of_int i, HOLogic.realT)) $

  1897 		   powerproduct2term(is, v)

  1898 	else Const ("Groups.times_class.times", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $

  1899 		   (Const ("uminus", HOLogic.realT --> HOLogic.realT) $

  1900 		   Free ((str_of_int o abs) i, HOLogic.realT)) $

  1901 		   powerproduct2term(is, vs);---------------------------*)

  1902 fun monom2term ((i, is) : mv_monom, vs) =

  1903     if list_is_null is

  1904     then Free (str_of_int i, HOLogic.realT)

  1905     else if i = 1

  1906     then powerproduct2term (is, vs)

  1907     else Const ("Groups.times_class.times", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $

  1908 	       (Free (str_of_int i, HOLogic.realT)) $

  1909 	       powerproduct2term (is, vs);

  1911 (*. converts the internal polynomial representation into an Isabelle term.*)

  1912 fun poly2term' ([] : mv_poly, vs) = Free(str_of_int 0, HOLogic.realT)

  1913   | poly2term' ([(c, e) : mv_monom], vs) = monom2term ((c, e), vs)

  1914   | poly2term' ((c, e) :: ces, vs) =

  1915     Const("Groups.plus_class.plus", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $

  1916          poly2term (ces, vs) $ monom2term ((c, e), vs)

  1917 and poly2term (xs, vs) = poly2term' (rev (sort (mv_geq LEX_) (xs)), vs);

  1920 (*. converts a monom into term representation .*)

  1921 (*. ignores the sign of the coefficients => use only for exp-poly functions .*)

  1922 fun monom2term2((c,e):mv_monom, v:string list) =

  1923     if c=0 then Free(str_of_int 0,HOLogic.realT)

  1924     else

  1925 	(

  1926 	 if list_is_null(e) then

  1927 	     (

  1928 	      Free(str_of_int (abs(c)),HOLogic.realT)

  1929 	      )

  1930 	 else

  1931 	     (

  1932 	      if abs(c)=1 then

  1933 		  (

  1934 		   powerproduct2term(e,v)

  1935 		   )

  1936 	      else

  1937 		  (

  1938 		   Const("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1939 		   Free(str_of_int (abs(c)),HOLogic.realT)  $

  1940 		   powerproduct2term(e,v)

  1941 		   )

  1942 		  )

  1943 	     );

  1945 (*. converts the expanded polynomial representation into the term representation .*)

  1946 fun exp2term' ([]:mv_poly,vars) =  Free(str_of_int 0,HOLogic.realT)

  1947   | exp2term' ([(c,e)],vars) =     monom2term((c,e),vars)

  1948   | exp2term' ((c1,e1)::others,vars) =

  1949     if c1<0 then

  1950 	Const("Groups.minus_class.minus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1951 	exp2term'(others,vars) $

  1952 	(

  1953 	 monom2term2((c1,e1),vars)

  1954 	 )

  1955     else

  1956 	Const("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1957 	exp2term'(others,vars) $

  1958 	(

  1959 	 monom2term2((c1,e1),vars)

  1960 	 );

  1962 (*. sorts the powerproduct by lexicographic termorder and converts them into

  1963     a term in polynomial representation .*)

  1964 fun poly2expanded (xs,vars) = exp2term'(rev(sort (mv_geq LEX_) (xs)),vars);

  1966 (*. converts a polynomial into expanded form .*)

  1967 fun polynomial2expanded t =

  1968     (let

  1969 	val vars=(((map free2str) o vars) t);

  1970     in

  1971 	SOME (poly2expanded (the (term2poly t vars), vars))

  1972     end) handle _ => NONE;

  1974 (*. converts a polynomial into polynomial form .*)

  1975 fun expanded2polynomial t =

  1976     (let

  1977 	val vars=(((map free2str) o vars) t);

  1978     in

  1979 	SOME (poly2term (the (expanded2poly t vars), vars))

  1980     end) handle _ => NONE;

  1983 (*. calculates the greatest common divisor of numerator and denominator and seperates it from each .*)

  1984 fun step_cancel (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =

  1985     let

  1986 	val p1' = Unsynchronized.ref [];

  1987 	val p2' = Unsynchronized.ref [];

  1988 	val p3  = Unsynchronized.ref []

  1989 	val vars = rev(get_vars(p1) union get_vars(p2));

  1990     in

  1991 	(

  1992          p1':= sort (mv_geq LEX_) (the (term2poly p1 vars ));

  1993        	 p2':= sort (mv_geq LEX_) (the (term2poly p2 vars ));

  1994 	 p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));

  1995 	 if (!p3)=[(1,mv_null2(vars))] then

  1996 	     (

  1997 	      Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2

  1998 	      )

  1999 	 else

  2000 	     (

  2002 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));

  2003 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));

  2005 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then

  2006 	      (

  2007 	       Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2008 	       $

  2009 	       (

  2010 		Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2011 		poly2term(!p1',vars) $

  2012 		poly2term(!p3,vars)

  2013 		)

  2014 	       $

  2015 	       (

  2016 		Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2017 		poly2term(!p2',vars) $

  2018 		poly2term(!p3,vars)

  2019 		)

  2020 	       )

  2021 	      else

  2022 	      (

  2023 	       p1':=mv_skalar_mul(!p1',~1);

  2024 	       p2':=mv_skalar_mul(!p2',~1);

  2025 	       p3:=mv_skalar_mul(!p3,~1);

  2026 	       (

  2027 		Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2028 		$

  2029 		(

  2030 		 Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2031 		 poly2term(!p1',vars) $

  2032 		 poly2term(!p3,vars)

  2033 		 )

  2034 		$

  2035 		(

  2036 		 Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2037 		 poly2term(!p2',vars) $

  2038 		 poly2term(!p3,vars)

  2039 		 )

  2040 		)

  2041 	       )

  2042 	      )

  2043 	     )

  2044     end

  2045 | step_cancel _ = error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction");

  2047 (*. calculates the greatest common divisor of numerator and denominator and divides each through it .*)

  2048 fun direct_cancel (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =

  2049     let

  2050 	val p1' = Unsynchronized.ref [];

  2051 	val p2' = Unsynchronized.ref [];

  2052 	val p3  = Unsynchronized.ref []

  2053 	val vars = rev(get_vars(p1) union get_vars(p2));

  2054     in

  2055 	(

  2056 	 p1':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p1 vars )),LEX_));

  2057 	 p2':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p2 vars )),LEX_));

  2058 	 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));

  2060 	 if (!p3)=[(1,mv_null2(vars))] then

  2061 	     (

  2062 	      (Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])

  2063 	      )

  2064 	 else

  2065 	     (

  2066 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));

  2067 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));

  2068 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then

  2069 	      (

  2070 	       (

  2071 		Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2072 		$

  2073 		(

  2074 		 poly2term((!p1'),vars)

  2075 		 )

  2076 		$

  2077 		(

  2078 		 poly2term((!p2'),vars)

  2079 		 )

  2080 		)

  2081 	       ,

  2082 	       if mv_grad(!p3)>0 then

  2083 		   [

  2084 		    (

  2085 		     Const ("HOL.Not",[bool]--->bool) $

  2086 		     (

  2087 		      Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $

  2088 		      poly2term((!p3),vars) $

  2089 		      Free("0",HOLogic.realT)

  2090 		      )

  2091 		     )

  2092 		    ]

  2093 	       else

  2094 		   []

  2095 		   )

  2096 	      else

  2097 		  (

  2098 		   p1':=mv_skalar_mul(!p1',~1);

  2099 		   p2':=mv_skalar_mul(!p2',~1);

  2100 		   if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);

  2101 		       (

  2102 			Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2103 			$

  2104 			(

  2105 			 poly2term((!p1'),vars)

  2106 			 )

  2107 			$

  2108 			(

  2109 			 poly2term((!p2'),vars)

  2110 			 )

  2111 			,

  2112 			if mv_grad(!p3)>0 then

  2113 			    [

  2114 			     (

  2115 			      Const ("HOL.Not",[bool]--->bool) $

  2116 			      (

  2117 			       Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $

  2118 			       poly2term((!p3),vars) $

  2119 			       Free("0",HOLogic.realT)

  2120 			       )

  2121 			      )

  2122 			     ]

  2123 			else

  2124 			    []

  2125 			    )

  2126 		       )

  2127 		  )

  2128 	     )

  2129     end

  2130   | direct_cancel _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");

  2132 (*. same es direct_cancel, this time for expanded forms (input+output).*)

  2133 fun direct_cancel_expanded (t as Const ("Fields.inverse_class.divide",_) $ p1 $ p2) =

  2134     let

  2135 	val p1' = Unsynchronized.ref [];

  2136 	val p2' = Unsynchronized.ref [];

  2137 	val p3  = Unsynchronized.ref []

  2138 	val vars = rev(get_vars(p1) union get_vars(p2));

  2139     in

  2140 	(

  2141 	 p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_));

  2142 	 p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_));

  2143 	 p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2'));

  2145 	 if (!p3)=[(1,mv_null2(vars))] then

  2146 	     (

  2147 	      (Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])

  2148 	      )

  2149 	 else

  2150 	     (

  2151 	      p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_)));

  2152 	      p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_)));

  2153 	      if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then

  2154 	      (

  2155 	       (

  2156 		Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2157 		$

  2158 		(

  2159 		 poly2expanded((!p1'),vars)

  2160 		 )

  2161 		$

  2162 		(

  2163 		 poly2expanded((!p2'),vars)

  2164 		 )

  2165 		)

  2166 	       ,

  2167 	       if mv_grad(!p3)>0 then

  2168 		   [

  2169 		    (

  2170 		     Const ("HOL.Not",[bool]--->bool) $

  2171 		     (

  2172 		      Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $

  2173 		      poly2expanded((!p3),vars) $

  2174 		      Free("0",HOLogic.realT)

  2175 		      )

  2176 		     )

  2177 		    ]

  2178 	       else

  2179 		   []

  2180 		   )

  2181 	      else

  2182 		  (

  2183 		   p1':=mv_skalar_mul(!p1',~1);

  2184 		   p2':=mv_skalar_mul(!p2',~1);

  2185 		   if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1);

  2186 		       (

  2187 			Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2188 			$

  2189 			(

  2190 			 poly2expanded((!p1'),vars)

  2191 			 )

  2192 			$

  2193 			(

  2194 			 poly2expanded((!p2'),vars)

  2195 			 )

  2196 			,

  2197 			if mv_grad(!p3)>0 then

  2198 			    [

  2199 			     (

  2200 			      Const ("HOL.Not",[bool]--->bool) $

  2201 			      (

  2202 			       Const("HOL.eq",[HOLogic.realT,HOLogic.realT]--->bool) $

  2203 			       poly2expanded((!p3),vars) $

  2204 			       Free("0",HOLogic.realT)

  2205 			       )

  2206 			      )

  2207 			     ]

  2208 			else

  2209 			    []

  2210 			    )

  2211 		       )

  2212 		  )

  2213 	     )

  2214     end

  2215   | direct_cancel_expanded _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");

  2218 (*. adds two fractions .*)

  2219 fun add_fract ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =

  2220     let

  2221 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);

  2222 	val t11'= Unsynchronized.ref  (the(term2poly t11 vars));

  2223 (* stopped Test_Isac.thy ...

  2224 val _= tracing"### add_fract: done t11"

  2225 *)

  2226 	val t12'= Unsynchronized.ref  (the(term2poly t12 vars));

  2227 (* stopped Test_Isac.thy ...

  2228 val _= tracing"### add_fract: done t12"

  2229 *)

  2230 	val t21'= Unsynchronized.ref  (the(term2poly t21 vars));

  2231 (* stopped Test_Isac.thy ...

  2232 val _= tracing"### add_fract: done t21"

  2233 *)

  2234 	val t22'= Unsynchronized.ref  (the(term2poly t22 vars));

  2235 (* stopped Test_Isac.thy ...

  2236 val _= tracing"### add_fract: done t22"

  2237 *)

  2238 	val den= Unsynchronized.ref  [];

  2239 	val nom= Unsynchronized.ref  [];

  2240 	val m1= Unsynchronized.ref  [];

  2241 	val m2= Unsynchronized.ref  [];

  2242     in

  2244 	(

  2245 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));

  2246 tracing"### add_fract: done sort mv_lcm";

  2247 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));

  2248 tracing"### add_fract: done sort mv_division t12";

  2249 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));

  2250 tracing"### add_fract: done sort mv_division t22";

  2251 	 nom :=sort (mv_geq LEX_)

  2252 		    (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),

  2253 				       mv_mul(!t21',!m2,LEX_),

  2254 				       LEX_),

  2255 				LEX_));

  2256 tracing"### add_fract: done sort mv_add";

  2257 	 (

  2258 	  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2259 	  $

  2260 	  (

  2261 	   poly2term((!nom),vars)

  2262 	   )

  2263 	  $

  2264 	  (

  2265 	   poly2term((!den),vars)

  2266 	   )

  2267 	  )

  2268 	 )

  2269     end

  2270   | add_fract (_,_) = error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call");

  2272 (*. adds two expanded fractions .*)

  2273 fun add_fract_exp ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =

  2274     let

  2275 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);

  2276 	val t11'= Unsynchronized.ref  (the(expanded2poly t11 vars));

  2277 	val t12'= Unsynchronized.ref  (the(expanded2poly t12 vars));

  2278 	val t21'= Unsynchronized.ref  (the(expanded2poly t21 vars));

  2279 	val t22'= Unsynchronized.ref  (the(expanded2poly t22 vars));

  2280 	val den= Unsynchronized.ref  [];

  2281 	val nom= Unsynchronized.ref  [];

  2282 	val m1= Unsynchronized.ref  [];

  2283 	val m2= Unsynchronized.ref  [];

  2284     in

  2286 	(

  2287 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));

  2288 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));

  2289 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));

  2290 	 nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_));

  2291 	 (

  2292 	  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2293 	  $

  2294 	  (

  2295 	   poly2expanded((!nom),vars)

  2296 	   )

  2297 	  $

  2298 	  (

  2299 	   poly2expanded((!den),vars)

  2300 	   )

  2301 	  )

  2302 	 )

  2303     end

  2304   | add_fract_exp (_,_) = error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call");

  2306 (*. adds a list of terms .*)

  2307 fun add_list_of_fractions []= (Free("0",HOLogic.realT),[])

  2308   | add_list_of_fractions [x]= direct_cancel x

  2309   | add_list_of_fractions (x::y::xs) =

  2310     let

  2311 	val (t1a,rest1)=direct_cancel(x);

  2312 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(x)";

  2313 	val (t2a,rest2)=direct_cancel(y);

  2314 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(y)";

  2315 	val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs));

  2316 val _= tracing"### add_list_of_fractions xs: has done add_list_of_fraction xs";

  2317 	val (t4a,rest4)=direct_cancel(t3a);

  2318 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(t3a)";

  2319 	val rest=rest1 union rest2 union rest3 union rest4;

  2320     in

  2321 	(tracing"### add_list_of_fractions in";

  2322 	 (

  2323 	 (t4a,rest)

  2324 	 )

  2325 	 )

  2326     end;

  2328 (*. adds a list of expanded terms .*)

  2329 fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])

  2330   | add_list_of_fractions_exp [x]= direct_cancel_expanded x

  2331   | add_list_of_fractions_exp (x::y::xs) =

  2332     let

  2333 	val (t1a,rest1)=direct_cancel_expanded(x);

  2334 	val (t2a,rest2)=direct_cancel_expanded(y);

  2335 	val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));

  2336 	val (t4a,rest4)=direct_cancel_expanded(t3a);

  2337 	val rest=rest1 union rest2 union rest3 union rest4;

  2338     in

  2339 	(

  2340 	 (t4a,rest)

  2341 	 )

  2342     end;

  2344 (*. calculates the lcm of a list of mv_poly .*)

  2345 fun calc_lcm ([x],var)= (x,var)

  2346   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2348 (*. converts a list of terms to a list of mv_poly .*)

  2349 fun t2d([],_)=[]

  2350   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2352 (*. same as t2d, this time for expanded forms .*)

  2353 fun t2d_exp([],_)=[]

  2354   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2356 (*. converts a list of fract terms to a list of their denominators .*)

  2357 fun termlist2denominators [] = ([],[])

  2358   | termlist2denominators xs =

  2359     let

  2360 	val xxs= Unsynchronized.ref  xs;

  2361 	val var= Unsynchronized.ref  [];

  2362     in

  2363 	var:=[];

  2364 	while length(!xxs)>0 do

  2365 	    (

  2366 	     let

  2367 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2368 	     in

  2369 		 (

  2370 		  xxs:=tl(!xxs);

  2371 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2372 		  )

  2373 	     end

  2374 	     );

  2375 	    (t2d(xs,!var),!var)

  2376     end;

  2378 (*. calculates the lcm of a list of mv_poly .*)

  2379 fun calc_lcm ([x],var)= (x,var)

  2380   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2382 (*. converts a list of terms to a list of mv_poly .*)

  2383 fun t2d([],_)=[]

  2384   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2386 (*. same as t2d, this time for expanded forms .*)

  2387 fun t2d_exp([],_)=[]

  2388   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2390 (*. converts a list of fract terms to a list of their denominators .*)

  2391 fun termlist2denominators [] = ([],[])

  2392   | termlist2denominators xs =

  2393     let

  2394 	val xxs= Unsynchronized.ref  xs;

  2395 	val var= Unsynchronized.ref  [];

  2396     in

  2397 	var:=[];

  2398 	while length(!xxs)>0 do

  2399 	    (

  2400 	     let

  2401 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2402 	     in

  2403 		 (

  2404 		  xxs:=tl(!xxs);

  2405 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2406 		  )

  2407 	     end

  2408 	     );

  2409 	    (t2d(xs,!var),!var)

  2410     end;

  2412 (*. same as termlist2denminators, this time for expanded forms .*)

  2413 fun termlist2denominators_exp [] = ([],[])

  2414   | termlist2denominators_exp xs =

  2415     let

  2416 	val xxs= Unsynchronized.ref  xs;

  2417 	val var= Unsynchronized.ref  [];

  2418     in

  2419 	var:=[];

  2420 	while length(!xxs)>0 do

  2421 	    (

  2422 	     let

  2423 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2424 	     in

  2425 		 (

  2426 		  xxs:=tl(!xxs);

  2427 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2428 		  )

  2429 	     end

  2430 	     );

  2431 	    (t2d_exp(xs,!var),!var)

  2432     end;

  2434 (*. reduces all fractions to the least common denominator .*)

  2435 fun com_den(x::xs,denom,den,var)=

  2436     let

  2437 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2438 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2439 	val p3= #1(mv_division(denom,p2,LEX_));

  2440 	val p1var=get_vars(p1');

  2441     in

  2442 	if length(xs)>0 then

  2443 	    if p3=[(1,mv_null2(var))] then

  2444 		(

  2445 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2446 		 $

  2447 		 (

  2448 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2449 		  $

  2450 		  poly2term(the (term2poly p1' p1var),p1var)

  2451 		  $

  2452 		  den

  2453 		  )

  2454 		 $

  2455 		 #1(com_den(xs,denom,den,var))

  2456 		,

  2457 		[]

  2458 		)

  2459 	    else

  2460 		(

  2461 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2462 		 $

  2463 		 (

  2464 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2465 		  $

  2466 		  (

  2467 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2468 		   poly2term(the (term2poly p1' p1var),p1var) $

  2469 		   poly2term(p3,var)

  2470 		   )

  2471 		  $

  2472 		  (

  2473 		   den

  2474 		   )

  2475 		  )

  2476 		 $

  2477 		 #1(com_den(xs,denom,den,var))

  2478 		,

  2479 		[]

  2480 		)

  2481 	else

  2482 	    if p3=[(1,mv_null2(var))] then

  2483 		(

  2484 		 (

  2485 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2486 		  $

  2487 		  poly2term(the (term2poly p1' p1var),p1var)

  2488 		  $

  2489 		  den

  2490 		  )

  2491 		 ,

  2492 		 []

  2493 		 )

  2494 	     else

  2495 		 (

  2496 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2497 		  $

  2498 		  (

  2499 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2500 		   poly2term(the (term2poly p1' p1var),p1var) $

  2501 		   poly2term(p3,var)

  2502 		   )

  2503 		  $

  2504 		  den

  2505 		  ,

  2506 		  []

  2507 		  )

  2508     end;

  2510 (*. same as com_den, this time for expanded forms .*)

  2511 fun com_den_exp(x::xs,denom,den,var)=

  2512     let

  2513 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2514 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2515 	val p3= #1(mv_division(denom,p2,LEX_));

  2516 	val p1var=get_vars(p1');

  2517     in

  2518 	if length(xs)>0 then

  2519 	    if p3=[(1,mv_null2(var))] then

  2520 		(

  2521 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2522 		 $

  2523 		 (

  2524 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2525 		  $

  2526 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2527 		  $

  2528 		  den

  2529 		  )

  2530 		 $

  2531 		 #1(com_den_exp(xs,denom,den,var))

  2532 		,

  2533 		[]

  2534 		)

  2535 	    else

  2536 		(

  2537 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2538 		 $

  2539 		 (

  2540 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2541 		  $

  2542 		  (

  2543 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2544 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2545 		   poly2expanded(p3,var)

  2546 		   )

  2547 		  $

  2548 		  (

  2549 		   den

  2550 		   )

  2551 		  )

  2552 		 $

  2553 		 #1(com_den_exp(xs,denom,den,var))

  2554 		,

  2555 		[]

  2556 		)

  2557 	else

  2558 	    if p3=[(1,mv_null2(var))] then

  2559 		(

  2560 		 (

  2561 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2562 		  $

  2563 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2564 		  $

  2565 		  den

  2566 		  )

  2567 		 ,

  2568 		 []

  2569 		 )

  2570 	     else

  2571 		 (

  2572 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2573 		  $

  2574 		  (

  2575 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2576 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2577 		   poly2expanded(p3,var)

  2578 		   )

  2579 		  $

  2580 		  den

  2581 		  ,

  2582 		  []

  2583 		  )

  2584     end;

  2586 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2587 -------------------------------------------------------------

  2588 (* WN0210???SK brauch ma des überhaupt *)

  2589 fun com_den2(x::xs,denom,den,var)=

  2590     let

  2591 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2592 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2593 	val p3= #1(mv_division(denom,p2,LEX_));

  2594 	val p1var=get_vars(p1');

  2595     in

  2596 	if length(xs)>0 then

  2597 	    if p3=[(1,mv_null2(var))] then

  2598 		(

  2599 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2600 		 poly2term(the(term2poly p1' p1var),p1var) $

  2601 		 com_den2(xs,denom,den,var)

  2602 		)

  2603 	    else

  2604 		(

  2605 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2606 		 (

  2607 		   let

  2608 		       val p3'=poly2term(p3,var);

  2609 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2610 		   in

  2611 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2612 		   end

  2613 		  ) $

  2614 		 com_den2(xs,denom,den,var)

  2615 		)

  2616 	else

  2617 	    if p3=[(1,mv_null2(var))] then

  2618 		(

  2619 		 poly2term(the(term2poly p1' p1var),p1var)

  2620 		 )

  2621 	     else

  2622 		 (

  2623 		   let

  2624 		       val p3'=poly2term(p3,var);

  2625 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2626 		   in

  2627 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2628 		   end

  2629 		  )

  2630     end;

  2632 (* WN0210???SK brauch ma des überhaupt *)

  2633 fun com_den_exp2(x::xs,denom,den,var)=

  2634     let

  2635 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2636 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2637 	val p3= #1(mv_division(denom,p2,LEX_));

  2638 	val p1var=get_vars p1';

  2639     in

  2640 	if length(xs)>0 then

  2641 	    if p3=[(1,mv_null2(var))] then

  2642 		(

  2643 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2644 		 poly2expanded(the (expanded2poly p1' p1var),p1var) $

  2645 		 com_den_exp2(xs,denom,den,var)

  2646 		)

  2647 	    else

  2648 		(

  2649 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2650 		 (

  2651 		   let

  2652 		       val p3'=poly2expanded(p3,var);

  2653 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2654 		   in

  2655 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2656 		   end

  2657 		  ) $

  2658 		 com_den_exp2(xs,denom,den,var)

  2659 		)

  2660 	else

  2661 	    if p3=[(1,mv_null2(var))] then

  2662 		(

  2663 		 poly2expanded(the (expanded2poly p1' p1var),p1var)

  2664 		 )

  2665 	     else

  2666 		 (

  2667 		   let

  2668 		       val p3'=poly2expanded(p3,var);

  2669 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2670 		   in

  2671 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2672 		   end

  2673 		  )

  2674     end;

  2675 ---------------------------------------------------------*)

  2678 (*. searches for an element y of a list ys, which has an gcd not 1 with x .*)

  2679 fun exists_gcd (x,[]) = false

  2680   | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then  exists_gcd (x,ys)

  2681 			   else true;

  2683 (*. divides each element of the list xs with y .*)

  2684 fun list_div ([],y) = []

  2685   | list_div (x::xs,y) =

  2686     let

  2687 	val (d,r)=mv_division(x,y,LEX_);

  2688     in

  2689 	if r=[] then

  2690 	    d::list_div(xs,y)

  2691 	else x::list_div(xs,y)

  2692     end;

  2694 (*. checks if x is in the list ys .*)

  2695 fun in_list (x,[]) = false

  2696   | in_list (x,y::ys) = if x=y then true

  2697 			else in_list(x,ys);

  2699 (*. deletes all equal elements of the list xs .*)

  2700 fun kill_equal [] = []

  2701   | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)

  2702 			 else x::kill_equal(xs);

  2704 (*. searches for new factors .*)

  2705 fun new_factors [] = []

  2706   | new_factors (list:mv_poly list):mv_poly list =

  2707     let

  2708 	val l = kill_equal list;

  2709 	val len = length(l);

  2710     in

  2711 	if len>=2 then

  2712 	    (

  2713 	     let

  2714 		 val x::y::xs=l;

  2715 		 val gcd=mv_gcd x y;

  2716 	     in

  2717 		 if gcd=[(1,mv_null2(#2(hd(x))))] then

  2718 		     (

  2719 		      if exists_gcd(x,xs) then new_factors (y::xs @ [x])

  2720 		      else x::new_factors(y::xs)

  2721 	             )

  2722 		 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))

  2723 	     end

  2724 	     )

  2725 	else

  2726 	    if len=1 then [hd(l)]

  2727 	    else []

  2728     end;

  2730 (*. gets the factors of a list .*)

  2731 fun get_factors x = new_factors x;

  2733 (*. multiplies the elements of the list .*)

  2734 fun multi_list [] = []

  2735   | multi_list (x::xs) = if xs=[] then x

  2736 			 else mv_mul(x,multi_list xs,LEX_);

  2738 (*. makes a term out of the elements of the list (polynomial representation) .*)

  2739 fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2740   | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)

  2741 			       else

  2742 				   (

  2743 				    Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2744 				    poly2term(sort (mv_geq LEX_) (x),vars) $

  2745 				    make_term(xs,vars)

  2746 				    );

  2748 (*. factorizes the denominator (polynomial representation) .*)

  2749 fun factorize_den (l,den,vars) =

  2750     let

  2751 	val factor_list=kill_equal( (get_factors l));

  2752 	val mlist=multi_list(factor_list);

  2753 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2754     in

  2755 	if rest=[] then

  2756 	    (

  2757 	     if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)

  2758 	     else make_term(last::factor_list,vars)

  2759 	     )

  2760 	else error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")

  2761     end;

  2763 (*. makes a term out of the elements of the list (expanded polynomial representation) .*)

  2764 fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2765   | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)

  2766 			       else

  2767 				   (

  2768 				    Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2769 				    poly2expanded(sort (mv_geq LEX_) (x),vars) $

  2770 				    make_exp(xs,vars)

  2771 				    );

  2773 (*. factorizes the denominator (expanded polynomial representation) .*)

  2774 fun factorize_den_exp (l,den,vars) =

  2775     let

  2776 	val factor_list=kill_equal( (get_factors l));

  2777 	val mlist=multi_list(factor_list);

  2778 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2779     in

  2780 	if rest=[] then

  2781 	    (

  2782 	     if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)

  2783 	     else make_exp(last::factor_list,vars)

  2784 	     )

  2785 	else error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")

  2786     end;

  2788 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2789 (*. the nominators and denominators with the correct factor .*)

  2790 (*. (polynomial representation) .*)

  2791 fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)

  2792   | step_add_list_of_fractions [x]= error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")

  2793   | step_add_list_of_fractions (xs) =

  2794     let

  2795         val den_list=termlist2denominators (xs); (* list of denominators *)

  2796 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2797 	val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2798     in

  2799 	com_den(xs,denom,den,var)

  2800     end;

  2802 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2803 (*. the nominators and denominators with the correct factor .*)

  2804 (*. (expanded polynomial representation) .*)

  2805 fun step_add_list_of_fractions_exp []  = (Free("0",HOLogic.realT),[]:term list)

  2806   | step_add_list_of_fractions_exp [x] = error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")

  2807   | step_add_list_of_fractions_exp (xs)=

  2808     let

  2809         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2810 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2811 	val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2812     in

  2813 	com_den_exp(xs,denom,den,var)

  2814     end;

  2816 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2817 -------------------------------------------------------------

  2818 (* WN0210???SK brauch ma des überhaupt *)

  2819 fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)

  2820   | step_add_list_of_fractions2 [x]=(x,[])

  2821   | step_add_list_of_fractions2 (xs) =

  2822     let

  2823         val den_list=termlist2denominators (xs); (* list of denominators *)

  2824 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2825 	val den=factorize_den(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2826     in

  2827 	(

  2828 	 Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2829 	 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $

  2830 	 poly2term(denom,var)

  2831 	,

  2832 	[]

  2833 	)

  2834     end;

  2836 (* WN0210???SK brauch ma des überhaupt *)

  2837 fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)

  2838   | step_add_list_of_fractions2_exp [x]=(x,[])

  2839   | step_add_list_of_fractions2_exp (xs) =

  2840     let

  2841         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2842 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2843 	val den=factorize_den_exp(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2844     in

  2845 	(

  2846 	 Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2847 	 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $

  2848 	 poly2expanded(denom,var)

  2849 	,

  2850 	[]

  2851 	)

  2852     end;

  2853 ---------------------------------------------- *)

  2856 (* converts a term, which contains several terms seperated by +, into a list of these terms .*)

  2857 fun term2list (t as (Const("Fields.inverse_class.divide",_) $ _ $ _)) = [t]

  2858   | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =

  2859       [Const ("Fields.inverse_class.divide",

  2860         [HOLogic.realT,HOLogic.realT] ---> HOLogic.realT) $

  2861 	  t $ Free("1",HOLogic.realT)

  2862      ]

  2863   | term2list (t as (Free(_,_))) =

  2864     [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2865 	  t $  Free("1",HOLogic.realT)

  2866      ]

  2867   | term2list (t as (Const("Groups.times_class.times",_) $ _ $ _)) =

  2868     [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2869 	  t $ Free("1",HOLogic.realT)

  2870      ]

  2871   | term2list (Const("Groups.plus_class.plus",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)

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

  2873     error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")

  2874   | term2list _ = error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");

  2876 (*.factors out the gcd of nominator and denominator:

  2877    a/b = (a' * gcd)/(b' * gcd),  a,b,gcd  are poly[2].*)

  2878 fun factout_p_  (thy:theory) t = SOME (step_cancel t,[]:term list);

  2880 (*.cancels a single fraction with normalform [2]

  2881    resulting in a canceled fraction [2], see factout_ .*)

  2882 fun cancel_p_ (thy:theory) t = (*WN.2.6.03 no rewrite -> NONE !*)

  2883     (let val (t',asm) = direct_cancel(*_expanded ... corrected MG.21.8.03*) t

  2884      in if t = t' then NONE else SOME (t',asm)

  2885      end) handle _ => NONE;

  2886 (*.the same as above with normalform [3]

  2887   val cancel_ :

  2888       theory ->        (*10.02 unused                                    *)

  2889       term -> 	       (*fraction in normalform [3]                      *)

  2890       (term * 	       (*fraction in normalform [3]                      *)

  2891        term list)      (*casual asumptions in normalform [3]             *)

  2892 	  option       (*NONE: the function is not applicable            *).*)

  2894 (*.transforms sums of at least 2 fractions [3] to

  2895    sums with the least common multiple as nominator.*)

  2896 fun common_nominator_p_ (thy:theory) t =

  2897 ((*tracing("### common_nominator_p_ called");*)

  2898     SOME (step_add_list_of_fractions(term2list(t))) handle _ => NONE

  2899 );

  2901 (*.add 2 or more fractions

  2902 val add_fraction_p_ :

  2903       theory ->        (*10.02 unused                                    *)

  2904       term -> 	       (*2 or more fractions with normalform [2]         *)

  2905       (term * 	       (*one fraction with normalform [2]                *)

  2906        term list)      (*casual assumptions in normalform [2] WN0210???SK  *)

  2907 	  option       (*NONE: the function is not applicable            *).*)

  2908 fun add_fraction_p_ (thy:theory) t =

  2909 (tracing("### add_fraction_p_ called");

  2910     (let val ts = term2list t

  2911      in if 1 < length ts

  2912 	then SOME (add_list_of_fractions ts)

  2913 	else NONE (*error ("RATIONALS_ADD_EXCEPTION: nothing to add")*)

  2914      end) handle _ => NONE

  2915 );

  2917 (*. brings the term into a normal form .*)

  2918 fun norm_rational_ (thy:theory) t =

  2919     SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;

  2920 fun norm_expanded_rat_ (thy:theory) t =

  2921     SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;

  2924 (*.evaluates conditions in calculate_Rational.*)

  2925 (*make local with FIXX@ME result:term *term list*)

  2926 val calc_rat_erls = prep_rls(

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

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

  2929 	 rules =

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

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

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

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

  2934 	  ],

  2935 	 scr = EmptyScr});

  2938 (*.simplifies expressions with numerals;

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

  2940    need to have constants to be commuted together respectively.*)

  2941 val calculate_Rational = prep_rls (merge_rls "calculate_Rational"

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

  2943 	    erls = calc_rat_erls, srls = Erls,

  2944 	    calc = [], errpatts = [],

  2945 	    rules =

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  2983 	       ], scr = EmptyScr})

  2984 	  calculate_Poly);

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

  2987 fun eval_is_expanded (thmid:string) _

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

  2989     if is_expanded arg

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

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

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

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

  2994   | eval_is_expanded _ _ _ _ = NONE;

  2996 val rational_erls =

  2997     merge_rls "rational_erls" calculate_Rational

  2998 	      (append_rls "is_expanded" Atools_erls

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

  3000 			   ]);

  3003 (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:

  3004  =================================================================

  3005  A[2] 'cancel_p': .

  3006  A[3] 'cancel': .

  3007  B[2] 'common_nominator_p': transforms summands in a term [2]

  3008          to fractions with the (least) common multiple as nominator.

  3009  B[3] 'norm_rational': normalizes arbitrary algebraic terms (without

  3010          radicals and transzendental functions) to one canceled fraction,

  3011 	 nominator and denominator in polynomial form.

  3013 In order to meet isac's requirements for interactive and stepwise calculation,

  3014 each 'reverse-rewerite-set' consists of an initialization for the interpreter

  3015 state and of 4 functions, each of which employs rewriting as much as possible.

  3016 The signature of these functions are the same in each 'reverse-rewrite-set'

  3017 respectively.*)

  3019 (* ************************************************************************* *)

  3021 local(*. cancel_p

  3022 ------------------------

  3023 cancels a single fraction consisting of two (uni- or multivariate)

  3024 polynomials WN0609???SK[2] into another such a fraction; examples:

  3026 	   a^2 + -1*b^2         a + b

  3027         -------------------- = ---------

  3028 	a^2 + -2*a*b + b^2     a + -1*b

  3030         a^2    a

  3031         --- = ---

  3032          a     1

  3034 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)

  3035 (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)

  3037 val {rules, rew_ord=(_,ro),...} =

  3038     rep_rls (assoc_rls "make_polynomial");

  3039 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3040            see rational.sml --- investigate rulesets for cancel_p ---*)

  3041 val {rules, rew_ord=(_,ro),...} =

  3042     rep_rls (assoc_rls "rev_rew_p");

  3044 (*.init_state = fn : term -> istate

  3045 initialzies the state of the script interpreter. The state is:

  3047 type rrlsstate =      (*state for reverse rewriting*)

  3048      (term *          (*the current formula*)

  3049       term *          (*the final term*)

  3050       rule list       (*'reverse rule list' (#)*)

  3051 	    list *    (*may be serveral, eg. in norm_rational*)

  3052       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3053        (term *        (*... rewrite with ...*)

  3054 	term list))   (*... assumptions*)

  3055 	  list);      (*derivation from given term to normalform

  3056 		       in reverse order with sym_thm;

  3057                        (#) could be extracted from here by (map #1)*).*)

  3058 (* val {rules, rew_ord=(_,ro),...} =

  3059        rep_rls (assoc_rls "rev_rew_p")        (*USE ALWAYS, SEE val cancel_p*);

  3060    val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*);

  3061    val t = t;

  3062    *)

  3063 fun init_state thy eval_rls ro t =

  3064     let val SOME (t',_) = factout_p_ thy t

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

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

  3067         val der = der @ [(Thm ("real_mult_div_cancel2",

  3068 			       num_str @{thm real_mult_div_cancel2}),

  3069 			  (t'',asm))]

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

  3071 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3072 				      "sym_real_mult_0",

  3073 				      "sym_real_mult_1"

  3074 				      (*..insufficient,eg.make_Polynomial*)])rs

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

  3077 (*.locate_rule = fn : rule list -> term -> rule

  3078 		      -> (rule * (term * term list) option) list.

  3079   checks a rule R for being a cancel-rule, and if it is,

  3080   then return the list of rules (+ the terms they are rewriting to)

  3081   which need to be applied before R should be applied.

  3082   precondition: the rule is applicable to the argument-term.

  3083 arguments:

  3084   rule list: the reverse rule list

  3085   -> term  : ... to which the rule shall be applied

  3086   -> rule  : ... to be applied to term

  3087 value:

  3088   -> (rule           : a rule rewriting to ...

  3089       * (term        : ... the resulting term ...

  3090          * term list): ... with the assumptions ( //#0).

  3091       ) list         : there may be several such rules;

  3092 		       the list is empty, if the rule has nothing to do

  3093 		       with cancelation.*)

  3094 (* val () = ();

  3095    *)

  3096 fun locate_rule thy eval_rls ro [rs] t r =

  3097     if (id_of_thm r) mem (map (id_of_thm)) rs

  3098     then let val ropt =

  3099 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3100 	 in case ropt of

  3101 		SOME ta => [(r, ta)]

  3102 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3103 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3104 			 []) end

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

  3106   | locate_rule _ _ _ _ _ _ =

  3107     error ("locate_rule: doesnt match rev-sets in istate");

  3109 (*.next_rule = fn : rule list -> term -> rule option

  3110   for a given term return the next rules to be done for cancelling.

  3111 arguments:

  3112   rule list     : the reverse rule list

  3113   term          : the term for which ...

  3114 value:

  3115   -> rule option: ... this rule is appropriate for cancellation;

  3116 		  there may be no such rule (if the term is canceled already.*)

  3117 (* val thy = thy;

  3118    val Rrls {rew_ord=(_,ro),...} = cancel;

  3119    val ([rs],t) = (rss,f);

  3120    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3122    val (thy, [rs]) = (thy, revsets);

  3123    val Rrls {rew_ord=(_,ro),...} = cancel;

  3124    nex [rs] t;

  3125    *)

  3126 fun next_rule thy eval_rls ro [rs] t =

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

  3128     in case der of

  3129 (* val (_,r,_)::_ = der;

  3130    *)

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

  3132 	 | _ => NONE

  3133     end

  3134   | next_rule _ _ _ _ _ =

  3135     error ("next_rule: doesnt match rev-sets in istate");

  3137 (*.val attach_form = f : rule list -> term -> term

  3138 			 -> (rule * (term * term list)) list

  3139   checks an input term TI, if it may belong to a current cancellation, by

  3140   trying to derive it from the given term TG.

  3141 arguments:

  3142   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3143   -> term: TI, the next one input by the user

  3144 value:

  3145   -> (rule           : the rule to be applied in order to reach TI

  3146       * (term        : ... obtained by applying the rule ...

  3147          * term list): ... and the respective assumptions.

  3148       ) list         : there may be several such rules;

  3149                        the list is empty, if the users term does not belong

  3150 		       to a cancellation of the term last agreed upon.*)

  3151 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

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

  3154 in

  3156 val cancel_p =

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

  3158 	  rew_ord=("ord_make_polynomial",

  3159 		   ord_make_polynomial false thy),

  3160 	  erls = rational_erls,

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

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

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

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

  3165 	  errpatts = [],

  3166 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3167 		     normal_form = cancel_p_ thy,

  3168 		     locate_rule = locate_rule thy Atools_erls ro,

  3169 		     next_rule   = next_rule thy Atools_erls ro,

  3170 		     attach_form = attach_form}}

  3171 end;(*local*)

  3173 local(*.ad [2] 'common_nominator_p'

  3174 ---------------------------------

  3175 FIXME Beschreibung .*)

  3178 val {rules=rules,rew_ord=(_,ro),...} =

  3179     rep_rls (assoc_rls "make_polynomial");

  3180 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3181            see rational.sml --- investigate rulesets for cancel_p ---*)

  3182 val {rules, rew_ord=(_,ro),...} =

  3183     rep_rls (assoc_rls "rev_rew_p");

  3184 val thy = thy;

  3187 (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option

  3188   as defined above*)

  3190 (*.init_state = fn : term -> istate

  3191 initialzies the state of the interactive interpreter. The state is:

  3193 type rrlsstate =      (*state for reverse rewriting*)

  3194      (term *          (*the current formula*)

  3195       term *          (*the final term*)

  3196       rule list       (*'reverse rule list' (#)*)

  3197 	    list *    (*may be serveral, eg. in norm_rational*)

  3198       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3199        (term *        (*... rewrite with ...*)

  3200 	term list))   (*... assumptions*)

  3201 	  list);      (*derivation from given term to normalform

  3202 		       in reverse order with sym_thm;

  3203                        (#) could be extracted from here by (map #1)*).*)

  3204 fun init_state thy eval_rls ro t =

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

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

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

  3208         val der = der @ [(Thm ("real_mult_div_cancel2",

  3209 			       num_str @{thm real_mult_div_cancel2}),

  3210 			  (t'',asm))]

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

  3212 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3213 				      "sym_real_mult_0",

  3214 				      "sym_real_mult_1"]) rs;

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

  3217 (* use"knowledge/Rational.ML";

  3218    *)

  3220 (*.locate_rule = fn : rule list -> term -> rule

  3221 		      -> (rule * (term * term list) option) list.

  3222   checks a rule R for being a cancel-rule, and if it is,

  3223   then return the list of rules (+ the terms they are rewriting to)

  3224   which need to be applied before R should be applied.

  3225   precondition: the rule is applicable to the argument-term.

  3226 arguments:

  3227   rule list: the reverse rule list

  3228   -> term  : ... to which the rule shall be applied

  3229   -> rule  : ... to be applied to term

  3230 value:

  3231   -> (rule           : a rule rewriting to ...

  3232       * (term        : ... the resulting term ...

  3233          * term list): ... with the assumptions ( //#0).

  3234       ) list         : there may be several such rules;

  3235 		       the list is empty, if the rule has nothing to do

  3236 		       with cancelation.*)

  3237 (* val () = ();

  3238    *)

  3239 fun locate_rule thy eval_rls ro [rs] t r =

  3240     if (id_of_thm r) mem (map (id_of_thm)) rs

  3241     then let val ropt =

  3242 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3243 	 in case ropt of

  3244 		SOME ta => [(r, ta)]

  3245 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3246 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3247 			 []) end

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

  3249   | locate_rule _ _ _ _ _ _ =

  3250     error ("locate_rule: doesnt match rev-sets in istate");

  3252 (*.next_rule = fn : rule list -> term -> rule option

  3253   for a given term return the next rules to be done for cancelling.

  3254 arguments:

  3255   rule list     : the reverse rule list

  3256   term          : the term for which ...

  3257 value:

  3258   -> rule option: ... this rule is appropriate for cancellation;

  3259 		  there may be no such rule (if the term is canceled already.*)

  3260 (* val thy = thy;

  3261    val Rrls {rew_ord=(_,ro),...} = cancel;

  3262    val ([rs],t) = (rss,f);

  3263    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3265    val (thy, [rs]) = (thy, revsets);

  3266    val Rrls {rew_ord=(_,ro),...} = cancel;

  3267    nex [rs] t;

  3268    *)

  3269 fun next_rule thy eval_rls ro [rs] t =

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

  3271     in case der of

  3272 (* val (_,r,_)::_ = der;

  3273    *)

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

  3275 	 | _ => NONE

  3276     end

  3277   | next_rule _ _ _ _ _ =

  3278     error ("next_rule: doesnt match rev-sets in istate");

  3280 (*.val attach_form = f : rule list -> term -> term

  3281 			 -> (rule * (term * term list)) list

  3282   checks an input term TI, if it may belong to a current cancellation, by

  3283   trying to derive it from the given term TG.

  3284 arguments:

  3285   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3286   -> term: TI, the next one input by the user

  3287 value:

  3288   -> (rule           : the rule to be applied in order to reach TI

  3289       * (term        : ... obtained by applying the rule ...

  3290          * term list): ... and the respective assumptions.

  3291       ) list         : there may be several such rules;

  3292                        the list is empty, if the users term does not belong

  3293 		       to a cancellation of the term last agreed upon.*)

  3294 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

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

  3297 (* if each pat* matches with the current term, the Rrls is applied

  3298    (there are no preconditions to be checked, they are @{term True}) *)

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

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

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

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

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

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

  3305 in

  3307 (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;

  3308   besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht

  3309   dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)

  3310 val common_nominator_p =

  3311     Rrls {id = "common_nominator_p", prepat=prepat,

  3312 	  rew_ord=("ord_make_polynomial",

  3313 		   ord_make_polynomial false thy),

  3314 	  erls = rational_erls,

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

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

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

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

  3319 	  errpatts = [],

  3320 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3321 		     normal_form = add_fraction_p_ thy,(*FIXME.WN0211*)

  3322 		     locate_rule = locate_rule thy Atools_erls ro,

  3323 		     next_rule   = next_rule thy Atools_erls ro,

  3324 		     attach_form = attach_form}}

  3325 end;(*local*)

  3327 end; (*struct*)

  3329 *}

  3330 ML {*

  3331 open RationalI;

  3332 (*##*)

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

  3335 fun is_ratpolyexp (Free _) = true

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

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

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

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

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

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

  3342                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

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

  3344                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

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

  3346                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

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

  3348                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

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

  3350                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3351   | is_ratpolyexp _ = false;

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

  3354 fun eval_is_ratpolyexp (thmid:string) _

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

  3356     if is_ratpolyexp arg

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

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

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

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

  3361   | eval_is_ratpolyexp _ _ _ _ = NONE;

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

  3364 fun eval_get_denominator (thmid:string) _

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

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

  3367                 denom)) thy =

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

  3369 	            Trueprop $ (mk_equality (t, denom)))

  3370   | eval_get_denominator _ _ _ _ = NONE;

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

  3373 fun eval_get_numerator (thmid:string) _

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

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

  3376             $denom )) thy =

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

  3378 	    Trueprop $ (mk_equality (t, num)))

  3379   | eval_get_numerator _ _ _ _ = NONE;

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

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

  3384 (*WN100906 removed in favour of discard_minus = discard_minus_ formerly

  3385    discard binary minus, shift unary minus into -1*;

  3386    unary minus before numerals are put into the numeral by parsing;

  3387    contains absolute minimum of thms for context in norm_Rational

  3388 *** val discard_minus = prep_rls(

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

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

  3391       rules =

  3392         [Thm ("real_diff_minus", num_str @{thm real_diff_minus}),

  3393 	           (*"a - b = a + -1 * b"*)

  3394 	         Thm ("sym_real_mult_minus1", num_str (@{thm real_mult_minus1} RS @{thm sym}))

  3395 	           (*- ?z = "-1 * ?z"*)],

  3396       scr = EmptyScr

  3397       }):rls;

  3398  *)

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

  3401 val powers_erls = prep_rls(

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

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

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

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

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

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

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

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

  3410 	       ],

  3411       scr = EmptyScr

  3412       }:rls);

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

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

  3415 val powers = prep_rls(

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  3443 	       ],

  3444       scr = EmptyScr

  3445       }:rls);

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

  3447 val rat_mult_divide = prep_rls(

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

  3449        rew_ord = ("dummy_ord",dummy_ord),

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

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

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

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

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

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

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

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

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

  3460 	       Thm ("divide_divide_eq_right",

  3461                      num_str @{thm divide_divide_eq_right}),

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

  3463 	       Thm ("divide_divide_eq_left",

  3464                      num_str @{thm divide_divide_eq_left}),

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

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

  3467 	       ],

  3468       scr = EmptyScr

  3469       }:rls);

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

  3472 val reduce_0_1_2 = prep_rls(

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

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

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

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

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

  3478 	       (*"1 * z = z"*)

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

  3480 	       "-1 * z = - z"*)

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

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

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

  3485 	       (*"0 * z = 0"*)

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

  3487 	       (*"0 + z = z"*)

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

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

  3491 	       Thm ("sym_real_mult_2",

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

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

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

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

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

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

  3499 	       ], scr = EmptyScr}:rls);

  3501 (*erls for calculate_Rational;

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

  3503 val norm_rat_erls = prep_rls(

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

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

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

  3507 	       ], scr = EmptyScr}:rls);

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

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

  3511 which is now replaced by MGs version below

  3512 vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)

  3513 val norm_Rational = prep_rls(

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

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

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

  3517 	       Rls_ discard_minus,

  3518 	       Rls_ powers,

  3519 	       Rls_ rat_mult_divide,

  3520 	       Rls_ expand,

  3521 	       Rls_ reduce_0_1_2,

  3522 	       (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)

  3523 	       Rls_ order_add_mult,

  3524 	       Rls_ collect_numerals,

  3525 	       Rls_ add_fractions_p,

  3526 	       Rls_ cancel_p

  3527 	       ],

  3528       scr = EmptyScr}:rls);

  3530 val norm_Rational_parenthesized = prep_rls(

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

  3532        rew_ord = ("dummy_ord", dummy_ord),

  3533       erls = Atools_erls, srls = Erls,

  3534       calc = [], errpatts = [],

  3535       rules = [Rls_  norm_Rational, (*from RL -- not the latest one*)

  3536 	       Rls_ discard_parentheses

  3537 	       ],

  3538       scr = EmptyScr}:rls);

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

  3541 val simplify_rational =

  3542     merge_rls "simplify_rational" expand_binoms

  3543     (append_rls "divide" calculate_Rational

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  3559 (*

  3560 ,

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

  3562 *)

  3563 		 ]);

  3564 *}

  3565 ML {*

  3567 (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)

  3569 (* ------------------------------------------------------------------ *)

  3570 (*                  Simplifier für beliebige Buchterme                *)

  3571 (* ------------------------------------------------------------------ *)

  3572 (*----------------------- norm_Rational_mg ---------------------------*)

  3573 (*. description of the simplifier see MG-DA.p.56ff .*)

  3574 (* ------------------------------------------------------------------- *)

  3575 val common_nominator_p_rls = prep_rls(

  3576   Rls {id = "common_nominator_p_rls", preconds = [],

  3577        rew_ord = ("dummy_ord",dummy_ord),

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

  3579 	 rules =

  3580 	 [Rls_ common_nominator_p

  3581 	  (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?

  3582 	    FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)

  3583 	  ],

  3584 	 scr = EmptyScr});

  3585 (* ------------------------------------------------------------------- *)

  3586 val cancel_p_rls = prep_rls(

  3587   Rls {id = "cancel_p_rls", preconds = [],

  3588        rew_ord = ("dummy_ord",dummy_ord),

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

  3590 	 rules =

  3591 	 [Rls_ cancel_p

  3592 	  (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)

  3593 	  ],

  3594 	 scr = EmptyScr});

  3595 (* -------------------------------------------------------------------- *)

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

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

  3598 val rat_mult_poly = prep_rls(

  3599   Rls {id = "rat_mult_poly", preconds = [],

  3600        rew_ord = ("dummy_ord",dummy_ord),

  3601 	 erls =  append_rls "e_rls-is_polyexp" e_rls

  3602 	         [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

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

  3604 	 rules =

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

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

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

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

  3609 	  ],

  3610 	 scr = EmptyScr});

  3612 (* ------------------------------------------------------------------ *)

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

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

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

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

  3617     ... WN0609???MG.*)

  3618 val rat_mult_div_pow = prep_rls(

  3619   Rls {id = "rat_mult_div_pow", preconds = [],

  3620        rew_ord = ("dummy_ord",dummy_ord),

  3621        erls = e_rls,

  3622        (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls

  3623 			[Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3624          with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get

  3625 	 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.

  3626          thus we decided to go on with this flaw*)

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

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

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

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

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

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

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

  3635 	       Thm ("real_divide_divide1_mg",

  3636                      num_str @{thm real_divide_divide1_mg}),

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

  3638 	       Thm ("divide_divide_eq_right",

  3639                      num_str @{thm divide_divide_eq_right}),

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

  3641 	       Thm ("divide_divide_eq_left",

  3642                      num_str @{thm divide_divide_eq_left}),

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

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

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

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

  3648 	       ],

  3649       scr = EmptyScr}:rls);

  3650 (* ------------------------------------------------------------------ *)

  3651 val rat_reduce_1 = prep_rls(

  3652   Rls {id = "rat_reduce_1", preconds = [],

  3653        rew_ord = ("dummy_ord",dummy_ord),

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

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

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

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

  3658 		(*"1 * z = z"*)

  3659 		],

  3660        scr = EmptyScr}:rls);

  3661 (* ------------------------------------------------------------------ *)

  3662 (*. looping part of norm_Rational(*_mg*) .*)

  3663 val norm_Rational_rls = prep_rls(

  3664    Rls {id = "norm_Rational_rls", preconds = [],

  3665        rew_ord = ("dummy_ord",dummy_ord),

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

  3667        rules = [Rls_ common_nominator_p_rls,

  3668 		Rls_ rat_mult_div_pow,

  3669 		Rls_ make_rat_poly_with_parentheses,

  3670 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)

  3671 		Rls_ rat_reduce_1

  3672 		],

  3673        scr = EmptyScr}:rls);

  3674 (* ------------------------------------------------------------------ *)

  3675 (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)

  3676  just be renaming:*)

  3677 val norm_Rational(*_mg*) = prep_rls(

  3678    Seq {id = "norm_Rational"(*_mg*), preconds = [],

  3679        rew_ord = ("dummy_ord",dummy_ord),

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

  3681        rules = [Rls_ discard_minus,

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

  3683 		Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)

  3684 		Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)

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

  3686 		Rls_ discard_parentheses1 (* mult only                       *)

  3687 		],

  3688        scr = EmptyScr}:rls);

  3689 (* ------------------------------------------------------------------ *)

  3691 *}

  3692 ML {*

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

  3694   [("calculate_Rational", calculate_Rational),

  3695    ("calc_rat_erls",calc_rat_erls),

  3696    ("rational_erls", rational_erls),

  3697    ("cancel_p", cancel_p),

  3698    ("common_nominator_p", common_nominator_p),

  3699    ("common_nominator_p_rls", common_nominator_p_rls),

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

  3701    ("powers_erls", powers_erls),

  3702    ("powers", powers),

  3703    ("rat_mult_divide", rat_mult_divide),

  3704    ("reduce_0_1_2", reduce_0_1_2),

  3705    ("rat_reduce_1", rat_reduce_1),

  3706    ("norm_rat_erls", norm_rat_erls),

  3707    ("norm_Rational", norm_Rational),

  3708    ("norm_Rational_rls", norm_Rational_rls),

  3709    ("norm_Rational_parenthesized", norm_Rational_parenthesized),

  3710    ("rat_mult_poly", rat_mult_poly),

  3711    ("rat_mult_div_pow", rat_mult_div_pow),

  3712    ("cancel_p_rls", cancel_p_rls)

  3713    ]);

  3715 calclist':= overwritel (!calclist',

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

  3717     ]);

  3719 (** problems **)

  3721 store_pbt

  3722  (prep_pbt thy "pbl_simp_rat" [] e_pblID

  3723  (["rational","simplification"],

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

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

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

  3727   ],

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

  3729   SOME "Simplify t_t",

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

  3732 (** methods **)

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

  3735   of norm_Rational in order to ease repair on inform*)

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

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

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

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

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

  3741 	   prls = append_rls "simplification_of_rationals_prls" e_rls

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

  3743 				   crls = e_rls, errpats = [],

  3744    nrls = norm_Rational_rls},

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

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

  3747    "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^

  3748    "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^

  3749    "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^

  3750    "    (Repeat                                                     " ^

  3751    "     ((Try (Rewrite_Set common_nominator_p_rls False) @@        " ^

  3752    "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^

  3753    "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^

  3754    "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^

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

  3756    "    Try (Rewrite_Set discard_parentheses1 False))               " ^

  3757    "   t_t)"));

  3759 *}

  3760 ML {*"test"*}

  3761 end (* theory*)

author	Walther Neuper <neuper@ist.tugraz.at>
	Mon, 26 Aug 2013 15:37:48 +0200
changeset 52091	fb5a2c25e3f0
parent 52089	7a740eedefc0
child 52092	fbd18c58a894
permissions	-rwxr-xr-x