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

     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 mk_asms baseT ts =

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

   295   in map (mk_noteq_0 baseT) as' end

   297 fun check_fraction t =

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

   299   in SOME (numerator, denominator) end

   300   handle Bind => NONE

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

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

   304 fun factout_p_ (thy: theory) t =

   305   let val opt = check_fraction t

   306   in

   307     case opt of

   308       NONE => NONE

   309     | SOME (numerator, denominator) =>

   310       let

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

   312         val baseT = type_of numerator

   313         val expT = HOLogic.realT

   314       in

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

   316           (SOME a, SOME b) =>

   317             let

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

   319               val es = replicate (length vs) 0

   320             in

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

   322               then NONE

   323               else

   324                 let

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

   326                   val ct = term_of_poly baseT expT vs c

   327                   val t' =

   328                     HOLogic.mk_binop "Fields.inverse_class.divide"

   329                       (HOLogic.mk_binop "Groups.times_class.times"

   330                         (term_of_poly baseT expT vs a', ct),

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

   332                   val asm = mk_asms baseT [b't, ct]

   333                 in SOME (t', asm) end

   334             end

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

   336       end

   337   end

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

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

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

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

   343   yields

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

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

   346     \<or> NONE *)

   347 fun cancel_p_ (thy: theory) t =

   348   let val opt = check_fraction t

   349   in

   350     case opt of

   351       NONE => NONE

   352     | SOME (numerator, denominator) =>

   353       let

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

   355         val baseT = type_of numerator

   356         val expT = HOLogic.realT

   357       in

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

   359           (SOME a, SOME b) =>

   360             let

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

   362               val es = replicate (length vs) 0

   363             in

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

   365               then NONE

   366               else

   367                 let

   368                   val bt' = term_of_poly baseT expT vs b'

   369                   val ct = term_of_poly baseT expT vs c

   370                   val t' =

   371                     HOLogic.mk_binop "Fields.inverse_class.divide"

   372                       (term_of_poly baseT expT vs a', bt')

   373                   val asm = mk_asms baseT [bt']

   374                 in SOME (t', asm) end

   375             end

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

   377       end

   378   end

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

   381 fun norm_frac_sum

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

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

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

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

   386   | norm_frac_sum

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

   388       nofrac $

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

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

   391   | norm_frac_sum

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

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

   394       nofrac)

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

   396   | norm_frac_sum _ = NONE

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

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

   400 fun common_nominator_p_ (thy: theory) t =

   401   let val opt = norm_frac_sum t

   402   in

   403     case opt of

   404       NONE => NONE

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

   406       let

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

   408         val baseT = type_of n1

   409         val expT = HOLogic.realT

   410       in

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

   412           (SOME a, SOME b) =>

   413             let

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

   415               val d1' = term_of_poly baseT expT vs a'

   416               val d2' = term_of_poly baseT expT vs b'

   417               val c' = term_of_poly baseT expT vs c

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

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

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

   421                 --------------------------------------------------------------------------*)

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

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

   424               (*--------------------------------------------------------------------------*)

   425               val t' =

   426                 HOLogic.mk_binop "Groups.plus_class.plus"

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

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

   429                   HOLogic.mk_binop "Fields.inverse_class.divide"

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

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

   432             in SOME (t', asm) end

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

   434       end

   435   end

   437 (* add fractions

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

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

   440 fun add_fraction_p_ (thy: theory) t =

   441   let val opt = norm_frac_sum t

   442   in

   443     case opt of

   444       NONE => NONE

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

   446       let

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

   448         val baseT = type_of n1

   449         val expT = HOLogic.realT

   450       in

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

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

   453             let

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

   455               val nomin = term_of_poly baseT expT vs

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

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

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

   459               val asm = mk_asms baseT [denom]

   460             in SOME (t', asm) end

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

   462       end

   463   end

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

   466 fun xs union [] = xs

   467   | [] union ys = ys

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

   470 (*. gcd of integers .*)

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

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

   473 		  else gcd_int b (a mod b);

   475 (*. univariate polynomials (uv) .*)

   476 (*. univariate polynomials are represented as a list

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

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

   479 type uv_poly = int list;

   481 (*. adds two uv polynomials .*)

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

   483   | uv_mod_add_poly (p1,[]) = p1

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

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

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

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

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

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

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

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

   495 fun uv_mod_deg ([]:uv_poly) = 0

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

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

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

   500 fun uv_mod_mod2(x,p)=

   501     let

   502 	val y=(x mod p);

   503     in

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

   505 	    (

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

   507 	     )

   508     end;

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

   511 fun uv_mod_list_modp [] p = []

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

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

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

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

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

   519 (*. only for uv_mod_divides .*)

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

   521    integer zu klein  *)

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

   523     error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

   524   | uv_mod_pdiv p1 [x] =

   525     let

   526 	val xs= Unsynchronized.ref  [];

   527     in

   528 	if x<>0 then

   529 	    (

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

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

   532 	     )

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

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

   535     end

   536   | uv_mod_pdiv p1 p2 =

   537     let

   538 	val n= uv_mod_deg(p2);

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

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

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

   542 	val lc2=hd(p2');

   543 	val q= Unsynchronized.ref  [];

   544 	val c= Unsynchronized.ref  0;

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

   546     in

   547 	(

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

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

   550 	 else

   551 	     (

   552 	      if (!m)<n then

   553 		  (

   554 		   output:=([0],p1)

   555 		   )

   556 	      else

   557 		  (

   558 		   while (!m)>=n do

   559 		       (

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

   561 			if !c<>0 then

   562 			    (

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

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

   565 			     m:=uv_mod_deg(!p1')

   566 			     )

   567 			else m:=0

   568 			);

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

   570 		   )

   571 	      );

   572 	     !output

   573 	 )

   574     end;

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

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

   578     let

   579 	val n=uv_mod_deg(p2);

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

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

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

   583 	val lc2=hd(p2');

   584 	val q= Unsynchronized.ref  [];

   585 	val c= Unsynchronized.ref  0;

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

   587     in

   588 	(

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

   590 	 else

   591 	     (

   592 	      if (!m)<n then

   593 		  (

   594 		   output:=([0],p1)

   595 		   )

   596 	      else

   597 		  (

   598 		   while (!m)>=n do

   599 		       (

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

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

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

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

   604 			m:=(!m)-1

   605 			);

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

   608 		       (

   609 			p1':=tl(!p1')

   610 			);

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

   613 		   )

   614 	      );

   615 	     !output:uv_poly * uv_poly

   616 	 )

   617     end;

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

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

   621   | uv_mod_prest [] p2 = []:uv_poly

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

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

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

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

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

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

   630 fun uv_mod_cont ([]:uv_poly) = 0

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

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

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

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

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

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

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

   640   | uv_mod_pp p =

   641     let

   642 	val c= Unsynchronized.ref  0;

   643     in

   644 	(

   645 	 c:=uv_mod_cont(p);

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

   648 	 else uv_mod_div_list(p,!c)

   649 	)

   650     end;

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

   653 fun uv_mod_lc ([]:uv_poly) = 0

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

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

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

   658     let

   659 	val f = Unsynchronized.ref  [];

   660 	val f'= Unsynchronized.ref  p2;

   661 	val fi= Unsynchronized.ref  [];

   662     in

   663 	(

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

   665  	 while uv_mod_deg(!f')>0 do

   666 	     (

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

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

   669 		  (

   670 		   fi:=(!f');

   671 		   f:=(!fi)::(!f)

   672 		   )

   673 	      else ()

   674 	      );

   675 	      (!f)

   677 	 )

   678     end;

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

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

   682   | uv_mod_gcd_modp p1 [] p= p1

   683   | uv_mod_gcd_modp p1 p2 p=

   684     let

   685 	val p1'= Unsynchronized.ref [];

   686 	val p2'= Unsynchronized.ref [];

   687 	val pc= Unsynchronized.ref [];

   688 	val g= Unsynchronized.ref  [];

   689 	val d= Unsynchronized.ref  0;

   690 	val prs= Unsynchronized.ref  [];

   691     in

   692 	(

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

   694 	     (

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

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

   697 	      )

   698 	 else

   699 	     (

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

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

   702 	      );

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

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

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

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

   709 	 else pc:=hd(!prs);

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

   712 	 !g

   713 	 )

   714     end;

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

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

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

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

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

   723 fun uv_mod_add_qu [] = 0.0

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

   726 (*. calculates the euklidean norm .*)

   727 fun uv_mod_norm ([]:uv_poly) = 0.0

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

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

   731 fun uv_mod_multi a b = a * b;

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

   734 fun uv_mod_prim(x,[])= false

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

   736 		else false

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

   738 			then

   739 			    if uv_mod_prim(x,ys) then true

   740 			    else false

   741 		    else false;

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

   744 fun uv_mod_nextprime(g,p)=

   745     let

   746 	val list= Unsynchronized.ref  [2];

   747 	val exit= Unsynchronized.ref  0;

   748 	val i = Unsynchronized.ref 2

   749     in

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

   751 	    (

   752 	     if uv_mod_prim(!i,!list) then

   753 		 (

   754 		  if (p mod !i <> 0)

   755 		      then

   756 			  (

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

   758 			   i:= (!i)+1

   759 			   )

   760 		  else i:=(!i)+1

   761 		  )

   762 	     else i:= (!i)+1

   763 		 );

   764 	    i:=(p+1);

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

   766 	    (

   767 	     if uv_mod_prim(!i,!list) then

   768 		 (

   769 		  if (g mod !i <> 0)

   770 		      then

   771 			  (

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

   773 			   exit:= (!i)

   774 			   )

   775 		  else i:=(!i)+1

   776 		      )

   777 	     else i:= (!i)+1

   778 		 );

   779 	    !exit

   780     end;

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

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

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

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

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

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

   790 (*. chinese remainder algorithm .*)

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

   792     let

   793 	val c= Unsynchronized.ref  0;

   794 	val r1'= Unsynchronized.ref  0;

   795 	val d= Unsynchronized.ref  0;

   796 	val a= Unsynchronized.ref  0;

   797     in

   798 	(

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

   800 	     (

   801 	      c:=(!c)+1

   802 	      );

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

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

   805 	 !r1'+(!d)*m1

   806 	 )

   807     end;

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

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

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

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

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

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

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

   817     let

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

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

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

   821 	val temp= Unsynchronized.ref  [];

   822 	val cp= Unsynchronized.ref  [];

   823 	val qp= Unsynchronized.ref  [];

   824 	val q= Unsynchronized.ref [];

   825 	val pn= Unsynchronized.ref  0;

   826 	val d= Unsynchronized.ref  0;

   827 	val g1= Unsynchronized.ref  0;

   828 	val p= Unsynchronized.ref  0;

   829 	val m= Unsynchronized.ref  0;

   830 	val exit= Unsynchronized.ref  0;

   831 	val i= Unsynchronized.ref  1;

   832     in

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

   834 	else

   835 	    (

   836 	     temp:= !p1;

   837 	     p1:= !p2;

   838 	     p2:= !temp

   839 	     );

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

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

   844 	p:=4;

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

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

   848 	  Real.fromInt(!d) *

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

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

   852 	while (!exit=0) do

   853 	    (

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

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

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

   857 		 (

   858 		  i:=1;

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

   860 		      (

   861 		       i:=(!i)+1

   862 		       );

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

   864 		  )

   865 	     else ();

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

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

   871 	     pn:=(!p);

   872 	     q:=(!qp);

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

   875 		 (

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

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

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

   879  		      (

   880  		       i:=1;

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

   882  			   (

   883  			    i:=(!i)+1

   884 		           );

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

   886  		      )

   887  		  else ();

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

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

   891  		      (

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

   893  		       )

   894  		  else

   895  		      (

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

   897  			   (

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

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

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

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

   902 		 	    )

   903 		       else

   904 			   (

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

   906 				(

   907 				 pn:= !p;

   908 				 q:= !qp

   909 				 )

   910 			    else ()

   911 			    )

   912 		       )

   913 		  );

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

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

   916 	     );

   917 	    uv_mod_smul_poly(!q,c):uv_poly

   918     end;

   920 (*. multivariate polynomials .*)

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

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

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

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

   926 (*. global variables .*)

   927 (*. order indicators .*)

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

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

   931 (*. datatypes for internal representation.*)

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

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

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

   936 type mv_poly = mv_monom list;

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

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

   940 fun mv_mg_hlp([]) = EQUAL

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

   942 		    else if x>0 then GREATER

   943 			 else mv_mg_hlp(list);

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

   946 fun mv_addlist([]) = 0

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

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

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

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

   952     if order=LEX_ then

   953 	(

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

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

   956 	     )

   957     else

   958 	if order=GGO_ then

   959 	    (

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

   961 	     else

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

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

   964 	     )

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

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

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

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

   970 let

   971     val temp= Unsynchronized.ref  EQUAL;

   972 in

   973     if order=LEX_ then

   974 	(

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

   976 	     error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   977 	 else

   978 	     (

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

   980 	      if !temp=EQUAL then

   981 		  ( if x1=x2 then EQUAL

   982 		    else if x1>x2 then GREATER

   983 			 else LESS

   984 			     )

   985 	      else (!temp)

   986 	      )

   987 	     )

   988     else

   989 	if order=GGO_ then

   990 	    (

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

   992 	      error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   993 	     else

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

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

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

   997 		     )

   998 	else error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")

   999 end;

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

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

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

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

  1006 (*. leading power product .*)

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

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

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

  1011 (*. leading monomial .*)

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

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

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

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

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

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

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

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

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

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

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

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

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

  1029   | mv_lc(p1,order)  =

  1030     let

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

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

  1033 	val lc= Unsynchronized.ref  [];

  1034     in

  1035 	(

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

  1037 	     (

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

  1039 	      p1o:=tl(!p1o)

  1040 	      );

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

  1042 	 mv_rev_to(!lc)

  1043 	 )

  1044     end;

  1046 (*. compares two powerproducts .*)

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

  1049 (*. help function for mv_add .*)

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

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

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

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

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

  1055   | mv_madd(p1,p2,order) =

  1056     (

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

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

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

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

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

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

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

  1064 	      )

  1066 (*. adds two multivariate polynomials .*)

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

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

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

  1071 (*. monom multiplication .*)

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

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

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

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

  1078 (*. zeros a list .*)

  1079 fun mv_null2([])=[]

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

  1082 (*. multiplies two multivariate polynomials .*)

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

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

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

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

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

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

  1090 fun mv_getmax([])=0

  1091   | mv_getmax(x::p1)= let

  1092 		       val m=mv_getmax(p1);

  1093 		   in

  1094 		       if m>x then m

  1095 		       else x

  1096 		   end;

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

  1098 fun mv_grad([]:mv_poly) = 0

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

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

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

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

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

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

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

  1109 fun mv_is_negativ([])=false

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

  1112 (*. division of monomials .*)

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

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

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

  1116     let

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

  1118 	val pp= Unsynchronized.ref  [];

  1119     in

  1120 	(

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

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

  1123 	     if #1(p2)<>0 then

  1124 		 (

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

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

  1127 		  else (!c,!pp)

  1128 		      )

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

  1130 		 )

  1131     end;

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

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

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

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

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

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

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

  1142   | mv_division(f,g,order)=

  1143     let

  1144 	val r= Unsynchronized.ref  [];

  1145 	val q= Unsynchronized.ref  [];

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

  1147 	val k= Unsynchronized.ref  0;

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

  1149 	val exit= Unsynchronized.ref  0;

  1150     in

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

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

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

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

  1155 	else

  1156 	    (

  1157 	     exit:=0;

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

  1159 		 (

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

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

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

  1163 		      (

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

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

  1166 		       )

  1167 		  else exit:=1;

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

  1169 		  else (exit:=1)

  1170 		  );

  1171 		 (rev(!q),!r)

  1172 		 )

  1173     end;

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

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

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

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

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

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

  1183 (*. multivariate case .*)

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

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

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

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

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

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

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

  1194 fun mv_deg([]:mv_poly) = 0

  1195   | mv_deg(p1)=

  1196     let

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

  1198     in

  1199 	if length(p1')=0 then 0

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

  1201     end;

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

  1204 fun mv_deg2([]:mv_poly) = 0

  1205   | mv_deg2(p1)=

  1206     let

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

  1208     in

  1209 	if length(p1')=0 then 0

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

  1211     end;

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

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

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

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

  1218 fun uv_content2([]:mv_poly) = 0

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

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

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

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

  1224     let

  1225 	val count= Unsynchronized.ref  0;

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

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

  1228 	val list= Unsynchronized.ref  [];

  1229     in

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

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

  1232 	     else

  1233 		 (

  1234 		  count:=0;

  1235 		  while (!count)<=max do

  1236 		      (

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

  1238 			   (

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

  1240 			    help:=tl(!help)

  1241 			    )

  1242 		       else

  1243 			   (

  1244 			    list:= 0::(!list)

  1245 			    );

  1246 		       count := (!count) + 1

  1247 		       );

  1248 		      (!list)

  1249 		      )

  1250     end;

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

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

  1254   | uv_to_poly p1 =

  1255     let

  1256 	val count= Unsynchronized.ref  0;

  1257 	val help= Unsynchronized.ref  p1;

  1258 	val list= Unsynchronized.ref  [];

  1259     in

  1260 	while length(!help)>0 do

  1261 	    (

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

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

  1264 	     count:=(!count)+1;

  1265 	     help:=tl(!help)

  1266 	     );

  1267 	    (!list)

  1268     end;

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

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

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

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

  1274     let

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

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

  1277     in

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

  1279     end

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

  1281     let

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

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

  1284     in

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

  1286     end

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

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

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

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

  1292     let

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

  1294 	val t= Unsynchronized.ref  [];

  1295 	val y= Unsynchronized.ref  [];

  1296 	val n= Unsynchronized.ref  1;

  1297 	val n1= Unsynchronized.ref [];

  1298     in

  1299 	(

  1300 	 while length(!x)>1 do

  1301 	     (

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

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

  1304 		   else n1:=[];

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

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

  1307 	      x:=tl(!x)

  1308 	      );

  1309 	 (!y)

  1310 	 )

  1311     end;

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

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

  1315   | mv_newton_add [x:mv_poly] t = x

  1316   | mv_newton_add (pl:mv_poly list) t =

  1317     let

  1318 	val expos= Unsynchronized.ref  [];

  1319 	val pll= Unsynchronized.ref  pl;

  1320     in

  1321 	(

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

  1324 	     (

  1325 	      pll:=tl(!pll)

  1326 	      );

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

  1328 	 mv_add(hd(pl),

  1329 		mv_mul(

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

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

  1332 		       LEX_

  1333 		       ),

  1334 		LEX_)

  1335 	 )

  1336     end;

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

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

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

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

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

  1343   | mv_newton pl =

  1344     let

  1345 	val c= Unsynchronized.ref  pl;

  1346 	val c1= Unsynchronized.ref  [];

  1347 	val n=length(pl);

  1348 	val k= Unsynchronized.ref  1;

  1349 	val i= Unsynchronized.ref  n;

  1350 	val ppl= Unsynchronized.ref  [];

  1351     in

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

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

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

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

  1356 	    (

  1357 	     k:=(!k)+1;

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

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

  1360 	     ppl:= !c;

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

  1362 	     );

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

  1364 	c1:=rev(!c1);

  1365 	ppl:= !c1;

  1366 	mv_newton_add (!c1) 1

  1367     end;

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

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

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

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

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

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

  1376     let

  1377 	val y= Unsynchronized.ref  xs;

  1378 	val x= Unsynchronized.ref  [];

  1379     in

  1380 	(

  1381 	 x:=e;

  1382 	 while length(!y)>0 do

  1383 	     (

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

  1385 	      y:=tl(!y)

  1386 	      );

  1387 	 !x

  1388 	 )

  1389     end;

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

  1392 fun list_is_null [] = true

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

  1395 (* check if main variable is zero*)

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

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

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

  1400   | mv_content(p1) =

  1401     let

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

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

  1404 	val result= Unsynchronized.ref  [];

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

  1406     in

  1407 	(

  1408 	 if length(!list)>1 then

  1409 	     (

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

  1411 		  (

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

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

  1415 		   else list:=tl(!list)

  1417 		       );

  1418 		  if length(!list)>0 then

  1419 		   (

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

  1421 		    )

  1422 		  else list:=(!result);

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

  1424 		  (!list)

  1425 		  )

  1426 	 else

  1427 	     (

  1428 	      mv_null3(!list)

  1429 	      )

  1430 	     )

  1431     end

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

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

  1435   | mv_pp(p1) = let

  1436 		    val cont= Unsynchronized.ref  [];

  1437 		    val pp= Unsynchronized.ref [];

  1438 		in

  1439 		    cont:=mv_content(p1);

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

  1441 		    if !pp=[]

  1442 			then error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")

  1443 		    else (!pp)

  1444 		end

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

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

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

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

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

  1451      let

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

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

  1454      in

  1455 	(

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

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

  1458         )

  1459     end

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

  1461 	(

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

  1463 	)

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

  1465 	(

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

  1467         )

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

  1469     let

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

  1471 	val cont =

  1472 		  (

  1473                    if main_zero(mv_content(p1')) andalso

  1474                      (main_zero(mv_content(p2'))) then

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

  1476                    else

  1477                      mv_gcd (mv_content(p1')) (mv_content(p2'))

  1478                   );

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

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

  1481 	val gcd= Unsynchronized.ref  [];

  1482 	val candidate= Unsynchronized.ref  [];

  1483 	val interpolation_list= Unsynchronized.ref  [];

  1484 	val delta= Unsynchronized.ref  [];

  1485         val p1r = Unsynchronized.ref [];

  1486         val p2r = Unsynchronized.ref [];

  1487         val p1r' = Unsynchronized.ref [];

  1488         val p2r' = Unsynchronized.ref [];

  1489 	val factor= Unsynchronized.ref  [];

  1490 	val r= Unsynchronized.ref  0;

  1491 	val gcd_r= Unsynchronized.ref  [];

  1492 	val d= Unsynchronized.ref  0;

  1493 	val exit= Unsynchronized.ref  0;

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

  1495     in

  1496 	(

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

  1498 	     (

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

  1500 	      )

  1501 	 else

  1502 	     (

  1503 	      while !exit=0 do

  1504 		  (

  1505 		   r:=(!r)+1;

  1506                    p1r := mv_lc(p1,LEX_);

  1507 		   p2r := mv_lc(p2,LEX_);

  1508                    if main_zero(!p1r) andalso

  1509                       main_zero(!p2r)

  1510                    then

  1511                        (

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

  1513                        )

  1514                    else

  1515                        (

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

  1517                        );

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

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

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

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

  1522                    then

  1523                        (

  1524 		       )

  1525 		   else

  1526 		       (

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

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

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

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

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

  1532 			if (!d < !current_degree) then

  1533 			    (

  1534 			     current_degree:= !d;

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

  1536 			     )

  1537 			else

  1538 			    (

  1539 			     if (!d = !current_degree) then

  1540 				 (

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

  1542 				  )

  1543 			     else ()

  1544 				 )

  1545 			    );

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

  1547 			  (

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

  1549 			   if !candidate=[] then ()

  1550 			   else

  1551 			       (

  1552 				candidate:=mv_pp(!candidate);

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

  1554 				    (

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

  1556 				     exit:=1

  1557 				     )

  1558 				else ()

  1559 				    );

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

  1561 			       )

  1562 		      else ()

  1563 			  )

  1564 	     );

  1565 	     (!gcd):mv_poly

  1566 	     )

  1567     end;

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

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

  1572     (

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

  1574      );

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

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

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

  1581 fun count_neg ([]:mv_poly) = 0

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

  1583 			  else count_neg xs;

  1585 (*. help function for is_polynomial

  1586     checks the order of the operators .*)

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

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

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

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

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

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

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

  1594   | test_polynomial _ v = false;

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

  1597 fun is_polynomial t = test_polynomial t " ";

  1599 (*. help function for is_expanded

  1600     checks the order of the operators .*)

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

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

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

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

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

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

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

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

  1609   | test_exp  _ v = false;

  1612 (*. help function for check_coeff:

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

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

  1615     let

  1616 	val x= Unsynchronized.ref  NONE;

  1617 	val len= Unsynchronized.ref  0;

  1618 	val vl= Unsynchronized.ref  [];

  1619 	val vh= Unsynchronized.ref  [];

  1620 	val i= Unsynchronized.ref  0;

  1621     in

  1622 	if is_numeral str then

  1623 	    (

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

  1625 		 )

  1626 	else (* variable *)

  1627 	    (

  1628 	     len:=length(v);

  1629 	     vh:=v;

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

  1631 		 (

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

  1633 		      (

  1634 		       vl:=1::(!vl)

  1635 		       )

  1636 		  else

  1637 		      (

  1638 		       vl:=0::(!vl)

  1639 		       );

  1640 		      vh:=tl(!vh);

  1641 		      i:=(!i)+1

  1642 		      );

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

  1644 	    )

  1645     end

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

  1647     let

  1648 	val t1pp= Unsynchronized.ref  [];

  1649 	val t2pp= Unsynchronized.ref  [];

  1650 	val t1c= Unsynchronized.ref  0;

  1651 	val t2c= Unsynchronized.ref  0;

  1652     in

  1653 	(

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

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

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

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

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

  1661 	 )

  1662     end

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

  1664     let

  1665 	val x= Unsynchronized.ref  NONE;

  1666 	val len= Unsynchronized.ref  0;

  1667 	val vl= Unsynchronized.ref  [];

  1668 	val vh= Unsynchronized.ref  [];

  1669 	val vtemp= Unsynchronized.ref  [];

  1670 	val i= Unsynchronized.ref  0;

  1671     in

  1672     (

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

  1674 	 (

  1675 	  len:=length(v);

  1676 	  vh:=v;

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

  1679 	      (

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

  1681 		   (

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

  1683 		    )

  1684 	       else

  1685 		   (

  1686 		    vl:=0::(!vl)

  1687 		    );

  1688 		   vh:=tl(!vh);

  1689 		   i:=(!i)+1

  1690 		   );

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

  1692 	      )

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

  1694 	 )

  1695     end

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

  1697     (

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

  1699 	 )

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

  1701     (

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

  1703 	 )

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

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

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

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

  1709     else false;

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

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

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

  1715 fun mk_monom v' p vs =

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

  1717     in map (conv p) vs end;

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

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

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

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

  1723     if is_numeral str

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

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

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

  1728     let

  1729 	val x= Unsynchronized.ref  NONE;

  1730 	val len= Unsynchronized.ref  0;

  1731 	val vl= Unsynchronized.ref  [];

  1732 	val vh= Unsynchronized.ref  [];

  1733 	val i= Unsynchronized.ref  0;

  1734     in

  1735 	if is_numeral str then

  1736 	    (

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

  1738 		 )

  1739 	else (* variable *)

  1740 	    (

  1741 	     len:=length v;

  1742 	     vh:= v;

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

  1744 		 (

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

  1746 		      (

  1747 		       vl:=1::(!vl)

  1748 		       )

  1749 		  else

  1750 		      (

  1751 		       vl:=0::(!vl)

  1752 		       );

  1753 		      vh:=tl(!vh);

  1754 		      i:=(!i)+1

  1755 		      );

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

  1757 	    )

  1758     end

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

  1760     let

  1761 	val t1pp= Unsynchronized.ref  [];

  1762 	val t2pp= Unsynchronized.ref  [];

  1763 	val t1c= Unsynchronized.ref  0;

  1764 	val t2c= Unsynchronized.ref  0;

  1765     in

  1766 	(

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

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

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

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

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

  1773 	 handle _ => NONE

  1775 	 )

  1776     end

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

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

  1779     let

  1780 	val x= Unsynchronized.ref  NONE;

  1781 	val len= Unsynchronized.ref  0;

  1782 	val vl= Unsynchronized.ref  [];

  1783 	val vh= Unsynchronized.ref  [];

  1784 	val vtemp= Unsynchronized.ref  [];

  1785 	val i= Unsynchronized.ref  0;

  1786     in

  1787     (

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

  1789 	 (

  1790 	  len:=length(v);

  1791 	  vh:=v;

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

  1794 	      (

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

  1796 		   (

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

  1798 		    )

  1799 	       else

  1800 		   (

  1801 		    vl:=0::(!vl)

  1802 		    );

  1803 		   vh:=tl(!vh);

  1804 		   i:=(!i)+1

  1805 		   );

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

  1807 	      )

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

  1809 	 )

  1810     end

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

  1812     (

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

  1814 	 )

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

  1816     (

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

  1818 	 )

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

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

  1822     term2poly

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

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

  1825     			       DIE DU MIR LETZTES MAL GEGEBEN HAST*)

  1826     	 mv_monom list        (*internal representation           *)

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

  1828 .*)

  1829 fun term2poly (t:term) v =

  1830      if is_polynomial t then term2poly' t v

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

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

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

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

  1837 fun expanded2poly (t:term) v =

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

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

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

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

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

  1845 fun powerproduct2term(xs,v) =

  1846     let

  1847 	val xss= Unsynchronized.ref  xs;

  1848 	val vv= Unsynchronized.ref  v;

  1849     in

  1850 	(

  1851 	 while hd(!xss)=0 do

  1852 	     (

  1853 	      xss:=tl(!xss);

  1854 	      vv:=tl(!vv)

  1855 	      );

  1857 	 if list_is_null(tl(!xss)) then

  1858 	     (

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

  1860 	      else

  1861 		  (

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

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

  1864 		   )

  1865 	      )

  1866 	 else

  1867 	     (

  1868 	      if hd(!xss)=1 then

  1869 		  (

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

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

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

  1873 		   )

  1874 	      else

  1875 		  (

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

  1877 		   (

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

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

  1880 		    ) $

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

  1882 		   )

  1883 	      )

  1884 	 )

  1885     end;

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

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

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

  1890     else

  1891 	(

  1892 	 if list_is_null(e) then

  1893 	     (

  1894 	      Free(str_of_int c,HOLogic.realT)

  1895 	      )

  1896 	 else

  1897 	     (

  1898 	      if c=1 then

  1899 		  (

  1900 		   powerproduct2term(e,v)

  1901 		   )

  1902 	      else

  1903 		  (

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

  1905 		   Free(str_of_int c,HOLogic.realT)  $

  1906 		   powerproduct2term(e,v)

  1907 		   )

  1908 		  )

  1909 	     );*)

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

  1913     if list_is_null is

  1914     then

  1915 	if i >= 0

  1916 	then Free (str_of_int i, HOLogic.realT)

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

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

  1919     else

  1920 	if i > 0

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

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

  1923 		   powerproduct2term(is, v)

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

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

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

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

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

  1929     if list_is_null is

  1930     then Free (str_of_int i, HOLogic.realT)

  1931     else if i = 1

  1932     then powerproduct2term (is, vs)

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

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

  1935 	       powerproduct2term (is, vs);

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

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

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

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

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

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

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

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

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

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

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

  1950     else

  1951 	(

  1952 	 if list_is_null(e) then

  1953 	     (

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

  1955 	      )

  1956 	 else

  1957 	     (

  1958 	      if abs(c)=1 then

  1959 		  (

  1960 		   powerproduct2term(e,v)

  1961 		   )

  1962 	      else

  1963 		  (

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

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

  1966 		   powerproduct2term(e,v)

  1967 		   )

  1968 		  )

  1969 	     );

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

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

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

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

  1975     if c1<0 then

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

  1977 	exp2term'(others,vars) $

  1978 	(

  1979 	 monom2term2((c1,e1),vars)

  1980 	 )

  1981     else

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

  1983 	exp2term'(others,vars) $

  1984 	(

  1985 	 monom2term2((c1,e1),vars)

  1986 	 );

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

  1989     a term in polynomial representation .*)

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

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

  1993 fun polynomial2expanded t =

  1994     (let

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

  1996     in

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

  1998     end) handle _ => NONE;

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

  2001 fun expanded2polynomial t =

  2002     (let

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

  2004     in

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

  2006     end) handle _ => NONE;

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

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

  2011     let

  2012 	val p1' = Unsynchronized.ref [];

  2013 	val p2' = Unsynchronized.ref [];

  2014 	val p3  = Unsynchronized.ref []

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

  2016     in

  2017 	(

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

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

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

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

  2022 	     (

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

  2024 	      )

  2025 	 else

  2026 	     (

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

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

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

  2032 	      (

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

  2034 	       $

  2035 	       (

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

  2037 		poly2term(!p1',vars) $

  2038 		poly2term(!p3,vars)

  2039 		)

  2040 	       $

  2041 	       (

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

  2043 		poly2term(!p2',vars) $

  2044 		poly2term(!p3,vars)

  2045 		)

  2046 	       )

  2047 	      else

  2048 	      (

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

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

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

  2052 	       (

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

  2054 		$

  2055 		(

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

  2057 		 poly2term(!p1',vars) $

  2058 		 poly2term(!p3,vars)

  2059 		 )

  2060 		$

  2061 		(

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

  2063 		 poly2term(!p2',vars) $

  2064 		 poly2term(!p3,vars)

  2065 		 )

  2066 		)

  2067 	       )

  2068 	      )

  2069 	     )

  2070     end

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

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

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

  2075     let

  2076 	val p1' = Unsynchronized.ref [];

  2077 	val p2' = Unsynchronized.ref [];

  2078 	val p3  = Unsynchronized.ref []

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

  2080     in

  2081 	(

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

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

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

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

  2087 	     (

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

  2089 	      )

  2090 	 else

  2091 	     (

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

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

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

  2095 	      (

  2096 	       (

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

  2098 		$

  2099 		(

  2100 		 poly2term((!p1'),vars)

  2101 		 )

  2102 		$

  2103 		(

  2104 		 poly2term((!p2'),vars)

  2105 		 )

  2106 		)

  2107 	       ,

  2108 	       if mv_grad(!p3)>0 then

  2109 		   [

  2110 		    (

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

  2112 		     (

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

  2114 		      poly2term((!p3),vars) $

  2115 		      Free("0",HOLogic.realT)

  2116 		      )

  2117 		     )

  2118 		    ]

  2119 	       else

  2120 		   []

  2121 		   )

  2122 	      else

  2123 		  (

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

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

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

  2127 		       (

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

  2129 			$

  2130 			(

  2131 			 poly2term((!p1'),vars)

  2132 			 )

  2133 			$

  2134 			(

  2135 			 poly2term((!p2'),vars)

  2136 			 )

  2137 			,

  2138 			if mv_grad(!p3)>0 then

  2139 			    [

  2140 			     (

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

  2142 			      (

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

  2144 			       poly2term((!p3),vars) $

  2145 			       Free("0",HOLogic.realT)

  2146 			       )

  2147 			      )

  2148 			     ]

  2149 			else

  2150 			    []

  2151 			    )

  2152 		       )

  2153 		  )

  2154 	     )

  2155     end

  2156   | direct_cancel _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");

  2158 (*. same es direct_cancel, this time for expanded forms (input+output).*)

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

  2160     let

  2161 	val p1' = Unsynchronized.ref [];

  2162 	val p2' = Unsynchronized.ref [];

  2163 	val p3  = Unsynchronized.ref []

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

  2165     in

  2166 	(

  2167 	 p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_));

  2168 	 p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_));

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

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

  2172 	     (

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

  2174 	      )

  2175 	 else

  2176 	     (

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

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

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

  2180 	      (

  2181 	       (

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

  2183 		$

  2184 		(

  2185 		 poly2expanded((!p1'),vars)

  2186 		 )

  2187 		$

  2188 		(

  2189 		 poly2expanded((!p2'),vars)

  2190 		 )

  2191 		)

  2192 	       ,

  2193 	       if mv_grad(!p3)>0 then

  2194 		   [

  2195 		    (

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

  2197 		     (

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

  2199 		      poly2expanded((!p3),vars) $

  2200 		      Free("0",HOLogic.realT)

  2201 		      )

  2202 		     )

  2203 		    ]

  2204 	       else

  2205 		   []

  2206 		   )

  2207 	      else

  2208 		  (

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

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

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

  2212 		       (

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

  2214 			$

  2215 			(

  2216 			 poly2expanded((!p1'),vars)

  2217 			 )

  2218 			$

  2219 			(

  2220 			 poly2expanded((!p2'),vars)

  2221 			 )

  2222 			,

  2223 			if mv_grad(!p3)>0 then

  2224 			    [

  2225 			     (

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

  2227 			      (

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

  2229 			       poly2expanded((!p3),vars) $

  2230 			       Free("0",HOLogic.realT)

  2231 			       )

  2232 			      )

  2233 			     ]

  2234 			else

  2235 			    []

  2236 			    )

  2237 		       )

  2238 		  )

  2239 	     )

  2240     end

  2241   | direct_cancel_expanded _ = error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction");

  2244 (*. adds two fractions .*)

  2245 fun add_fract ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =

  2246     let

  2247 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);

  2248 	val t11'= Unsynchronized.ref  (the(term2poly t11 vars));

  2249 (* stopped Test_Isac.thy ...

  2250 val _= tracing"### add_fract: done t11"

  2251 *)

  2252 	val t12'= Unsynchronized.ref  (the(term2poly t12 vars));

  2253 (* stopped Test_Isac.thy ...

  2254 val _= tracing"### add_fract: done t12"

  2255 *)

  2256 	val t21'= Unsynchronized.ref  (the(term2poly t21 vars));

  2257 (* stopped Test_Isac.thy ...

  2258 val _= tracing"### add_fract: done t21"

  2259 *)

  2260 	val t22'= Unsynchronized.ref  (the(term2poly t22 vars));

  2261 (* stopped Test_Isac.thy ...

  2262 val _= tracing"### add_fract: done t22"

  2263 *)

  2264 	val den= Unsynchronized.ref  [];

  2265 	val nom= Unsynchronized.ref  [];

  2266 	val m1= Unsynchronized.ref  [];

  2267 	val m2= Unsynchronized.ref  [];

  2268     in

  2270 	(

  2271 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));

  2272 tracing"### add_fract: done sort mv_lcm";

  2273 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));

  2274 tracing"### add_fract: done sort mv_division t12";

  2275 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));

  2276 tracing"### add_fract: done sort mv_division t22";

  2277 	 nom :=sort (mv_geq LEX_)

  2278 		    (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),

  2279 				       mv_mul(!t21',!m2,LEX_),

  2280 				       LEX_),

  2281 				LEX_));

  2282 tracing"### add_fract: done sort mv_add";

  2283 	 (

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

  2285 	  $

  2286 	  (

  2287 	   poly2term((!nom),vars)

  2288 	   )

  2289 	  $

  2290 	  (

  2291 	   poly2term((!den),vars)

  2292 	   )

  2293 	  )

  2294 	 )

  2295     end

  2296   | add_fract (_,_) = error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call");

  2298 (*. adds two expanded fractions .*)

  2299 fun add_fract_exp ((Const("Fields.inverse_class.divide",_) $ t11 $ t12),(Const("Fields.inverse_class.divide",_) $ t21 $ t22)) =

  2300     let

  2301 	val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22);

  2302 	val t11'= Unsynchronized.ref  (the(expanded2poly t11 vars));

  2303 	val t12'= Unsynchronized.ref  (the(expanded2poly t12 vars));

  2304 	val t21'= Unsynchronized.ref  (the(expanded2poly t21 vars));

  2305 	val t22'= Unsynchronized.ref  (the(expanded2poly t22 vars));

  2306 	val den= Unsynchronized.ref  [];

  2307 	val nom= Unsynchronized.ref  [];

  2308 	val m1= Unsynchronized.ref  [];

  2309 	val m2= Unsynchronized.ref  [];

  2310     in

  2312 	(

  2313 	 den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22'));

  2314 	 m1  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_)));

  2315 	 m2  :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_)));

  2316 	 nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_));

  2317 	 (

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

  2319 	  $

  2320 	  (

  2321 	   poly2expanded((!nom),vars)

  2322 	   )

  2323 	  $

  2324 	  (

  2325 	   poly2expanded((!den),vars)

  2326 	   )

  2327 	  )

  2328 	 )

  2329     end

  2330   | add_fract_exp (_,_) = error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call");

  2332 (*. adds a list of terms .*)

  2333 fun add_list_of_fractions []= (Free("0",HOLogic.realT),[])

  2334   | add_list_of_fractions [x]= direct_cancel x

  2335   | add_list_of_fractions (x::y::xs) =

  2336     let

  2337 	val (t1a,rest1)=direct_cancel(x);

  2338 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(x)";

  2339 	val (t2a,rest2)=direct_cancel(y);

  2340 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(y)";

  2341 	val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs));

  2342 val _= tracing"### add_list_of_fractions xs: has done add_list_of_fraction xs";

  2343 	val (t4a,rest4)=direct_cancel(t3a);

  2344 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(t3a)";

  2345 	val rest=rest1 union rest2 union rest3 union rest4;

  2346     in

  2347 	(tracing"### add_list_of_fractions in";

  2348 	 (

  2349 	 (t4a,rest)

  2350 	 )

  2351 	 )

  2352     end;

  2354 (*. adds a list of expanded terms .*)

  2355 fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])

  2356   | add_list_of_fractions_exp [x]= direct_cancel_expanded x

  2357   | add_list_of_fractions_exp (x::y::xs) =

  2358     let

  2359 	val (t1a,rest1)=direct_cancel_expanded(x);

  2360 	val (t2a,rest2)=direct_cancel_expanded(y);

  2361 	val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));

  2362 	val (t4a,rest4)=direct_cancel_expanded(t3a);

  2363 	val rest=rest1 union rest2 union rest3 union rest4;

  2364     in

  2365 	(

  2366 	 (t4a,rest)

  2367 	 )

  2368     end;

  2370 (*. calculates the lcm of a list of mv_poly .*)

  2371 fun calc_lcm ([x],var)= (x,var)

  2372   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2374 (*. converts a list of terms to a list of mv_poly .*)

  2375 fun t2d([],_)=[]

  2376   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2378 (*. same as t2d, this time for expanded forms .*)

  2379 fun t2d_exp([],_)=[]

  2380   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2382 (*. converts a list of fract terms to a list of their denominators .*)

  2383 fun termlist2denominators [] = ([],[])

  2384   | termlist2denominators xs =

  2385     let

  2386 	val xxs= Unsynchronized.ref  xs;

  2387 	val var= Unsynchronized.ref  [];

  2388     in

  2389 	var:=[];

  2390 	while length(!xxs)>0 do

  2391 	    (

  2392 	     let

  2393 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2394 	     in

  2395 		 (

  2396 		  xxs:=tl(!xxs);

  2397 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2398 		  )

  2399 	     end

  2400 	     );

  2401 	    (t2d(xs,!var),!var)

  2402     end;

  2404 (*. calculates the lcm of a list of mv_poly .*)

  2405 fun calc_lcm ([x],var)= (x,var)

  2406   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2408 (*. converts a list of terms to a list of mv_poly .*)

  2409 fun t2d([],_)=[]

  2410   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2412 (*. same as t2d, this time for expanded forms .*)

  2413 fun t2d_exp([],_)=[]

  2414   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2416 (*. converts a list of fract terms to a list of their denominators .*)

  2417 fun termlist2denominators [] = ([],[])

  2418   | termlist2denominators xs =

  2419     let

  2420 	val xxs= Unsynchronized.ref  xs;

  2421 	val var= Unsynchronized.ref  [];

  2422     in

  2423 	var:=[];

  2424 	while length(!xxs)>0 do

  2425 	    (

  2426 	     let

  2427 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2428 	     in

  2429 		 (

  2430 		  xxs:=tl(!xxs);

  2431 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2432 		  )

  2433 	     end

  2434 	     );

  2435 	    (t2d(xs,!var),!var)

  2436     end;

  2438 (*. same as termlist2denminators, this time for expanded forms .*)

  2439 fun termlist2denominators_exp [] = ([],[])

  2440   | termlist2denominators_exp xs =

  2441     let

  2442 	val xxs= Unsynchronized.ref  xs;

  2443 	val var= Unsynchronized.ref  [];

  2444     in

  2445 	var:=[];

  2446 	while length(!xxs)>0 do

  2447 	    (

  2448 	     let

  2449 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2450 	     in

  2451 		 (

  2452 		  xxs:=tl(!xxs);

  2453 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2454 		  )

  2455 	     end

  2456 	     );

  2457 	    (t2d_exp(xs,!var),!var)

  2458     end;

  2460 (*. reduces all fractions to the least common denominator .*)

  2461 fun com_den(x::xs,denom,den,var)=

  2462     let

  2463 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2464 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2465 	val p3= #1(mv_division(denom,p2,LEX_));

  2466 	val p1var=get_vars(p1');

  2467     in

  2468 	if length(xs)>0 then

  2469 	    if p3=[(1,mv_null2(var))] then

  2470 		(

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

  2472 		 $

  2473 		 (

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

  2475 		  $

  2476 		  poly2term(the (term2poly p1' p1var),p1var)

  2477 		  $

  2478 		  den

  2479 		  )

  2480 		 $

  2481 		 #1(com_den(xs,denom,den,var))

  2482 		,

  2483 		[]

  2484 		)

  2485 	    else

  2486 		(

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

  2488 		 $

  2489 		 (

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

  2491 		  $

  2492 		  (

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

  2494 		   poly2term(the (term2poly p1' p1var),p1var) $

  2495 		   poly2term(p3,var)

  2496 		   )

  2497 		  $

  2498 		  (

  2499 		   den

  2500 		   )

  2501 		  )

  2502 		 $

  2503 		 #1(com_den(xs,denom,den,var))

  2504 		,

  2505 		[]

  2506 		)

  2507 	else

  2508 	    if p3=[(1,mv_null2(var))] then

  2509 		(

  2510 		 (

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

  2512 		  $

  2513 		  poly2term(the (term2poly p1' p1var),p1var)

  2514 		  $

  2515 		  den

  2516 		  )

  2517 		 ,

  2518 		 []

  2519 		 )

  2520 	     else

  2521 		 (

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

  2523 		  $

  2524 		  (

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

  2526 		   poly2term(the (term2poly p1' p1var),p1var) $

  2527 		   poly2term(p3,var)

  2528 		   )

  2529 		  $

  2530 		  den

  2531 		  ,

  2532 		  []

  2533 		  )

  2534     end;

  2536 (*. same as com_den, this time for expanded forms .*)

  2537 fun com_den_exp(x::xs,denom,den,var)=

  2538     let

  2539 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2540 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2541 	val p3= #1(mv_division(denom,p2,LEX_));

  2542 	val p1var=get_vars(p1');

  2543     in

  2544 	if length(xs)>0 then

  2545 	    if p3=[(1,mv_null2(var))] then

  2546 		(

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

  2548 		 $

  2549 		 (

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

  2551 		  $

  2552 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2553 		  $

  2554 		  den

  2555 		  )

  2556 		 $

  2557 		 #1(com_den_exp(xs,denom,den,var))

  2558 		,

  2559 		[]

  2560 		)

  2561 	    else

  2562 		(

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

  2564 		 $

  2565 		 (

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

  2567 		  $

  2568 		  (

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

  2570 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2571 		   poly2expanded(p3,var)

  2572 		   )

  2573 		  $

  2574 		  (

  2575 		   den

  2576 		   )

  2577 		  )

  2578 		 $

  2579 		 #1(com_den_exp(xs,denom,den,var))

  2580 		,

  2581 		[]

  2582 		)

  2583 	else

  2584 	    if p3=[(1,mv_null2(var))] then

  2585 		(

  2586 		 (

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

  2588 		  $

  2589 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2590 		  $

  2591 		  den

  2592 		  )

  2593 		 ,

  2594 		 []

  2595 		 )

  2596 	     else

  2597 		 (

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

  2599 		  $

  2600 		  (

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

  2602 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2603 		   poly2expanded(p3,var)

  2604 		   )

  2605 		  $

  2606 		  den

  2607 		  ,

  2608 		  []

  2609 		  )

  2610     end;

  2612 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2613 -------------------------------------------------------------

  2614 (* WN0210???SK brauch ma des überhaupt *)

  2615 fun com_den2(x::xs,denom,den,var)=

  2616     let

  2617 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2618 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2619 	val p3= #1(mv_division(denom,p2,LEX_));

  2620 	val p1var=get_vars(p1');

  2621     in

  2622 	if length(xs)>0 then

  2623 	    if p3=[(1,mv_null2(var))] then

  2624 		(

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

  2626 		 poly2term(the(term2poly p1' p1var),p1var) $

  2627 		 com_den2(xs,denom,den,var)

  2628 		)

  2629 	    else

  2630 		(

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

  2632 		 (

  2633 		   let

  2634 		       val p3'=poly2term(p3,var);

  2635 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2636 		   in

  2637 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2638 		   end

  2639 		  ) $

  2640 		 com_den2(xs,denom,den,var)

  2641 		)

  2642 	else

  2643 	    if p3=[(1,mv_null2(var))] then

  2644 		(

  2645 		 poly2term(the(term2poly p1' p1var),p1var)

  2646 		 )

  2647 	     else

  2648 		 (

  2649 		   let

  2650 		       val p3'=poly2term(p3,var);

  2651 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2652 		   in

  2653 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2654 		   end

  2655 		  )

  2656     end;

  2658 (* WN0210???SK brauch ma des überhaupt *)

  2659 fun com_den_exp2(x::xs,denom,den,var)=

  2660     let

  2661 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2662 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2663 	val p3= #1(mv_division(denom,p2,LEX_));

  2664 	val p1var=get_vars p1';

  2665     in

  2666 	if length(xs)>0 then

  2667 	    if p3=[(1,mv_null2(var))] then

  2668 		(

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

  2670 		 poly2expanded(the (expanded2poly p1' p1var),p1var) $

  2671 		 com_den_exp2(xs,denom,den,var)

  2672 		)

  2673 	    else

  2674 		(

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

  2676 		 (

  2677 		   let

  2678 		       val p3'=poly2expanded(p3,var);

  2679 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2680 		   in

  2681 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2682 		   end

  2683 		  ) $

  2684 		 com_den_exp2(xs,denom,den,var)

  2685 		)

  2686 	else

  2687 	    if p3=[(1,mv_null2(var))] then

  2688 		(

  2689 		 poly2expanded(the (expanded2poly p1' p1var),p1var)

  2690 		 )

  2691 	     else

  2692 		 (

  2693 		   let

  2694 		       val p3'=poly2expanded(p3,var);

  2695 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2696 		   in

  2697 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2698 		   end

  2699 		  )

  2700     end;

  2701 ---------------------------------------------------------*)

  2704 (*. searches for an element y of a list ys, which has an gcd not 1 with x .*)

  2705 fun exists_gcd (x,[]) = false

  2706   | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then  exists_gcd (x,ys)

  2707 			   else true;

  2709 (*. divides each element of the list xs with y .*)

  2710 fun list_div ([],y) = []

  2711   | list_div (x::xs,y) =

  2712     let

  2713 	val (d,r)=mv_division(x,y,LEX_);

  2714     in

  2715 	if r=[] then

  2716 	    d::list_div(xs,y)

  2717 	else x::list_div(xs,y)

  2718     end;

  2720 (*. checks if x is in the list ys .*)

  2721 fun in_list (x,[]) = false

  2722   | in_list (x,y::ys) = if x=y then true

  2723 			else in_list(x,ys);

  2725 (*. deletes all equal elements of the list xs .*)

  2726 fun kill_equal [] = []

  2727   | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)

  2728 			 else x::kill_equal(xs);

  2730 (*. searches for new factors .*)

  2731 fun new_factors [] = []

  2732   | new_factors (list:mv_poly list):mv_poly list =

  2733     let

  2734 	val l = kill_equal list;

  2735 	val len = length(l);

  2736     in

  2737 	if len>=2 then

  2738 	    (

  2739 	     let

  2740 		 val x::y::xs=l;

  2741 		 val gcd=mv_gcd x y;

  2742 	     in

  2743 		 if gcd=[(1,mv_null2(#2(hd(x))))] then

  2744 		     (

  2745 		      if exists_gcd(x,xs) then new_factors (y::xs @ [x])

  2746 		      else x::new_factors(y::xs)

  2747 	             )

  2748 		 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))

  2749 	     end

  2750 	     )

  2751 	else

  2752 	    if len=1 then [hd(l)]

  2753 	    else []

  2754     end;

  2756 (*. gets the factors of a list .*)

  2757 fun get_factors x = new_factors x;

  2759 (*. multiplies the elements of the list .*)

  2760 fun multi_list [] = []

  2761   | multi_list (x::xs) = if xs=[] then x

  2762 			 else mv_mul(x,multi_list xs,LEX_);

  2764 (*. makes a term out of the elements of the list (polynomial representation) .*)

  2765 fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2766   | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)

  2767 			       else

  2768 				   (

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

  2770 				    poly2term(sort (mv_geq LEX_) (x),vars) $

  2771 				    make_term(xs,vars)

  2772 				    );

  2774 (*. factorizes the denominator (polynomial representation) .*)

  2775 fun factorize_den (l,den,vars) =

  2776     let

  2777 	val factor_list=kill_equal( (get_factors l));

  2778 	val mlist=multi_list(factor_list);

  2779 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2780     in

  2781 	if rest=[] then

  2782 	    (

  2783 	     if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)

  2784 	     else make_term(last::factor_list,vars)

  2785 	     )

  2786 	else error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")

  2787     end;

  2789 (*. makes a term out of the elements of the list (expanded polynomial representation) .*)

  2790 fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2791   | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)

  2792 			       else

  2793 				   (

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

  2795 				    poly2expanded(sort (mv_geq LEX_) (x),vars) $

  2796 				    make_exp(xs,vars)

  2797 				    );

  2799 (*. factorizes the denominator (expanded polynomial representation) .*)

  2800 fun factorize_den_exp (l,den,vars) =

  2801     let

  2802 	val factor_list=kill_equal( (get_factors l));

  2803 	val mlist=multi_list(factor_list);

  2804 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2805     in

  2806 	if rest=[] then

  2807 	    (

  2808 	     if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)

  2809 	     else make_exp(last::factor_list,vars)

  2810 	     )

  2811 	else error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")

  2812     end;

  2814 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2815 (*. the nominators and denominators with the correct factor .*)

  2816 (*. (polynomial representation) .*)

  2817 fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)

  2818   | step_add_list_of_fractions [x]= error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")

  2819   | step_add_list_of_fractions (xs) =

  2820     let

  2821         val den_list=termlist2denominators (xs); (* list of denominators *)

  2822 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2823 	val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2824     in

  2825 	com_den(xs,denom,den,var)

  2826     end;

  2828 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2829 (*. the nominators and denominators with the correct factor .*)

  2830 (*. (expanded polynomial representation) .*)

  2831 fun step_add_list_of_fractions_exp []  = (Free("0",HOLogic.realT),[]:term list)

  2832   | step_add_list_of_fractions_exp [x] = error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")

  2833   | step_add_list_of_fractions_exp (xs)=

  2834     let

  2835         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2836 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2837 	val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2838     in

  2839 	com_den_exp(xs,denom,den,var)

  2840     end;

  2842 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2843 -------------------------------------------------------------

  2844 (* WN0210???SK brauch ma des überhaupt *)

  2845 fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)

  2846   | step_add_list_of_fractions2 [x]=(x,[])

  2847   | step_add_list_of_fractions2 (xs) =

  2848     let

  2849         val den_list=termlist2denominators (xs); (* list of denominators *)

  2850 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2851 	val den=factorize_den(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2852     in

  2853 	(

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

  2855 	 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $

  2856 	 poly2term(denom,var)

  2857 	,

  2858 	[]

  2859 	)

  2860     end;

  2862 (* WN0210???SK brauch ma des überhaupt *)

  2863 fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)

  2864   | step_add_list_of_fractions2_exp [x]=(x,[])

  2865   | step_add_list_of_fractions2_exp (xs) =

  2866     let

  2867         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2868 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2869 	val den=factorize_den_exp(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2870     in

  2871 	(

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

  2873 	 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $

  2874 	 poly2expanded(denom,var)

  2875 	,

  2876 	[]

  2877 	)

  2878     end;

  2879 ---------------------------------------------- *)

  2882 (* converts a term, which contains several terms seperated by +, into a list of these terms .*)

  2883 fun term2list (t as (Const("Fields.inverse_class.divide",_) $ _ $ _)) = [t]

  2884   | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =

  2885       [Const ("Fields.inverse_class.divide",

  2886         [HOLogic.realT,HOLogic.realT] ---> HOLogic.realT) $

  2887 	  t $ Free("1",HOLogic.realT)

  2888      ]

  2889   | term2list (t as (Free(_,_))) =

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

  2891 	  t $  Free("1",HOLogic.realT)

  2892      ]

  2893   | term2list (t as (Const("Groups.times_class.times",_) $ _ $ _)) =

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

  2895 	  t $ Free("1",HOLogic.realT)

  2896      ]

  2897   | term2list (Const("Groups.plus_class.plus",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)

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

  2899     error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")

  2900   | term2list _ = error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");

  2902 (*.factors out the gcd of nominator and denominator:

  2903    a/b = (a' * gcd)/(b' * gcd),  a,b,gcd  are poly[2].*)

  2905 (*. brings the term into a normal form .*)

  2906 fun norm_rational_ (thy:theory) t =

  2907     SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;

  2908 fun norm_expanded_rat_ (thy:theory) t =

  2909     SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;

  2912 (*.evaluates conditions in calculate_Rational.*)

  2913 (*make local with FIXX@ME result:term *term list*)

  2914 val calc_rat_erls = prep_rls(

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

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

  2917 	 rules =

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

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

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

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

  2922 	  ],

  2923 	 scr = EmptyScr});

  2926 (*.simplifies expressions with numerals;

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

  2928    need to have constants to be commuted together respectively.*)

  2929 val calculate_Rational = prep_rls (merge_rls "calculate_Rational"

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

  2931 	    erls = calc_rat_erls, srls = Erls,

  2932 	    calc = [], errpatts = [],

  2933 	    rules =

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  2971 	       ], scr = EmptyScr})

  2972 	  calculate_Poly);

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

  2975 fun eval_is_expanded (thmid:string) _

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

  2977     if is_expanded arg

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

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

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

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

  2982   | eval_is_expanded _ _ _ _ = NONE;

  2984 val rational_erls =

  2985     merge_rls "rational_erls" calculate_Rational

  2986 	      (append_rls "is_expanded" Atools_erls

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

  2988 			   ]);

  2991 (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:

  2992  =================================================================

  2993  A[2] 'cancel_p': .

  2994  A[3] 'cancel': .

  2995  B[2] 'common_nominator_p': transforms summands in a term [2]

  2996          to fractions with the (least) common multiple as nominator.

  2997  B[3] 'norm_rational': normalizes arbitrary algebraic terms (without

  2998          radicals and transzendental functions) to one canceled fraction,

  2999 	 nominator and denominator in polynomial form.

  3001 In order to meet isac's requirements for interactive and stepwise calculation,

  3002 each 'reverse-rewerite-set' consists of an initialization for the interpreter

  3003 state and of 4 functions, each of which employs rewriting as much as possible.

  3004 The signature of these functions are the same in each 'reverse-rewrite-set'

  3005 respectively.*)

  3007 (* ************************************************************************* *)

  3009 local(*. cancel_p

  3010 ------------------------

  3011 cancels a single fraction consisting of two (uni- or multivariate)

  3012 polynomials WN0609???SK[2] into another such a fraction; examples:

  3014 	   a^2 + -1*b^2         a + b

  3015         -------------------- = ---------

  3016 	a^2 + -2*a*b + b^2     a + -1*b

  3018         a^2    a

  3019         --- = ---

  3020          a     1

  3022 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)

  3023 (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)

  3025 val {rules, rew_ord=(_,ro),...} =

  3026     rep_rls (assoc_rls "make_polynomial");

  3027 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3028            see rational.sml --- investigate rulesets for cancel_p ---*)

  3029 val {rules, rew_ord=(_,ro),...} =

  3030     rep_rls (assoc_rls "rev_rew_p");

  3032 (*.init_state = fn : term -> istate

  3033 initialzies the state of the script interpreter. The state is:

  3035 type rrlsstate =      (*state for reverse rewriting*)

  3036      (term *          (*the current formula*)

  3037       term *          (*the final term*)

  3038       rule list       (*'reverse rule list' (#)*)

  3039 	    list *    (*may be serveral, eg. in norm_rational*)

  3040       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3041        (term *        (*... rewrite with ...*)

  3042 	term list))   (*... assumptions*)

  3043 	  list);      (*derivation from given term to normalform

  3044 		       in reverse order with sym_thm;

  3045                        (#) could be extracted from here by (map #1)*).*)

  3046 (* val {rules, rew_ord=(_,ro),...} =

  3047        rep_rls (assoc_rls "rev_rew_p")        (*USE ALWAYS, SEE val cancel_p*);

  3048    val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*);

  3049    val t = t;

  3050    *)

  3051 fun init_state thy eval_rls ro t =

  3052     let val SOME (t',_) = factout_p_ thy t

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

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

  3055         val der = der @ [(Thm ("real_mult_div_cancel2",

  3056 			       num_str @{thm real_mult_div_cancel2}),

  3057 			  (t'',asm))]

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

  3059 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3060 				      "sym_real_mult_0",

  3061 				      "sym_real_mult_1"

  3062 				      (*..insufficient,eg.make_Polynomial*)])rs

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

  3065 (*.locate_rule = fn : rule list -> term -> rule

  3066 		      -> (rule * (term * term list) option) list.

  3067   checks a rule R for being a cancel-rule, and if it is,

  3068   then return the list of rules (+ the terms they are rewriting to)

  3069   which need to be applied before R should be applied.

  3070   precondition: the rule is applicable to the argument-term.

  3071 arguments:

  3072   rule list: the reverse rule list

  3073   -> term  : ... to which the rule shall be applied

  3074   -> rule  : ... to be applied to term

  3075 value:

  3076   -> (rule           : a rule rewriting to ...

  3077       * (term        : ... the resulting term ...

  3078          * term list): ... with the assumptions ( //#0).

  3079       ) list         : there may be several such rules;

  3080 		       the list is empty, if the rule has nothing to do

  3081 		       with cancelation.*)

  3082 (* val () = ();

  3083    *)

  3084 fun locate_rule thy eval_rls ro [rs] t r =

  3085     if (id_of_thm r) mem (map (id_of_thm)) rs

  3086     then let val ropt =

  3087 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3088 	 in case ropt of

  3089 		SOME ta => [(r, ta)]

  3090 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3091 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3092 			 []) end

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

  3094   | locate_rule _ _ _ _ _ _ =

  3095     error ("locate_rule: doesnt match rev-sets in istate");

  3097 (*.next_rule = fn : rule list -> term -> rule option

  3098   for a given term return the next rules to be done for cancelling.

  3099 arguments:

  3100   rule list     : the reverse rule list

  3101   term          : the term for which ...

  3102 value:

  3103   -> rule option: ... this rule is appropriate for cancellation;

  3104 		  there may be no such rule (if the term is canceled already.*)

  3105 (* val thy = thy;

  3106    val Rrls {rew_ord=(_,ro),...} = cancel;

  3107    val ([rs],t) = (rss,f);

  3108    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3110    val (thy, [rs]) = (thy, revsets);

  3111    val Rrls {rew_ord=(_,ro),...} = cancel;

  3112    nex [rs] t;

  3113    *)

  3114 fun next_rule thy eval_rls ro [rs] t =

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

  3116     in case der of

  3117 (* val (_,r,_)::_ = der;

  3118    *)

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

  3120 	 | _ => NONE

  3121     end

  3122   | next_rule _ _ _ _ _ =

  3123     error ("next_rule: doesnt match rev-sets in istate");

  3125 (*.val attach_form = f : rule list -> term -> term

  3126 			 -> (rule * (term * term list)) list

  3127   checks an input term TI, if it may belong to a current cancellation, by

  3128   trying to derive it from the given term TG.

  3129 arguments:

  3130   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3131   -> term: TI, the next one input by the user

  3132 value:

  3133   -> (rule           : the rule to be applied in order to reach TI

  3134       * (term        : ... obtained by applying the rule ...

  3135          * term list): ... and the respective assumptions.

  3136       ) list         : there may be several such rules;

  3137                        the list is empty, if the users term does not belong

  3138 		       to a cancellation of the term last agreed upon.*)

  3139 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

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

  3142 in

  3144 val cancel_p =

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

  3146 	  rew_ord=("ord_make_polynomial",

  3147 		   ord_make_polynomial false thy),

  3148 	  erls = rational_erls,

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

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

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

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

  3153 	  errpatts = [],

  3154 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3155 		     normal_form = cancel_p_ thy,

  3156 		     locate_rule = locate_rule thy Atools_erls ro,

  3157 		     next_rule   = next_rule thy Atools_erls ro,

  3158 		     attach_form = attach_form}}

  3159 end;(*local*)

  3161 local(*.ad [2] 'common_nominator_p'

  3162 ---------------------------------

  3163 FIXME Beschreibung .*)

  3166 val {rules=rules,rew_ord=(_,ro),...} =

  3167     rep_rls (assoc_rls "make_polynomial");

  3168 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3169            see rational.sml --- investigate rulesets for cancel_p ---*)

  3170 val {rules, rew_ord=(_,ro),...} =

  3171     rep_rls (assoc_rls "rev_rew_p");

  3172 val thy = thy;

  3175 (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option

  3176   as defined above*)

  3178 (*.init_state = fn : term -> istate

  3179 initialzies the state of the interactive interpreter. The state is:

  3181 type rrlsstate =      (*state for reverse rewriting*)

  3182      (term *          (*the current formula*)

  3183       term *          (*the final term*)

  3184       rule list       (*'reverse rule list' (#)*)

  3185 	    list *    (*may be serveral, eg. in norm_rational*)

  3186       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3187        (term *        (*... rewrite with ...*)

  3188 	term list))   (*... assumptions*)

  3189 	  list);      (*derivation from given term to normalform

  3190 		       in reverse order with sym_thm;

  3191                        (#) could be extracted from here by (map #1)*).*)

  3192 fun init_state thy eval_rls ro t =

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

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

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

  3196         val der = der @ [(Thm ("real_mult_div_cancel2",

  3197 			       num_str @{thm real_mult_div_cancel2}),

  3198 			  (t'',asm))]

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

  3200 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3201 				      "sym_real_mult_0",

  3202 				      "sym_real_mult_1"]) rs;

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

  3205 (* use"knowledge/Rational.ML";

  3206    *)

  3208 (*.locate_rule = fn : rule list -> term -> rule

  3209 		      -> (rule * (term * term list) option) list.

  3210   checks a rule R for being a cancel-rule, and if it is,

  3211   then return the list of rules (+ the terms they are rewriting to)

  3212   which need to be applied before R should be applied.

  3213   precondition: the rule is applicable to the argument-term.

  3214 arguments:

  3215   rule list: the reverse rule list

  3216   -> term  : ... to which the rule shall be applied

  3217   -> rule  : ... to be applied to term

  3218 value:

  3219   -> (rule           : a rule rewriting to ...

  3220       * (term        : ... the resulting term ...

  3221          * term list): ... with the assumptions ( //#0).

  3222       ) list         : there may be several such rules;

  3223 		       the list is empty, if the rule has nothing to do

  3224 		       with cancelation.*)

  3225 (* val () = ();

  3226    *)

  3227 fun locate_rule thy eval_rls ro [rs] t r =

  3228     if (id_of_thm r) mem (map (id_of_thm)) rs

  3229     then let val ropt =

  3230 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3231 	 in case ropt of

  3232 		SOME ta => [(r, ta)]

  3233 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3234 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3235 			 []) end

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

  3237   | locate_rule _ _ _ _ _ _ =

  3238     error ("locate_rule: doesnt match rev-sets in istate");

  3240 (*.next_rule = fn : rule list -> term -> rule option

  3241   for a given term return the next rules to be done for cancelling.

  3242 arguments:

  3243   rule list     : the reverse rule list

  3244   term          : the term for which ...

  3245 value:

  3246   -> rule option: ... this rule is appropriate for cancellation;

  3247 		  there may be no such rule (if the term is canceled already.*)

  3248 (* val thy = thy;

  3249    val Rrls {rew_ord=(_,ro),...} = cancel;

  3250    val ([rs],t) = (rss,f);

  3251    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3253    val (thy, [rs]) = (thy, revsets);

  3254    val Rrls {rew_ord=(_,ro),...} = cancel;

  3255    nex [rs] t;

  3256    *)

  3257 fun next_rule thy eval_rls ro [rs] t =

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

  3259     in case der of

  3260 (* val (_,r,_)::_ = der;

  3261    *)

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

  3263 	 | _ => NONE

  3264     end

  3265   | next_rule _ _ _ _ _ =

  3266     error ("next_rule: doesnt match rev-sets in istate");

  3268 (*.val attach_form = f : rule list -> term -> term

  3269 			 -> (rule * (term * term list)) list

  3270   checks an input term TI, if it may belong to a current cancellation, by

  3271   trying to derive it from the given term TG.

  3272 arguments:

  3273   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3274   -> term: TI, the next one input by the user

  3275 value:

  3276   -> (rule           : the rule to be applied in order to reach TI

  3277       * (term        : ... obtained by applying the rule ...

  3278          * term list): ... and the respective assumptions.

  3279       ) list         : there may be several such rules;

  3280                        the list is empty, if the users term does not belong

  3281 		       to a cancellation of the term last agreed upon.*)

  3282 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

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

  3285 (* if each pat* matches with the current term, the Rrls is applied

  3286    (there are no preconditions to be checked, they are @{term True}) *)

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

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

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

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

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

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

  3293 in

  3295 (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;

  3296   besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht

  3297   dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)

  3298 val common_nominator_p =

  3299     Rrls {id = "common_nominator_p", prepat=prepat,

  3300 	  rew_ord=("ord_make_polynomial",

  3301 		   ord_make_polynomial false thy),

  3302 	  erls = rational_erls,

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

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

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

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

  3307 	  errpatts = [],

  3308 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3309 		     normal_form = add_fraction_p_ thy,(*FIXME.WN0211*)

  3310 		     locate_rule = locate_rule thy Atools_erls ro,

  3311 		     next_rule   = next_rule thy Atools_erls ro,

  3312 		     attach_form = attach_form}}

  3313 end;(*local*)

  3315 end; (*struct*)

  3317 *}

  3318 ML {*

  3319 open RationalI;

  3320 (*##*)

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

  3323 fun is_ratpolyexp (Free _) = true

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

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

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

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

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

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

  3330                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3331   | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) =

  3332                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3333   | is_ratpolyexp (Const ("Groups.times_class.times",_) $ t1 $ t2) =

  3334                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3335   | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =

  3336                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3337   | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ t1 $ t2) =

  3338                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3339   | is_ratpolyexp _ = false;

  3341 (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)

  3342 fun eval_is_ratpolyexp (thmid:string) _

  3343 		       (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =

  3344     if is_ratpolyexp arg

  3345     then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  3346 	         Trueprop $ (mk_equality (t, @{term True})))

  3347     else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  3348 	         Trueprop $ (mk_equality (t, @{term False})))

  3349   | eval_is_ratpolyexp _ _ _ _ = NONE;

  3351 (*("get_denominator", ("Rational.get_denominator", eval_get_denominator ""))*)

  3352 fun eval_get_denominator (thmid:string) _

  3353 		      (t as Const ("Rational.get_denominator", _) $

  3354               (Const ("Fields.inverse_class.divide", _) $ num $

  3355                 denom)) thy =

  3356       SOME (mk_thmid thmid "" (term_to_string''' thy denom) "",

  3357 	            Trueprop $ (mk_equality (t, denom)))

  3358   | eval_get_denominator _ _ _ _ = NONE;

  3360 (*("get_numerator", ("Rational.get_numerator", eval_get_numerator ""))*)

  3361 fun eval_get_numerator (thmid:string) _

  3362       (t as Const ("Rational.get_numerator", _) $

  3363           (Const ("Fields.inverse_class.divide", _) $num

  3364             $denom )) thy =

  3365     SOME (mk_thmid thmid "" (term_to_string''' thy num) "",

  3366 	    Trueprop $ (mk_equality (t, num)))

  3367   | eval_get_numerator _ _ _ _ = NONE;

  3369 (*-------------------18.3.03 --> struct <-----------vvv--*)

  3370 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)

  3372 (*WN100906 removed in favour of discard_minus = discard_minus_ formerly

  3373    discard binary minus, shift unary minus into -1*;

  3374    unary minus before numerals are put into the numeral by parsing;

  3375    contains absolute minimum of thms for context in norm_Rational

  3376 *** val discard_minus = prep_rls(

  3377   Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3378       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3379       rules =

  3380         [Thm ("real_diff_minus", num_str @{thm real_diff_minus}),

  3381 	           (*"a - b = a + -1 * b"*)

  3382 	         Thm ("sym_real_mult_minus1", num_str (@{thm real_mult_minus1} RS @{thm sym}))

  3383 	           (*- ?z = "-1 * ?z"*)],

  3384       scr = EmptyScr

  3385       }):rls;

  3386  *)

  3388 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)

  3389 val powers_erls = prep_rls(

  3390   Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3391       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3392       rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),

  3393 	       Calc ("Atools.is'_even",eval_is_even "#is_even_"),

  3394 	       Calc ("Orderings.ord_class.less",eval_equ "#less_"),

  3395 	       Thm ("not_false", num_str @{thm not_false}),

  3396 	       Thm ("not_true", num_str @{thm not_true}),

  3397 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

  3398 	       ],

  3399       scr = EmptyScr

  3400       }:rls);

  3401 (*.all powers over + distributed; atoms over * collected, other distributed

  3402    contains absolute minimum of thms for context in norm_Rational .*)

  3403 val powers = prep_rls(

  3404   Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3405       erls = powers_erls, srls = Erls, calc = [], errpatts = [],

  3406       rules = [Thm ("realpow_multI", num_str @{thm realpow_multI}),

  3407 	       (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)

  3408 	       Thm ("realpow_pow",num_str @{thm realpow_pow}),

  3409 	       (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)

  3410 	       Thm ("realpow_oneI",num_str @{thm realpow_oneI}),

  3411 	       (*"r ^^^ 1 = r"*)

  3412 	       Thm ("realpow_minus_even",num_str @{thm realpow_minus_even}),

  3413 	       (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)

  3414 	       Thm ("realpow_minus_odd",num_str @{thm realpow_minus_odd}),

  3415 	       (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)

  3417 	       (*----- collect atoms over * -----*)

  3418 	       Thm ("realpow_two_atom",num_str @{thm realpow_two_atom}),

  3419 	       (*"r is_atom ==> r * r = r ^^^ 2"*)

  3420 	       Thm ("realpow_plus_1",num_str @{thm realpow_plus_1}),

  3421 	       (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)

  3422 	       Thm ("realpow_addI_atom",num_str @{thm realpow_addI_atom}),

  3423 	       (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)

  3425 	       (*----- distribute none-atoms -----*)

  3426 	       Thm ("realpow_def_atom",num_str @{thm realpow_def_atom}),

  3427 	       (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)

  3428 	       Thm ("realpow_eq_oneI",num_str @{thm realpow_eq_oneI}),

  3429 	       (*"1 ^^^ n = 1"*)

  3430 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

  3431 	       ],

  3432       scr = EmptyScr

  3433       }:rls);

  3434 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3435 val rat_mult_divide = prep_rls(

  3436   Rls {id = "rat_mult_divide", preconds = [],

  3437        rew_ord = ("dummy_ord",dummy_ord),

  3438       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3439       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),

  3440 	       (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3441 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

  3442 	       (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3443 	       otherwise inv.to a / b / c = ...*)

  3444 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

  3445 	       (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much

  3446 		     and does not commute a / b * c ^^^ 2 !*)

  3448 	       Thm ("divide_divide_eq_right",

  3449                      num_str @{thm divide_divide_eq_right}),

  3450 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3451 	       Thm ("divide_divide_eq_left",

  3452                      num_str @{thm divide_divide_eq_left}),

  3453 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3454 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")

  3455 	       ],

  3456       scr = EmptyScr

  3457       }:rls);

  3459 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3460 val reduce_0_1_2 = prep_rls(

  3461   Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

  3462       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3463       rules = [(*Thm ("divide_1",num_str @{thm divide_1}),

  3464 		 "?x / 1 = ?x" unnecess.for normalform*)

  3465 	       Thm ("mult_1_left",num_str @{thm mult_1_left}),

  3466 	       (*"1 * z = z"*)

  3467 	       (*Thm ("real_mult_minus1",num_str @{thm real_mult_minus1}),

  3468 	       "-1 * z = - z"*)

  3469 	       (*Thm ("real_minus_mult_cancel",num_str @{thm real_minus_mult_cancel}),

  3470 	       "- ?x * - ?y = ?x * ?y"*)

  3472 	       Thm ("mult_zero_left",num_str @{thm mult_zero_left}),

  3473 	       (*"0 * z = 0"*)

  3474 	       Thm ("add_0_left",num_str @{thm add_0_left}),

  3475 	       (*"0 + z = z"*)

  3476 	       (*Thm ("right_minus",num_str @{thm right_minus}),

  3477 	       "?z + - ?z = 0"*)

  3479 	       Thm ("sym_real_mult_2",

  3480                      num_str (@{thm real_mult_2} RS @{thm sym})),

  3481 	       (*"z1 + z1 = 2 * z1"*)

  3482 	       Thm ("real_mult_2_assoc",num_str @{thm real_mult_2_assoc}),

  3483 	       (*"z1 + (z1 + k) = 2 * z1 + k"*)

  3485 	       Thm ("divide_zero_left",num_str @{thm divide_zero_left})

  3486 	       (*"0 / ?x = 0"*)

  3487 	       ], scr = EmptyScr}:rls);

  3489 (*erls for calculate_Rational;

  3490   make local with FIXX@ME result:term *term list WN0609???SKMG*)

  3491 val norm_rat_erls = prep_rls(

  3492   Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3493       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3494       rules = [Calc ("Atools.is'_const",eval_const "#is_const_")

  3495 	       ], scr = EmptyScr}:rls);

  3497 (*.consists of rls containing the absolute minimum of thms.*)

  3498 (*040209: this version has been used by RL for his equations,

  3499 which is now replaced by MGs version below

  3500 vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)

  3501 val norm_Rational = prep_rls(

  3502   Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3503       erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3504       rules = [(*sequence given by operator precedence*)

  3505 	       Rls_ discard_minus,

  3506 	       Rls_ powers,

  3507 	       Rls_ rat_mult_divide,

  3508 	       Rls_ expand,

  3509 	       Rls_ reduce_0_1_2,

  3510 	       (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)

  3511 	       Rls_ order_add_mult,

  3512 	       Rls_ collect_numerals,

  3513 	       Rls_ add_fractions_p,

  3514 	       Rls_ cancel_p

  3515 	       ],

  3516       scr = EmptyScr}:rls);

  3518 val norm_Rational_parenthesized = prep_rls(

  3519   Seq {id = "norm_Rational_parenthesized", preconds = []:term list,

  3520        rew_ord = ("dummy_ord", dummy_ord),

  3521       erls = Atools_erls, srls = Erls,

  3522       calc = [], errpatts = [],

  3523       rules = [Rls_  norm_Rational, (*from RL -- not the latest one*)

  3524 	       Rls_ discard_parentheses

  3525 	       ],

  3526       scr = EmptyScr}:rls);

  3528 (*WN030318???SK: simplifies all but cancel and common_nominator*)

  3529 val simplify_rational =

  3530     merge_rls "simplify_rational" expand_binoms

  3531     (append_rls "divide" calculate_Rational

  3532 		[Thm ("divide_1",num_str @{thm divide_1}),

  3533 		 (*"?x / 1 = ?x"*)

  3534 		 Thm ("rat_mult",num_str @{thm rat_mult}),

  3535 		 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3536 		 Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

  3537 		 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3538 		 otherwise inv.to a / b / c = ...*)

  3539 		 Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

  3540 		 (*"?a / ?b * ?c = ?a * ?c / ?b"*)

  3541 		 Thm ("add_minus",num_str @{thm add_minus}),

  3542 		 (*"?a + ?b - ?b = ?a"*)

  3543 		 Thm ("add_minus1",num_str @{thm add_minus1}),

  3544 		 (*"?a - ?b + ?b = ?a"*)

  3545 		 Thm ("divide_minus1",num_str @{thm divide_minus1})

  3546 		 (*"?x / -1 = - ?x"*)

  3547 (*

  3548 ,

  3549 		 Thm ("",num_str @{thm })

  3550 *)

  3551 		 ]);

  3552 *}

  3553 ML {*

  3555 (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)

  3557 (* ------------------------------------------------------------------ *)

  3558 (*                  Simplifier für beliebige Buchterme                *)

  3559 (* ------------------------------------------------------------------ *)

  3560 (*----------------------- norm_Rational_mg ---------------------------*)

  3561 (*. description of the simplifier see MG-DA.p.56ff .*)

  3562 (* ------------------------------------------------------------------- *)

  3563 val common_nominator_p_rls = prep_rls(

  3564   Rls {id = "common_nominator_p_rls", preconds = [],

  3565        rew_ord = ("dummy_ord",dummy_ord),

  3566 	 erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3567 	 rules =

  3568 	 [Rls_ common_nominator_p

  3569 	  (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?

  3570 	    FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)

  3571 	  ],

  3572 	 scr = EmptyScr});

  3573 (* ------------------------------------------------------------------- *)

  3574 (* "Rls" causes repeated application of cancel_p to one and the same term *)

  3575 val cancel_p_rls = prep_rls(

  3576   Rls

  3577     {id = "cancel_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

  3578     erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3579     rules =

  3580       [Rls_ cancel_p (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)

  3581       ],

  3582 	  scr = EmptyScr});

  3583 (* -------------------------------------------------------------------- *)

  3584 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3585     used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)

  3586 val rat_mult_poly = prep_rls(

  3587   Rls {id = "rat_mult_poly", preconds = [],

  3588        rew_ord = ("dummy_ord",dummy_ord),

  3589 	 erls =  append_rls "e_rls-is_polyexp" e_rls

  3590 	         [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3591 	 srls = Erls, calc = [], errpatts = [],

  3592 	 rules =

  3593 	 [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

  3594 	  (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3595 	  Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})

  3596 	  (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3597 	  ],

  3598 	 scr = EmptyScr});

  3600 (* ------------------------------------------------------------------ *)

  3601 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3602     used in looping part norm_Rational_rls, see example DA-M02-main.p.60

  3603     .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls,

  3604     I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028

  3605     ... WN0609???MG.*)

  3606 val rat_mult_div_pow = prep_rls(

  3607   Rls {id = "rat_mult_div_pow", preconds = [],

  3608        rew_ord = ("dummy_ord",dummy_ord),

  3609        erls = e_rls,

  3610        (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls

  3611 			[Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3612          with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get

  3613 	 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.

  3614          thus we decided to go on with this flaw*)

  3615 		 srls = Erls, calc = [], errpatts = [],

  3616       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),

  3617 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3618 	       Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

  3619 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3620 	       Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}),

  3621 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3623 	       Thm ("real_divide_divide1_mg",

  3624                      num_str @{thm real_divide_divide1_mg}),

  3625 	       (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)

  3626 	       Thm ("divide_divide_eq_right",

  3627                      num_str @{thm divide_divide_eq_right}),

  3628 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3629 	       Thm ("divide_divide_eq_left",

  3630                      num_str @{thm divide_divide_eq_left}),

  3631 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3632 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e"),

  3634 	       Thm ("rat_power", num_str @{thm rat_power})

  3635 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

  3636 	       ],

  3637       scr = EmptyScr}:rls);

  3638 (* ------------------------------------------------------------------ *)

  3639 val rat_reduce_1 = prep_rls(

  3640   Rls {id = "rat_reduce_1", preconds = [],

  3641        rew_ord = ("dummy_ord",dummy_ord),

  3642        erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3643        rules = [Thm ("divide_1",num_str @{thm divide_1}),

  3644 		(*"?x / 1 = ?x"*)

  3645 		Thm ("mult_1_left",num_str @{thm mult_1_left})

  3646 		(*"1 * z = z"*)

  3647 		],

  3648        scr = EmptyScr}:rls);

  3649 (* ------------------------------------------------------------------ *)

  3650 (*. looping part of norm_Rational(*_mg*) .*)

  3651 val norm_Rational_rls = prep_rls(

  3652    Rls {id = "norm_Rational_rls", preconds = [],

  3653        rew_ord = ("dummy_ord",dummy_ord),

  3654        erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3655        rules = [Rls_ common_nominator_p_rls,

  3656 		Rls_ rat_mult_div_pow,

  3657 		Rls_ make_rat_poly_with_parentheses,

  3658 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)

  3659 		Rls_ rat_reduce_1

  3660 		],

  3661        scr = EmptyScr}:rls);

  3662 (* ------------------------------------------------------------------ *)

  3663 (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)

  3664  just be renaming:*)

  3665 val norm_Rational (*_mg*) = prep_rls(

  3666   Seq

  3667     {id = "norm_Rational"(*_mg*), preconds = [],

  3668     rew_ord = ("dummy_ord",dummy_ord),

  3669     erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3670     rules = [Rls_ discard_minus,

  3671       Rls_ rat_mult_poly,          (* removes double fractions like a/b/c    *)

  3672       Rls_ make_rat_poly_with_parentheses,           (*WN0510 also in(#)below*)

  3673       Rls_ cancel_p_rls,           (*FIXME.MG:cancel_p does NOT order sometim*)

  3674       Rls_ norm_Rational_rls,             (* the main rls, looping (#)       *)

  3675       Rls_ discard_parentheses1           (* mult only                       *)

  3676       ],

  3677     scr = EmptyScr}:rls);

  3678 "-------- rls norm_Rational --------------------------------------------------";

  3679 (* ------------------------------------------------------------------ *)

  3681 *}

  3682 ML {*

  3683 ruleset' := overwritelthy @{theory} (!ruleset',

  3684   [("calculate_Rational", calculate_Rational),

  3685    ("calc_rat_erls",calc_rat_erls),

  3686    ("rational_erls", rational_erls),

  3687    ("cancel_p", cancel_p),

  3688    ("common_nominator_p", common_nominator_p),

  3689    ("common_nominator_p_rls", common_nominator_p_rls),

  3690    (*WN120410 ("discard_minus", discard_minus), is defined in Poly.thy*)

  3691    ("powers_erls", powers_erls),

  3692    ("powers", powers),

  3693    ("rat_mult_divide", rat_mult_divide),

  3694    ("reduce_0_1_2", reduce_0_1_2),

  3695    ("rat_reduce_1", rat_reduce_1),

  3696    ("norm_rat_erls", norm_rat_erls),

  3697    ("norm_Rational", norm_Rational),

  3698    ("norm_Rational_rls", norm_Rational_rls),

  3699    ("norm_Rational_parenthesized", norm_Rational_parenthesized),

  3700    ("rat_mult_poly", rat_mult_poly),

  3701    ("rat_mult_div_pow", rat_mult_div_pow),

  3702    ("cancel_p_rls", cancel_p_rls)

  3703    ]);

  3705 calclist':= overwritel (!calclist',

  3706    [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))

  3707     ]);

  3709 (** problems **)

  3711 store_pbt

  3712  (prep_pbt thy "pbl_simp_rat" [] e_pblID

  3713  (["rational","simplification"],

  3714   [("#Given" ,["Term t_t"]),

  3715    ("#Where" ,["t_t is_ratpolyexp"]),

  3716    ("#Find"  ,["normalform n_n"])

  3717   ],

  3718   append_rls "e_rls" e_rls [(*for preds in where_*)],

  3719   SOME "Simplify t_t",

  3720   [["simplification","of_rationals"]]));

  3722 (** methods **)

  3724 (*WN061025 this methods script is copied from (auto-generated) script

  3725   of norm_Rational in order to ease repair on inform*)

  3726 store_met (prep_met thy "met_simp_rat" [] e_metID (["simplification","of_rationals"],

  3727 	  [("#Given" ,["Term t_t"]),

  3728 	   ("#Where" ,["t_t is_ratpolyexp"]),

  3729 	   ("#Find"  ,["normalform n_n"])],

  3730 	  {rew_ord'="tless_true", rls' = e_rls, calc = [], srls = e_rls,

  3731 	   prls = append_rls "simplification_of_rationals_prls" e_rls

  3732 				    [(*for preds in where_*) Calc ("Rational.is'_ratpolyexp", eval_is_ratpolyexp "")],

  3733 				   crls = e_rls, errpats = [],

  3734    nrls = norm_Rational_rls},

  3735    "Script SimplifyScript (t_t::real) =                             " ^

  3736    "  ((Try (Rewrite_Set discard_minus False) @@                   " ^

  3737    "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^

  3738    "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^

  3739    "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^

  3740    "    (Repeat                                                     " ^

  3741    "     ((Try (Rewrite_Set common_nominator_p_rls False) @@        " ^

  3742    "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^

  3743    "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^

  3744    "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^

  3745    "       Try (Rewrite_Set rat_reduce_1 False)))) @@               " ^

  3746    "    Try (Rewrite_Set discard_parentheses1 False))               " ^

  3747    "   t_t)"));

  3749 *}

  3750 ML {*"test"*}

  3751 end (* theory*)

author	Walther Neuper <neuper@ist.tugraz.at>
	Mon, 02 Sep 2013 15:17:34 +0200
changeset 52100	0831a4a6ec8a
parent 52096	ee2a5f066e44
child 52101	c3f399ce32af
permissions	-rwxr-xr-x