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

     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

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

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

   118   in

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

   120   end

   121 val is_expanded = is_poly

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

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

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

   126     [] @ term_of_es baseT expT vs es

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

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

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

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

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

   132     @ term_of_es baseT expT vs es

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

   135     let val es' = term_of_es baseT expT vs es

   136     in

   137       if c = 1

   138       then

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

   140         then Free (isastr_of_int c, baseT)

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

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

   143     end

   145 fun term_of_poly baseT expT vs p =

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

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

   148 *}

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

   151       (c) Stefan Karnel 2002

   152       Institute for Mathematics D and Institute for Software Technology,

   153       TU-Graz SS 2002 *}

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

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

   157     [2] 'polynomial' normalform (Polynom)

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

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

   160     Instead of

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

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

   163     we say:

   164       'fraction in normalform [2]'

   165       'fraction in normalform [3]'

   166     or

   167       'fraction [2]'

   168       'fraction [3]'.

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

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

   171 *}

   173 ML {*

   174 val thy = @{theory};

   176 signature RATIONALI =

   177 sig

   178   type mv_monom

   179   type mv_poly

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

   181   val calculate_Rational : rls

   182   val calc_rat_erls:rls

   183   val cancel_p : rls

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

   185   val check_fraction : term -> (term * term) option

   186   val check_frac_sum : term -> ((term * term) * (term * term)) option

   187   val common_nominator_p : rls

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

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

   190 			 (string * term) option

   191   val expanded2polynomial : term -> term option

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

   193   val is_expanded : term -> bool

   194   val is_polynomial : term -> bool

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

   197   val mv_lcm : mv_poly -> mv_poly -> mv_poly

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

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

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

   202                            roots into a simple and canceled fraction

   203                            with normalform [2].*)

   204 *)

   205 (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME

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

   207                            roots into a simple and canceled fraction

   208                            with normalform [2].*)

   209 *)

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

   211   val polynomial2expanded : term -> term option

   212   val rational_erls :

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

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

   216 end (* sig *)

   218 (*.**************************************************************************

   219 survey on the functions

   220 ~~~~~~~~~~~~~~~~~~~~~~~

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

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

   223  factout_p_         :                  | factout_          :

   224  cancel_p_          :                  | cancel_           :

   225                     :cancel_p          |                   :cancel

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

   227  common_nominator_p_:                  | common_nominator_ :

   228                     :common_nominator_p|                   :common_nominator

   229  add_fraction_p_    :                  | add_fraction_     :

   230 --------------------:------------------+-------------------:-----------------

   231 ???SK                 :norm_rational_p   |                   :norm_rational

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

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

   237 conversion from Isabelle-term to internal representation

   238 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   240 ... BITTE FORTSETZEN ...

   242 polynomial2expanded = ...

   243 expanded2polynomial = ...

   245 remark: polynomial2expanded o expanded2polynomial = I,

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

   248 functions for greatest common divisor and canceling

   249 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   250 ################################################################################

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

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

   253 ################################################################################

   254 mv_gcd

   255 factout_

   256 factout_p_

   257 cancel_

   258 cancel_p_

   260 functions for least common multiple and addition of fractions

   261 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   262 mv_lcm

   263 common_nominator_

   264 common_nominator_p_

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

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

   268 functions for normalform of rationals

   269 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

   271           schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?

   273 norm_rational_

   274 norm_expanded_rat_

   276 **************************************************************************.*)

   279 (*##*)

   280 structure RationalI : RATIONALI =

   281 struct

   282 (*##*)

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

   285 fun x mem [] = false

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

   290 val is_expanded = is_poly

   291 val is_polynomial = is_poly

   293 fun mk_noteq_0 baseT t =

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

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

   297 fun mk_asms baseT ts =

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

   299   in map (mk_noteq_0 baseT) as' end

   301 fun check_fraction t =

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

   303   in SOME (numerator, denominator) end

   304   handle Bind => NONE

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

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

   308 fun factout_p_ (thy: theory) t =

   309   let val opt = check_fraction t

   310   in

   311     case opt of

   312       NONE => NONE

   313     | SOME (numerator, denominator) =>

   314       let

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

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

   317         val baseT = type_of numerator

   318         val expT = HOLogic.realT

   319       in

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

   321           (SOME a, SOME b) =>

   322             let

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

   324               val es = replicate (length vs) 0

   325             in

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

   327               then NONE

   328               else

   329                 let

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

   331                   val ct = term_of_poly baseT expT vs c

   332                   val t' =

   333                     HOLogic.mk_binop "Fields.inverse_class.divide"

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

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

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

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

   338                 in SOME (t', asm) end

   339             end

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

   341       end

   342   end

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

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

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

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

   348   yields

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

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

   351     \<or> NONE *)

   352 fun cancel_p_ (_: theory) t =

   353   let val opt = check_fraction t

   354   in

   355     case opt of

   356       NONE => NONE

   357     | SOME (numerator, denominator) =>

   358       let

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

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

   361         val baseT = type_of numerator

   362         val expT = HOLogic.realT

   363       in

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

   365           (SOME a, SOME b) =>

   366             let

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

   368               val es = replicate (length vs) 0

   369             in

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

   371               then NONE

   372               else

   373                 let

   374                   val bt' = term_of_poly baseT expT vs b'

   375                   val ct = term_of_poly baseT expT vs c

   376                   val t' =

   377                     HOLogic.mk_binop "Fields.inverse_class.divide"

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

   379                   val asm = mk_asms baseT [bt']

   380                 in SOME (t', asm) end

   381             end

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

   383       end

   384   end

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

   387 fun check_frac_sum

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

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

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

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

   392   | check_frac_sum

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

   394       nofrac $

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

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

   397   | check_frac_sum

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

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

   400       nofrac)

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

   402   | check_frac_sum _ = NONE

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

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

   406 fun common_nominator_p_ (_: theory) t =

   407   let val opt = check_frac_sum t

   408   in

   409     case opt of

   410       NONE => NONE

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

   412       let

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

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

   415       in

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

   417           (SOME a, SOME b) =>

   418             let

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

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

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

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

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

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

   425               val denom =

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

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

   428                 else

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

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

   431               val t' =

   432                 HOLogic.mk_binop "Groups.plus_class.plus"

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

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

   435                   HOLogic.mk_binop "Fields.inverse_class.divide"

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

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

   438             in SOME (t', asm) end

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

   440       end

   441   end

   443 (* add fractions

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

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

   446 fun add_fraction_p_ (thy: theory) t =

   447   case check_frac_sum t of

   448     NONE => NONE

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

   450     let

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

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

   453     in

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

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

   456           let

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

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

   459             val nomin = term_of_poly baseT expT vs

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

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

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

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

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

   465     end

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

   468 fun xs union [] = xs

   469   | [] union ys = ys

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

   472 (*. gcd of integers .*)

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

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

   475 		  else gcd_int b (a mod b);

   477 (*. univariate polynomials (uv) .*)

   478 (*. univariate polynomials are represented as a list

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

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

   481 type uv_poly = int list;

   483 (*. adds two uv polynomials .*)

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

   485   | uv_mod_add_poly (p1,[]) = p1

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

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

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

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

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

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

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

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

   497 fun uv_mod_deg ([]:uv_poly) = 0

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

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

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

   502 fun uv_mod_mod2(x,p)=

   503     let

   504 	val y=(x mod p);

   505     in

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

   507 	    (

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

   509 	     )

   510     end;

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

   513 fun uv_mod_list_modp [] p = []

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

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

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

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

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

   521 (*. only for uv_mod_divides .*)

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

   523    integer zu klein  *)

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

   525     error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

   526   | uv_mod_pdiv p1 [x] =

   527     let

   528 	val xs= Unsynchronized.ref  [];

   529     in

   530 	if x<>0 then

   531 	    (

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

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

   534 	     )

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

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

   537     end

   538   | uv_mod_pdiv p1 p2 =

   539     let

   540 	val n= uv_mod_deg(p2);

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

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

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

   544 	val lc2=hd(p2');

   545 	val q= Unsynchronized.ref  [];

   546 	val c= Unsynchronized.ref  0;

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

   548     in

   549 	(

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

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

   552 	 else

   553 	     (

   554 	      if (!m)<n then

   555 		  (

   556 		   output:=([0],p1)

   557 		   )

   558 	      else

   559 		  (

   560 		   while (!m)>=n do

   561 		       (

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

   563 			if !c<>0 then

   564 			    (

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

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

   567 			     m:=uv_mod_deg(!p1')

   568 			     )

   569 			else m:=0

   570 			);

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

   572 		   )

   573 	      );

   574 	     !output

   575 	 )

   576     end;

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

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

   580     let

   581 	val n=uv_mod_deg(p2);

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

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

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

   585 	val lc2=hd(p2');

   586 	val q= Unsynchronized.ref  [];

   587 	val c= Unsynchronized.ref  0;

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

   589     in

   590 	(

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

   592 	 else

   593 	     (

   594 	      if (!m)<n then

   595 		  (

   596 		   output:=([0],p1)

   597 		   )

   598 	      else

   599 		  (

   600 		   while (!m)>=n do

   601 		       (

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

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

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

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

   606 			m:=(!m)-1

   607 			);

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

   610 		       (

   611 			p1':=tl(!p1')

   612 			);

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

   615 		   )

   616 	      );

   617 	     !output:uv_poly * uv_poly

   618 	 )

   619     end;

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

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

   623   | uv_mod_prest [] p2 = []:uv_poly

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

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

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

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

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

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

   632 fun uv_mod_cont ([]:uv_poly) = 0

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

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

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

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

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

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

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

   642   | uv_mod_pp p =

   643     let

   644 	val c= Unsynchronized.ref  0;

   645     in

   646 	(

   647 	 c:=uv_mod_cont(p);

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

   650 	 else uv_mod_div_list(p,!c)

   651 	)

   652     end;

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

   655 fun uv_mod_lc ([]:uv_poly) = 0

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

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

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

   660     let

   661 	val f = Unsynchronized.ref  [];

   662 	val f'= Unsynchronized.ref  p2;

   663 	val fi= Unsynchronized.ref  [];

   664     in

   665 	(

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

   667  	 while uv_mod_deg(!f')>0 do

   668 	     (

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

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

   671 		  (

   672 		   fi:=(!f');

   673 		   f:=(!fi)::(!f)

   674 		   )

   675 	      else ()

   676 	      );

   677 	      (!f)

   679 	 )

   680     end;

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

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

   684   | uv_mod_gcd_modp p1 [] p= p1

   685   | uv_mod_gcd_modp p1 p2 p=

   686     let

   687 	val p1'= Unsynchronized.ref [];

   688 	val p2'= Unsynchronized.ref [];

   689 	val pc= Unsynchronized.ref [];

   690 	val g= Unsynchronized.ref  [];

   691 	val d= Unsynchronized.ref  0;

   692 	val prs= Unsynchronized.ref  [];

   693     in

   694 	(

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

   696 	     (

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

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

   699 	      )

   700 	 else

   701 	     (

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

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

   704 	      );

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

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

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

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

   711 	 else pc:=hd(!prs);

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

   714 	 !g

   715 	 )

   716     end;

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

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

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

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

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

   725 fun uv_mod_add_qu [] = 0.0

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

   728 (*. calculates the euklidean norm .*)

   729 fun uv_mod_norm ([]:uv_poly) = 0.0

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

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

   733 fun uv_mod_multi a b = a * b;

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

   736 fun uv_mod_prim(x,[])= false

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

   738 		else false

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

   740 			then

   741 			    if uv_mod_prim(x,ys) then true

   742 			    else false

   743 		    else false;

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

   746 fun uv_mod_nextprime(g,p)=

   747     let

   748 	val list= Unsynchronized.ref  [2];

   749 	val exit= Unsynchronized.ref  0;

   750 	val i = Unsynchronized.ref 2

   751     in

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

   753 	    (

   754 	     if uv_mod_prim(!i,!list) then

   755 		 (

   756 		  if (p mod !i <> 0)

   757 		      then

   758 			  (

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

   760 			   i:= (!i)+1

   761 			   )

   762 		  else i:=(!i)+1

   763 		  )

   764 	     else i:= (!i)+1

   765 		 );

   766 	    i:=(p+1);

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

   768 	    (

   769 	     if uv_mod_prim(!i,!list) then

   770 		 (

   771 		  if (g mod !i <> 0)

   772 		      then

   773 			  (

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

   775 			   exit:= (!i)

   776 			   )

   777 		  else i:=(!i)+1

   778 		      )

   779 	     else i:= (!i)+1

   780 		 );

   781 	    !exit

   782     end;

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

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

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

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

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

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

   792 (*. chinese remainder algorithm .*)

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

   794     let

   795 	val c= Unsynchronized.ref  0;

   796 	val r1'= Unsynchronized.ref  0;

   797 	val d= Unsynchronized.ref  0;

   798 	val a= Unsynchronized.ref  0;

   799     in

   800 	(

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

   802 	     (

   803 	      c:=(!c)+1

   804 	      );

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

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

   807 	 !r1'+(!d)*m1

   808 	 )

   809     end;

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

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

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

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

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

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

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

   819     let

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

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

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

   823 	val temp= Unsynchronized.ref  [];

   824 	val cp= Unsynchronized.ref  [];

   825 	val qp= Unsynchronized.ref  [];

   826 	val q= Unsynchronized.ref [];

   827 	val pn= Unsynchronized.ref  0;

   828 	val d= Unsynchronized.ref  0;

   829 	val g1= Unsynchronized.ref  0;

   830 	val p= Unsynchronized.ref  0;

   831 	val m= Unsynchronized.ref  0;

   832 	val exit= Unsynchronized.ref  0;

   833 	val i= Unsynchronized.ref  1;

   834     in

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

   836 	else

   837 	    (

   838 	     temp:= !p1;

   839 	     p1:= !p2;

   840 	     p2:= !temp

   841 	     );

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

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

   846 	p:=4;

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

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

   850 	  Real.fromInt(!d) *

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

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

   854 	while (!exit=0) do

   855 	    (

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

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

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

   859 		 (

   860 		  i:=1;

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

   862 		      (

   863 		       i:=(!i)+1

   864 		       );

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

   866 		  )

   867 	     else ();

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

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

   873 	     pn:=(!p);

   874 	     q:=(!qp);

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

   877 		 (

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

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

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

   881  		      (

   882  		       i:=1;

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

   884  			   (

   885  			    i:=(!i)+1

   886 		           );

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

   888  		      )

   889  		  else ();

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

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

   893  		      (

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

   895  		       )

   896  		  else

   897  		      (

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

   899  			   (

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

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

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

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

   904 		 	    )

   905 		       else

   906 			   (

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

   908 				(

   909 				 pn:= !p;

   910 				 q:= !qp

   911 				 )

   912 			    else ()

   913 			    )

   914 		       )

   915 		  );

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

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

   918 	     );

   919 	    uv_mod_smul_poly(!q,c):uv_poly

   920     end;

   922 (*. multivariate polynomials .*)

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

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

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

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

   928 (*. global variables .*)

   929 (*. order indicators .*)

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

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

   933 (*. datatypes for internal representation.*)

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

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

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

   938 type mv_poly = mv_monom list;

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

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

   942 fun mv_mg_hlp([]) = EQUAL

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

   944 		    else if x>0 then GREATER

   945 			 else mv_mg_hlp(list);

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

   948 fun mv_addlist([]) = 0

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

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

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

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

   954     if order=LEX_ then

   955 	(

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

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

   958 	     )

   959     else

   960 	if order=GGO_ then

   961 	    (

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

   963 	     else

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

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

   966 	     )

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

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

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

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

   972 let

   973     val temp= Unsynchronized.ref  EQUAL;

   974 in

   975     if order=LEX_ then

   976 	(

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

   978 	     error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   979 	 else

   980 	     (

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

   982 	      if !temp=EQUAL then

   983 		  ( if x1=x2 then EQUAL

   984 		    else if x1>x2 then GREATER

   985 			 else LESS

   986 			     )

   987 	      else (!temp)

   988 	      )

   989 	     )

   990     else

   991 	if order=GGO_ then

   992 	    (

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

   994 	      error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   995 	     else

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

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

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

   999 		     )

  1000 	else error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")

  1001 end;

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

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

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

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

  1008 (*. leading power product .*)

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

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

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

  1013 (*. leading monomial .*)

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

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

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

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

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

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

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

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

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

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

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

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

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

  1031   | mv_lc(p1,order)  =

  1032     let

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

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

  1035 	val lc= Unsynchronized.ref  [];

  1036     in

  1037 	(

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

  1039 	     (

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

  1041 	      p1o:=tl(!p1o)

  1042 	      );

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

  1044 	 mv_rev_to(!lc)

  1045 	 )

  1046     end;

  1048 (*. compares two powerproducts .*)

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

  1051 (*. help function for mv_add .*)

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

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

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

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

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

  1057   | mv_madd(p1,p2,order) =

  1058     (

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

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

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

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

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

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

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

  1066 	      )

  1068 (*. adds two multivariate polynomials .*)

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

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

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

  1073 (*. monom multiplication .*)

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

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

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

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

  1080 (*. zeros a list .*)

  1081 fun mv_null2([])=[]

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

  1084 (*. multiplies two multivariate polynomials .*)

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

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

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

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

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

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

  1092 fun mv_getmax([])=0

  1093   | mv_getmax(x::p1)= let

  1094 		       val m=mv_getmax(p1);

  1095 		   in

  1096 		       if m>x then m

  1097 		       else x

  1098 		   end;

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

  1100 fun mv_grad([]:mv_poly) = 0

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

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

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

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

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

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

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

  1111 fun mv_is_negativ([])=false

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

  1114 (*. division of monomials .*)

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

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

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

  1118     let

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

  1120 	val pp= Unsynchronized.ref  [];

  1121     in

  1122 	(

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

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

  1125 	     if #1(p2)<>0 then

  1126 		 (

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

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

  1129 		  else (!c,!pp)

  1130 		      )

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

  1132 		 )

  1133     end;

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

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

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

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

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

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

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

  1144   | mv_division(f,g,order)=

  1145     let

  1146 	val r= Unsynchronized.ref  [];

  1147 	val q= Unsynchronized.ref  [];

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

  1149 	val k= Unsynchronized.ref  0;

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

  1151 	val exit= Unsynchronized.ref  0;

  1152     in

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

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

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

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

  1157 	else

  1158 	    (

  1159 	     exit:=0;

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

  1161 		 (

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

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

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

  1165 		      (

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

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

  1168 		       )

  1169 		  else exit:=1;

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

  1171 		  else (exit:=1)

  1172 		  );

  1173 		 (rev(!q),!r)

  1174 		 )

  1175     end;

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

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

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

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

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

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

  1185 (*. multivariate case .*)

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

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

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

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

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

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

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

  1196 fun mv_deg([]:mv_poly) = 0

  1197   | mv_deg(p1)=

  1198     let

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

  1200     in

  1201 	if length(p1')=0 then 0

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

  1203     end;

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

  1206 fun mv_deg2([]:mv_poly) = 0

  1207   | mv_deg2(p1)=

  1208     let

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

  1210     in

  1211 	if length(p1')=0 then 0

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

  1213     end;

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

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

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

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

  1220 fun uv_content2([]:mv_poly) = 0

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

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

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

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

  1226     let

  1227 	val count= Unsynchronized.ref  0;

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

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

  1230 	val list= Unsynchronized.ref  [];

  1231     in

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

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

  1234 	     else

  1235 		 (

  1236 		  count:=0;

  1237 		  while (!count)<=max do

  1238 		      (

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

  1240 			   (

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

  1242 			    help:=tl(!help)

  1243 			    )

  1244 		       else

  1245 			   (

  1246 			    list:= 0::(!list)

  1247 			    );

  1248 		       count := (!count) + 1

  1249 		       );

  1250 		      (!list)

  1251 		      )

  1252     end;

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

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

  1256   | uv_to_poly p1 =

  1257     let

  1258 	val count= Unsynchronized.ref  0;

  1259 	val help= Unsynchronized.ref  p1;

  1260 	val list= Unsynchronized.ref  [];

  1261     in

  1262 	while length(!help)>0 do

  1263 	    (

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

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

  1266 	     count:=(!count)+1;

  1267 	     help:=tl(!help)

  1268 	     );

  1269 	    (!list)

  1270     end;

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

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

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

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

  1276     let

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

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

  1279     in

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

  1281     end

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

  1283     let

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

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

  1286     in

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

  1288     end

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

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

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

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

  1294     let

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

  1296 	val t= Unsynchronized.ref  [];

  1297 	val y= Unsynchronized.ref  [];

  1298 	val n= Unsynchronized.ref  1;

  1299 	val n1= Unsynchronized.ref [];

  1300     in

  1301 	(

  1302 	 while length(!x)>1 do

  1303 	     (

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

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

  1306 		   else n1:=[];

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

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

  1309 	      x:=tl(!x)

  1310 	      );

  1311 	 (!y)

  1312 	 )

  1313     end;

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

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

  1317   | mv_newton_add [x:mv_poly] t = x

  1318   | mv_newton_add (pl:mv_poly list) t =

  1319     let

  1320 	val expos= Unsynchronized.ref  [];

  1321 	val pll= Unsynchronized.ref  pl;

  1322     in

  1323 	(

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

  1326 	     (

  1327 	      pll:=tl(!pll)

  1328 	      );

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

  1330 	 mv_add(hd(pl),

  1331 		mv_mul(

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

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

  1334 		       LEX_

  1335 		       ),

  1336 		LEX_)

  1337 	 )

  1338     end;

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

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

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

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

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

  1345   | mv_newton pl =

  1346     let

  1347 	val c= Unsynchronized.ref  pl;

  1348 	val c1= Unsynchronized.ref  [];

  1349 	val n=length(pl);

  1350 	val k= Unsynchronized.ref  1;

  1351 	val i= Unsynchronized.ref  n;

  1352 	val ppl= Unsynchronized.ref  [];

  1353     in

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

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

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

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

  1358 	    (

  1359 	     k:=(!k)+1;

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

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

  1362 	     ppl:= !c;

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

  1364 	     );

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

  1366 	c1:=rev(!c1);

  1367 	ppl:= !c1;

  1368 	mv_newton_add (!c1) 1

  1369     end;

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

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

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

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

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

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

  1378     let

  1379 	val y= Unsynchronized.ref  xs;

  1380 	val x= Unsynchronized.ref  [];

  1381     in

  1382 	(

  1383 	 x:=e;

  1384 	 while length(!y)>0 do

  1385 	     (

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

  1387 	      y:=tl(!y)

  1388 	      );

  1389 	 !x

  1390 	 )

  1391     end;

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

  1394 fun list_is_null [] = true

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

  1397 (* check if main variable is zero*)

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

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

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

  1402   | mv_content(p1) =

  1403     let

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

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

  1406 	val result= Unsynchronized.ref  [];

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

  1408     in

  1409 	(

  1410 	 if length(!list)>1 then

  1411 	     (

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

  1413 		  (

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

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

  1417 		   else list:=tl(!list)

  1419 		       );

  1420 		  if length(!list)>0 then

  1421 		   (

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

  1423 		    )

  1424 		  else list:=(!result);

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

  1426 		  (!list)

  1427 		  )

  1428 	 else

  1429 	     (

  1430 	      mv_null3(!list)

  1431 	      )

  1432 	     )

  1433     end

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

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

  1437   | mv_pp(p1) = let

  1438 		    val cont= Unsynchronized.ref  [];

  1439 		    val pp= Unsynchronized.ref [];

  1440 		in

  1441 		    cont:=mv_content(p1);

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

  1443 		    if !pp=[]

  1444 			then error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")

  1445 		    else (!pp)

  1446 		end

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

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

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

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

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

  1453      let

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

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

  1456      in

  1457 	(

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

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

  1460         )

  1461     end

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

  1463 	(

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

  1465 	)

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

  1467 	(

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

  1469         )

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

  1471     let

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

  1473 	val cont =

  1474 		  (

  1475                    if main_zero(mv_content(p1')) andalso

  1476                      (main_zero(mv_content(p2'))) then

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

  1478                    else

  1479                      mv_gcd (mv_content(p1')) (mv_content(p2'))

  1480                   );

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

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

  1483 	val gcd= Unsynchronized.ref  [];

  1484 	val candidate= Unsynchronized.ref  [];

  1485 	val interpolation_list= Unsynchronized.ref  [];

  1486 	val delta= Unsynchronized.ref  [];

  1487         val p1r = Unsynchronized.ref [];

  1488         val p2r = Unsynchronized.ref [];

  1489         val p1r' = Unsynchronized.ref [];

  1490         val p2r' = Unsynchronized.ref [];

  1491 	val factor= Unsynchronized.ref  [];

  1492 	val r= Unsynchronized.ref  0;

  1493 	val gcd_r= Unsynchronized.ref  [];

  1494 	val d= Unsynchronized.ref  0;

  1495 	val exit= Unsynchronized.ref  0;

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

  1497     in

  1498 	(

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

  1500 	     (

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

  1502 	      )

  1503 	 else

  1504 	     (

  1505 	      while !exit=0 do

  1506 		  (

  1507 		   r:=(!r)+1;

  1508                    p1r := mv_lc(p1,LEX_);

  1509 		   p2r := mv_lc(p2,LEX_);

  1510                    if main_zero(!p1r) andalso

  1511                       main_zero(!p2r)

  1512                    then

  1513                        (

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

  1515                        )

  1516                    else

  1517                        (

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

  1519                        );

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

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

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

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

  1524                    then

  1525                        (

  1526 		       )

  1527 		   else

  1528 		       (

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

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

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

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

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

  1534 			if (!d < !current_degree) then

  1535 			    (

  1536 			     current_degree:= !d;

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

  1538 			     )

  1539 			else

  1540 			    (

  1541 			     if (!d = !current_degree) then

  1542 				 (

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

  1544 				  )

  1545 			     else ()

  1546 				 )

  1547 			    );

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

  1549 			  (

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

  1551 			   if !candidate=[] then ()

  1552 			   else

  1553 			       (

  1554 				candidate:=mv_pp(!candidate);

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

  1556 				    (

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

  1558 				     exit:=1

  1559 				     )

  1560 				else ()

  1561 				    );

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

  1563 			       )

  1564 		      else ()

  1565 			  )

  1566 	     );

  1567 	     (!gcd):mv_poly

  1568 	     )

  1569     end;

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

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

  1574     (

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

  1576      );

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

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

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

  1583 fun count_neg ([]:mv_poly) = 0

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

  1585 			  else count_neg xs;

  1587 (*. help function for is_polynomial

  1588     checks the order of the operators .*)

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

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

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

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

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

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

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

  1596   | test_polynomial _ v = false;

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

  1599 fun is_polynomial t = test_polynomial t " ";

  1601 (*. help function for is_expanded

  1602     checks the order of the operators .*)

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

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

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

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

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

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

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

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

  1611   | test_exp  _ v = false;

  1614 (*. help function for check_coeff:

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

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

  1617     let

  1618 	val x= Unsynchronized.ref  NONE;

  1619 	val len= Unsynchronized.ref  0;

  1620 	val vl= Unsynchronized.ref  [];

  1621 	val vh= Unsynchronized.ref  [];

  1622 	val i= Unsynchronized.ref  0;

  1623     in

  1624 	if is_numeral str then

  1625 	    (

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

  1627 		 )

  1628 	else (* variable *)

  1629 	    (

  1630 	     len:=length(v);

  1631 	     vh:=v;

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

  1633 		 (

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

  1635 		      (

  1636 		       vl:=1::(!vl)

  1637 		       )

  1638 		  else

  1639 		      (

  1640 		       vl:=0::(!vl)

  1641 		       );

  1642 		      vh:=tl(!vh);

  1643 		      i:=(!i)+1

  1644 		      );

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

  1646 	    )

  1647     end

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

  1649     let

  1650 	val t1pp= Unsynchronized.ref  [];

  1651 	val t2pp= Unsynchronized.ref  [];

  1652 	val t1c= Unsynchronized.ref  0;

  1653 	val t2c= Unsynchronized.ref  0;

  1654     in

  1655 	(

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

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

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

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

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

  1663 	 )

  1664     end

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

  1666     let

  1667 	val x= Unsynchronized.ref  NONE;

  1668 	val len= Unsynchronized.ref  0;

  1669 	val vl= Unsynchronized.ref  [];

  1670 	val vh= Unsynchronized.ref  [];

  1671 	val vtemp= Unsynchronized.ref  [];

  1672 	val i= Unsynchronized.ref  0;

  1673     in

  1674     (

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

  1676 	 (

  1677 	  len:=length(v);

  1678 	  vh:=v;

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

  1681 	      (

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

  1683 		   (

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

  1685 		    )

  1686 	       else

  1687 		   (

  1688 		    vl:=0::(!vl)

  1689 		    );

  1690 		   vh:=tl(!vh);

  1691 		   i:=(!i)+1

  1692 		   );

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

  1694 	      )

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

  1696 	 )

  1697     end

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

  1699     (

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

  1701 	 )

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

  1703     (

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

  1705 	 )

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

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

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

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

  1711     else false;

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

  1714 fun mk_monom v' p vs =

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

  1716     in map (conv p) vs end;

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

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

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

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

  1722     if is_numeral str

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

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

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

  1727     let

  1728 	val x= Unsynchronized.ref  NONE;

  1729 	val len= Unsynchronized.ref  0;

  1730 	val vl= Unsynchronized.ref  [];

  1731 	val vh= Unsynchronized.ref  [];

  1732 	val i= Unsynchronized.ref  0;

  1733     in

  1734 	if is_numeral str then

  1735 	    (

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

  1737 		 )

  1738 	else (* variable *)

  1739 	    (

  1740 	     len:=length v;

  1741 	     vh:= v;

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

  1743 		 (

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

  1745 		      (

  1746 		       vl:=1::(!vl)

  1747 		       )

  1748 		  else

  1749 		      (

  1750 		       vl:=0::(!vl)

  1751 		       );

  1752 		      vh:=tl(!vh);

  1753 		      i:=(!i)+1

  1754 		      );

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

  1756 	    )

  1757     end

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

  1759     let

  1760 	val t1pp= Unsynchronized.ref  [];

  1761 	val t2pp= Unsynchronized.ref  [];

  1762 	val t1c= Unsynchronized.ref  0;

  1763 	val t2c= Unsynchronized.ref  0;

  1764     in

  1765 	(

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

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

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

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

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

  1772 	 handle _ => NONE

  1774 	 )

  1775     end

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

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

  1778     let

  1779 	val x= Unsynchronized.ref  NONE;

  1780 	val len= Unsynchronized.ref  0;

  1781 	val vl= Unsynchronized.ref  [];

  1782 	val vh= Unsynchronized.ref  [];

  1783 	val vtemp= Unsynchronized.ref  [];

  1784 	val i= Unsynchronized.ref  0;

  1785     in

  1786     (

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

  1788 	 (

  1789 	  len:=length(v);

  1790 	  vh:=v;

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

  1793 	      (

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

  1795 		   (

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

  1797 		    )

  1798 	       else

  1799 		   (

  1800 		    vl:=0::(!vl)

  1801 		    );

  1802 		   vh:=tl(!vh);

  1803 		   i:=(!i)+1

  1804 		   );

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

  1806 	      )

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

  1808 	 )

  1809     end

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

  1811     (

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

  1813 	 )

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

  1815     (

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

  1817 	 )

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

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

  1821     term2poly

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

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

  1824     			       DIE DU MIR LETZTES MAL GEGEBEN HAST*)

  1825     	 mv_monom list        (*internal representation           *)

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

  1827 .*)

  1828 fun term2poly (t:term) v =

  1829      if is_polynomial t then term2poly' t v

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

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

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

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

  1836 fun expanded2poly (t:term) v =

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

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

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

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

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

  1844 fun powerproduct2term(xs,v) =

  1845     let

  1846 	val xss= Unsynchronized.ref  xs;

  1847 	val vv= Unsynchronized.ref  v;

  1848     in

  1849 	(

  1850 	 while hd(!xss)=0 do

  1851 	     (

  1852 	      xss:=tl(!xss);

  1853 	      vv:=tl(!vv)

  1854 	      );

  1856 	 if list_is_null(tl(!xss)) then

  1857 	     (

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

  1859 	      else

  1860 		  (

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

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

  1863 		   )

  1864 	      )

  1865 	 else

  1866 	     (

  1867 	      if hd(!xss)=1 then

  1868 		  (

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

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

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

  1872 		   )

  1873 	      else

  1874 		  (

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

  1876 		   (

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

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

  1879 		    ) $

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

  1881 		   )

  1882 	      )

  1883 	 )

  1884     end;

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

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

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

  1889     else

  1890 	(

  1891 	 if list_is_null(e) then

  1892 	     (

  1893 	      Free(str_of_int c,HOLogic.realT)

  1894 	      )

  1895 	 else

  1896 	     (

  1897 	      if c=1 then

  1898 		  (

  1899 		   powerproduct2term(e,v)

  1900 		   )

  1901 	      else

  1902 		  (

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

  1904 		   Free(str_of_int c,HOLogic.realT)  $

  1905 		   powerproduct2term(e,v)

  1906 		   )

  1907 		  )

  1908 	     );*)

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

  1912     if list_is_null is

  1913     then

  1914 	if i >= 0

  1915 	then Free (str_of_int i, HOLogic.realT)

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

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

  1918     else

  1919 	if i > 0

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

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

  1922 		   powerproduct2term(is, v)

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

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

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

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

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

  1928     if list_is_null is

  1929     then Free (str_of_int i, HOLogic.realT)

  1930     else if i = 1

  1931     then powerproduct2term (is, vs)

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

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

  1934 	       powerproduct2term (is, vs);

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

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

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

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

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

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

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

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

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

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

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

  1949     else

  1950 	(

  1951 	 if list_is_null(e) then

  1952 	     (

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

  1954 	      )

  1955 	 else

  1956 	     (

  1957 	      if abs(c)=1 then

  1958 		  (

  1959 		   powerproduct2term(e,v)

  1960 		   )

  1961 	      else

  1962 		  (

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

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

  1965 		   powerproduct2term(e,v)

  1966 		   )

  1967 		  )

  1968 	     );

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

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

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

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

  1974     if c1<0 then

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

  1976 	exp2term'(others,vars) $

  1977 	(

  1978 	 monom2term2((c1,e1),vars)

  1979 	 )

  1980     else

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

  1982 	exp2term'(others,vars) $

  1983 	(

  1984 	 monom2term2((c1,e1),vars)

  1985 	 );

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

  1988     a term in polynomial representation .*)

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

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

  1992 fun polynomial2expanded t =

  1993     (let

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

  1995     in

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

  1997     end) handle _ => NONE;

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

  2000 fun expanded2polynomial t =

  2001     (let

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

  2003     in

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

  2005     end) handle _ => NONE;

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

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

  2010     let

  2011 	val p1' = Unsynchronized.ref [];

  2012 	val p2' = Unsynchronized.ref [];

  2013 	val p3  = Unsynchronized.ref []

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

  2015     in

  2016 	(

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

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

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

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

  2021 	     (

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

  2023 	      )

  2024 	 else

  2025 	     (

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

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

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

  2031 	      (

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

  2033 	       $

  2034 	       (

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

  2036 		poly2term(!p1',vars) $

  2037 		poly2term(!p3,vars)

  2038 		)

  2039 	       $

  2040 	       (

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

  2042 		poly2term(!p2',vars) $

  2043 		poly2term(!p3,vars)

  2044 		)

  2045 	       )

  2046 	      else

  2047 	      (

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

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

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

  2051 	       (

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

  2053 		$

  2054 		(

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

  2056 		 poly2term(!p1',vars) $

  2057 		 poly2term(!p3,vars)

  2058 		 )

  2059 		$

  2060 		(

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

  2062 		 poly2term(!p2',vars) $

  2063 		 poly2term(!p3,vars)

  2064 		 )

  2065 		)

  2066 	       )

  2067 	      )

  2068 	     )

  2069     end

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

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

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

  2074     let

  2075 	val p1' = Unsynchronized.ref [];

  2076 	val p2' = Unsynchronized.ref [];

  2077 	val p3  = Unsynchronized.ref []

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

  2079     in

  2080 	(

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

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

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

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

  2086 	     (

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

  2088 	      )

  2089 	 else

  2090 	     (

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

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

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

  2094 	      (

  2095 	       (

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

  2097 		$

  2098 		(

  2099 		 poly2term((!p1'),vars)

  2100 		 )

  2101 		$

  2102 		(

  2103 		 poly2term((!p2'),vars)

  2104 		 )

  2105 		)

  2106 	       ,

  2107 	       if mv_grad(!p3)>0 then

  2108 		   [

  2109 		    (

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

  2111 		     (

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

  2113 		      poly2term((!p3),vars) $

  2114 		      Free("0",HOLogic.realT)

  2115 		      )

  2116 		     )

  2117 		    ]

  2118 	       else

  2119 		   []

  2120 		   )

  2121 	      else

  2122 		  (

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

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

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

  2126 		       (

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

  2128 			$

  2129 			(

  2130 			 poly2term((!p1'),vars)

  2131 			 )

  2132 			$

  2133 			(

  2134 			 poly2term((!p2'),vars)

  2135 			 )

  2136 			,

  2137 			if mv_grad(!p3)>0 then

  2138 			    [

  2139 			     (

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

  2141 			      (

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

  2143 			       poly2term((!p3),vars) $

  2144 			       Free("0",HOLogic.realT)

  2145 			       )

  2146 			      )

  2147 			     ]

  2148 			else

  2149 			    []

  2150 			    )

  2151 		       )

  2152 		  )

  2153 	     )

  2154     end

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

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

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

  2159     let

  2160 	val p1' = Unsynchronized.ref [];

  2161 	val p2' = Unsynchronized.ref [];

  2162 	val p3  = Unsynchronized.ref []

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

  2164     in

  2165 	(

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

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

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

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

  2171 	     (

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

  2173 	      )

  2174 	 else

  2175 	     (

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

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

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

  2179 	      (

  2180 	       (

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

  2182 		$

  2183 		(

  2184 		 poly2expanded((!p1'),vars)

  2185 		 )

  2186 		$

  2187 		(

  2188 		 poly2expanded((!p2'),vars)

  2189 		 )

  2190 		)

  2191 	       ,

  2192 	       if mv_grad(!p3)>0 then

  2193 		   [

  2194 		    (

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

  2196 		     (

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

  2198 		      poly2expanded((!p3),vars) $

  2199 		      Free("0",HOLogic.realT)

  2200 		      )

  2201 		     )

  2202 		    ]

  2203 	       else

  2204 		   []

  2205 		   )

  2206 	      else

  2207 		  (

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

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

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

  2211 		       (

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

  2213 			$

  2214 			(

  2215 			 poly2expanded((!p1'),vars)

  2216 			 )

  2217 			$

  2218 			(

  2219 			 poly2expanded((!p2'),vars)

  2220 			 )

  2221 			,

  2222 			if mv_grad(!p3)>0 then

  2223 			    [

  2224 			     (

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

  2226 			      (

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

  2228 			       poly2expanded((!p3),vars) $

  2229 			       Free("0",HOLogic.realT)

  2230 			       )

  2231 			      )

  2232 			     ]

  2233 			else

  2234 			    []

  2235 			    )

  2236 		       )

  2237 		  )

  2238 	     )

  2239     end

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

  2243 (*. adds two fractions .*)

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

  2245     let

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

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

  2248 (* stopped Test_Isac.thy ...

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

  2250 *)

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

  2252 (* stopped Test_Isac.thy ...

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

  2254 *)

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

  2256 (* stopped Test_Isac.thy ...

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

  2258 *)

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

  2260 (* stopped Test_Isac.thy ...

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

  2262 *)

  2263 	val den= Unsynchronized.ref  [];

  2264 	val nom= Unsynchronized.ref  [];

  2265 	val m1= Unsynchronized.ref  [];

  2266 	val m2= Unsynchronized.ref  [];

  2267     in

  2269 	(

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

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

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

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

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

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

  2276 	 nom :=sort (mv_geq LEX_)

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

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

  2279 				       LEX_),

  2280 				LEX_));

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

  2282 	 (

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

  2284 	  $

  2285 	  (

  2286 	   poly2term((!nom),vars)

  2287 	   )

  2288 	  $

  2289 	  (

  2290 	   poly2term((!den),vars)

  2291 	   )

  2292 	  )

  2293 	 )

  2294     end

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

  2297 (*. adds two expanded fractions .*)

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

  2299     let

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

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

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

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

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

  2305 	val den= Unsynchronized.ref  [];

  2306 	val nom= Unsynchronized.ref  [];

  2307 	val m1= Unsynchronized.ref  [];

  2308 	val m2= Unsynchronized.ref  [];

  2309     in

  2311 	(

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

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

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

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

  2316 	 (

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

  2318 	  $

  2319 	  (

  2320 	   poly2expanded((!nom),vars)

  2321 	   )

  2322 	  $

  2323 	  (

  2324 	   poly2expanded((!den),vars)

  2325 	   )

  2326 	  )

  2327 	 )

  2328     end

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

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

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

  2333   | add_list_of_fractions [x]= direct_cancel x

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

  2335     let

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

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

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

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

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

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

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

  2343 val _= tracing"### add_list_of_fractions xs: has done direct_cancel(t3a)";

  2344 	val rest=rest1 union rest2 union rest3 union rest4;

  2345     in

  2346 	(tracing"### add_list_of_fractions in";

  2347 	 (

  2348 	 (t4a,rest)

  2349 	 )

  2350 	 )

  2351     end;

  2353 (*. adds a list of expanded terms .*)

  2354 fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[])

  2355   | add_list_of_fractions_exp [x]= direct_cancel_expanded x

  2356   | add_list_of_fractions_exp (x::y::xs) =

  2357     let

  2358 	val (t1a,rest1)=direct_cancel_expanded(x);

  2359 	val (t2a,rest2)=direct_cancel_expanded(y);

  2360 	val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs));

  2361 	val (t4a,rest4)=direct_cancel_expanded(t3a);

  2362 	val rest=rest1 union rest2 union rest3 union rest4;

  2363     in

  2364 	(

  2365 	 (t4a,rest)

  2366 	 )

  2367     end;

  2369 (*. calculates the lcm of a list of mv_poly .*)

  2370 fun calc_lcm ([x],var)= (x,var)

  2371   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2373 (*. converts a list of terms to a list of mv_poly .*)

  2374 fun t2d([],_)=[]

  2375   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2377 (*. same as t2d, this time for expanded forms .*)

  2378 fun t2d_exp([],_)=[]

  2379   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2381 (*. converts a list of fract terms to a list of their denominators .*)

  2382 fun termlist2denominators [] = ([],[])

  2383   | termlist2denominators xs =

  2384     let

  2385 	val xxs= Unsynchronized.ref  xs;

  2386 	val var= Unsynchronized.ref  [];

  2387     in

  2388 	var:=[];

  2389 	while length(!xxs)>0 do

  2390 	    (

  2391 	     let

  2392 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2393 	     in

  2394 		 (

  2395 		  xxs:=tl(!xxs);

  2396 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2397 		  )

  2398 	     end

  2399 	     );

  2400 	    (t2d(xs,!var),!var)

  2401     end;

  2403 (*. calculates the lcm of a list of mv_poly .*)

  2404 fun calc_lcm ([x],var)= (x,var)

  2405   | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var);

  2407 (*. converts a list of terms to a list of mv_poly .*)

  2408 fun t2d([],_)=[]

  2409   | t2d((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars);

  2411 (*. same as t2d, this time for expanded forms .*)

  2412 fun t2d_exp([],_)=[]

  2413   | t2d_exp((t as (Const("Fields.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars);

  2415 (*. converts a list of fract terms to a list of their denominators .*)

  2416 fun termlist2denominators [] = ([],[])

  2417   | termlist2denominators xs =

  2418     let

  2419 	val xxs= Unsynchronized.ref  xs;

  2420 	val var= Unsynchronized.ref  [];

  2421     in

  2422 	var:=[];

  2423 	while length(!xxs)>0 do

  2424 	    (

  2425 	     let

  2426 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2427 	     in

  2428 		 (

  2429 		  xxs:=tl(!xxs);

  2430 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2431 		  )

  2432 	     end

  2433 	     );

  2434 	    (t2d(xs,!var),!var)

  2435     end;

  2437 (*. same as termlist2denminators, this time for expanded forms .*)

  2438 fun termlist2denominators_exp [] = ([],[])

  2439   | termlist2denominators_exp xs =

  2440     let

  2441 	val xxs= Unsynchronized.ref  xs;

  2442 	val var= Unsynchronized.ref  [];

  2443     in

  2444 	var:=[];

  2445 	while length(!xxs)>0 do

  2446 	    (

  2447 	     let

  2448 		 val (t as Const ("Fields.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs);

  2449 	     in

  2450 		 (

  2451 		  xxs:=tl(!xxs);

  2452 		  var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var))

  2453 		  )

  2454 	     end

  2455 	     );

  2456 	    (t2d_exp(xs,!var),!var)

  2457     end;

  2459 (*. reduces all fractions to the least common denominator .*)

  2460 fun com_den(x::xs,denom,den,var)=

  2461     let

  2462 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2463 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2464 	val p3= #1(mv_division(denom,p2,LEX_));

  2465 	val p1var=get_vars(p1');

  2466     in

  2467 	if length(xs)>0 then

  2468 	    if p3=[(1,mv_null2(var))] then

  2469 		(

  2470 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2471 		 $

  2472 		 (

  2473 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2474 		  $

  2475 		  poly2term(the (term2poly p1' p1var),p1var)

  2476 		  $

  2477 		  den

  2478 		  )

  2479 		 $

  2480 		 #1(com_den(xs,denom,den,var))

  2481 		,

  2482 		[]

  2483 		)

  2484 	    else

  2485 		(

  2486 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2487 		 $

  2488 		 (

  2489 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2490 		  $

  2491 		  (

  2492 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2493 		   poly2term(the (term2poly p1' p1var),p1var) $

  2494 		   poly2term(p3,var)

  2495 		   )

  2496 		  $

  2497 		  (

  2498 		   den

  2499 		   )

  2500 		  )

  2501 		 $

  2502 		 #1(com_den(xs,denom,den,var))

  2503 		,

  2504 		[]

  2505 		)

  2506 	else

  2507 	    if p3=[(1,mv_null2(var))] then

  2508 		(

  2509 		 (

  2510 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2511 		  $

  2512 		  poly2term(the (term2poly p1' p1var),p1var)

  2513 		  $

  2514 		  den

  2515 		  )

  2516 		 ,

  2517 		 []

  2518 		 )

  2519 	     else

  2520 		 (

  2521 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2522 		  $

  2523 		  (

  2524 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2525 		   poly2term(the (term2poly p1' p1var),p1var) $

  2526 		   poly2term(p3,var)

  2527 		   )

  2528 		  $

  2529 		  den

  2530 		  ,

  2531 		  []

  2532 		  )

  2533     end;

  2535 (*. same as com_den, this time for expanded forms .*)

  2536 fun com_den_exp(x::xs,denom,den,var)=

  2537     let

  2538 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2539 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2540 	val p3= #1(mv_division(denom,p2,LEX_));

  2541 	val p1var=get_vars(p1');

  2542     in

  2543 	if length(xs)>0 then

  2544 	    if p3=[(1,mv_null2(var))] then

  2545 		(

  2546 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2547 		 $

  2548 		 (

  2549 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2550 		  $

  2551 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2552 		  $

  2553 		  den

  2554 		  )

  2555 		 $

  2556 		 #1(com_den_exp(xs,denom,den,var))

  2557 		,

  2558 		[]

  2559 		)

  2560 	    else

  2561 		(

  2562 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2563 		 $

  2564 		 (

  2565 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2566 		  $

  2567 		  (

  2568 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2569 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2570 		   poly2expanded(p3,var)

  2571 		   )

  2572 		  $

  2573 		  (

  2574 		   den

  2575 		   )

  2576 		  )

  2577 		 $

  2578 		 #1(com_den_exp(xs,denom,den,var))

  2579 		,

  2580 		[]

  2581 		)

  2582 	else

  2583 	    if p3=[(1,mv_null2(var))] then

  2584 		(

  2585 		 (

  2586 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2587 		  $

  2588 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2589 		  $

  2590 		  den

  2591 		  )

  2592 		 ,

  2593 		 []

  2594 		 )

  2595 	     else

  2596 		 (

  2597 		  Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2598 		  $

  2599 		  (

  2600 		   Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2601 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2602 		   poly2expanded(p3,var)

  2603 		   )

  2604 		  $

  2605 		  den

  2606 		  ,

  2607 		  []

  2608 		  )

  2609     end;

  2611 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2612 -------------------------------------------------------------

  2613 (* WN0210???SK brauch ma des überhaupt *)

  2614 fun com_den2(x::xs,denom,den,var)=

  2615     let

  2616 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2617 	val p2= sort (mv_geq LEX_) (the(term2poly p2' var));

  2618 	val p3= #1(mv_division(denom,p2,LEX_));

  2619 	val p1var=get_vars(p1');

  2620     in

  2621 	if length(xs)>0 then

  2622 	    if p3=[(1,mv_null2(var))] then

  2623 		(

  2624 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2625 		 poly2term(the(term2poly p1' p1var),p1var) $

  2626 		 com_den2(xs,denom,den,var)

  2627 		)

  2628 	    else

  2629 		(

  2630 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2631 		 (

  2632 		   let

  2633 		       val p3'=poly2term(p3,var);

  2634 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2635 		   in

  2636 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2637 		   end

  2638 		  ) $

  2639 		 com_den2(xs,denom,den,var)

  2640 		)

  2641 	else

  2642 	    if p3=[(1,mv_null2(var))] then

  2643 		(

  2644 		 poly2term(the(term2poly p1' p1var),p1var)

  2645 		 )

  2646 	     else

  2647 		 (

  2648 		   let

  2649 		       val p3'=poly2term(p3,var);

  2650 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2651 		   in

  2652 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2653 		   end

  2654 		  )

  2655     end;

  2657 (* WN0210???SK brauch ma des überhaupt *)

  2658 fun com_den_exp2(x::xs,denom,den,var)=

  2659     let

  2660 	val (t as Const ("Fields.inverse_class.divide",_) $ p1' $ p2')=x;

  2661 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

  2662 	val p3= #1(mv_division(denom,p2,LEX_));

  2663 	val p1var=get_vars p1';

  2664     in

  2665 	if length(xs)>0 then

  2666 	    if p3=[(1,mv_null2(var))] then

  2667 		(

  2668 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2669 		 poly2expanded(the (expanded2poly p1' p1var),p1var) $

  2670 		 com_den_exp2(xs,denom,den,var)

  2671 		)

  2672 	    else

  2673 		(

  2674 		 Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2675 		 (

  2676 		   let

  2677 		       val p3'=poly2expanded(p3,var);

  2678 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2679 		   in

  2680 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2681 		   end

  2682 		  ) $

  2683 		 com_den_exp2(xs,denom,den,var)

  2684 		)

  2685 	else

  2686 	    if p3=[(1,mv_null2(var))] then

  2687 		(

  2688 		 poly2expanded(the (expanded2poly p1' p1var),p1var)

  2689 		 )

  2690 	     else

  2691 		 (

  2692 		   let

  2693 		       val p3'=poly2expanded(p3,var);

  2694 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2695 		   in

  2696 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2697 		   end

  2698 		  )

  2699     end;

  2700 ---------------------------------------------------------*)

  2703 (*. searches for an element y of a list ys, which has an gcd not 1 with x .*)

  2704 fun exists_gcd (x,[]) = false

  2705   | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then  exists_gcd (x,ys)

  2706 			   else true;

  2708 (*. divides each element of the list xs with y .*)

  2709 fun list_div ([],y) = []

  2710   | list_div (x::xs,y) =

  2711     let

  2712 	val (d,r)=mv_division(x,y,LEX_);

  2713     in

  2714 	if r=[] then

  2715 	    d::list_div(xs,y)

  2716 	else x::list_div(xs,y)

  2717     end;

  2719 (*. checks if x is in the list ys .*)

  2720 fun in_list (x,[]) = false

  2721   | in_list (x,y::ys) = if x=y then true

  2722 			else in_list(x,ys);

  2724 (*. deletes all equal elements of the list xs .*)

  2725 fun kill_equal [] = []

  2726   | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)

  2727 			 else x::kill_equal(xs);

  2729 (*. searches for new factors .*)

  2730 fun new_factors [] = []

  2731   | new_factors (list:mv_poly list):mv_poly list =

  2732     let

  2733 	val l = kill_equal list;

  2734 	val len = length(l);

  2735     in

  2736 	if len>=2 then

  2737 	    (

  2738 	     let

  2739 		 val x::y::xs=l;

  2740 		 val gcd=mv_gcd x y;

  2741 	     in

  2742 		 if gcd=[(1,mv_null2(#2(hd(x))))] then

  2743 		     (

  2744 		      if exists_gcd(x,xs) then new_factors (y::xs @ [x])

  2745 		      else x::new_factors(y::xs)

  2746 	             )

  2747 		 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))

  2748 	     end

  2749 	     )

  2750 	else

  2751 	    if len=1 then [hd(l)]

  2752 	    else []

  2753     end;

  2755 (*. gets the factors of a list .*)

  2756 fun get_factors x = new_factors x;

  2758 (*. multiplies the elements of the list .*)

  2759 fun multi_list [] = []

  2760   | multi_list (x::xs) = if xs=[] then x

  2761 			 else mv_mul(x,multi_list xs,LEX_);

  2763 (*. makes a term out of the elements of the list (polynomial representation) .*)

  2764 fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2765   | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)

  2766 			       else

  2767 				   (

  2768 				    Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2769 				    poly2term(sort (mv_geq LEX_) (x),vars) $

  2770 				    make_term(xs,vars)

  2771 				    );

  2773 (*. factorizes the denominator (polynomial representation) .*)

  2774 fun factorize_den (l,den,vars) =

  2775     let

  2776 	val factor_list=kill_equal( (get_factors l));

  2777 	val mlist=multi_list(factor_list);

  2778 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2779     in

  2780 	if rest=[] then

  2781 	    (

  2782 	     if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)

  2783 	     else make_term(last::factor_list,vars)

  2784 	     )

  2785 	else error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")

  2786     end;

  2788 (*. makes a term out of the elements of the list (expanded polynomial representation) .*)

  2789 fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2790   | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)

  2791 			       else

  2792 				   (

  2793 				    Const ("Groups.times_class.times",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2794 				    poly2expanded(sort (mv_geq LEX_) (x),vars) $

  2795 				    make_exp(xs,vars)

  2796 				    );

  2798 (*. factorizes the denominator (expanded polynomial representation) .*)

  2799 fun factorize_den_exp (l,den,vars) =

  2800     let

  2801 	val factor_list=kill_equal( (get_factors l));

  2802 	val mlist=multi_list(factor_list);

  2803 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2804     in

  2805 	if rest=[] then

  2806 	    (

  2807 	     if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)

  2808 	     else make_exp(last::factor_list,vars)

  2809 	     )

  2810 	else error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")

  2811     end;

  2813 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2814 (*. the nominators and denominators with the correct factor .*)

  2815 (*. (polynomial representation) .*)

  2816 fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)

  2817   | step_add_list_of_fractions [x]= error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")

  2818   | step_add_list_of_fractions (xs) =

  2819     let

  2820         val den_list=termlist2denominators (xs); (* list of denominators *)

  2821 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2822 	val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2823     in

  2824 	com_den(xs,denom,den,var)

  2825     end;

  2827 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2828 (*. the nominators and denominators with the correct factor .*)

  2829 (*. (expanded polynomial representation) .*)

  2830 fun step_add_list_of_fractions_exp []  = (Free("0",HOLogic.realT),[]:term list)

  2831   | step_add_list_of_fractions_exp [x] = error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")

  2832   | step_add_list_of_fractions_exp (xs)=

  2833     let

  2834         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2835 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2836 	val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2837     in

  2838 	com_den_exp(xs,denom,den,var)

  2839     end;

  2841 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2842 -------------------------------------------------------------

  2843 (* WN0210???SK brauch ma des überhaupt *)

  2844 fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)

  2845   | step_add_list_of_fractions2 [x]=(x,[])

  2846   | step_add_list_of_fractions2 (xs) =

  2847     let

  2848         val den_list=termlist2denominators (xs); (* list of denominators *)

  2849 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2850 	val den=factorize_den(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2851     in

  2852 	(

  2853 	 Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2854 	 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $

  2855 	 poly2term(denom,var)

  2856 	,

  2857 	[]

  2858 	)

  2859     end;

  2861 (* WN0210???SK brauch ma des überhaupt *)

  2862 fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)

  2863   | step_add_list_of_fractions2_exp [x]=(x,[])

  2864   | step_add_list_of_fractions2_exp (xs) =

  2865     let

  2866         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2867 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2868 	val den=factorize_den_exp(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2869     in

  2870 	(

  2871 	 Const ("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2872 	 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $

  2873 	 poly2expanded(denom,var)

  2874 	,

  2875 	[]

  2876 	)

  2877     end;

  2878 ---------------------------------------------- *)

  2881 (* converts a term, which contains several terms seperated by +, into a list of these terms .*)

  2882 fun term2list (t as (Const("Fields.inverse_class.divide",_) $ _ $ _)) = [t]

  2883   | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =

  2884       [Const ("Fields.inverse_class.divide",

  2885         [HOLogic.realT,HOLogic.realT] ---> HOLogic.realT) $

  2886 	  t $ Free("1",HOLogic.realT)

  2887      ]

  2888   | term2list (t as (Free(_,_))) =

  2889     [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2890 	  t $  Free("1",HOLogic.realT)

  2891      ]

  2892   | term2list (t as (Const("Groups.times_class.times",_) $ _ $ _)) =

  2893     [Const("Fields.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2894 	  t $ Free("1",HOLogic.realT)

  2895      ]

  2896   | term2list (Const("Groups.plus_class.plus",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)

  2897   | term2list (Const("Groups.minus_class.minus",_) $ t1 $ t2) =

  2898     error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")

  2899   | term2list _ = error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");

  2901 (*.factors out the gcd of nominator and denominator:

  2902    a/b = (a' * gcd)/(b' * gcd),  a,b,gcd  are poly[2].*)

  2904 (*. brings the term into a normal form .*)

  2905 fun norm_rational_ (thy:theory) t =

  2906     SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;

  2907 fun norm_expanded_rat_ (thy:theory) t =

  2908     SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;

  2911 (*.evaluates conditions in calculate_Rational.*)

  2912 (*make local with FIXX@ME result:term *term list*)

  2913 val calc_rat_erls = prep_rls(

  2914   Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  2915 	 erls = e_rls, srls = Erls, calc = [], errpatts = [],

  2916 	 rules =

  2917 	 [Calc ("HOL.eq",eval_equal "#equal_"),

  2918 	  Calc ("Atools.is'_const",eval_const "#is_const_"),

  2919 	  Thm ("not_true",num_str @{thm not_true}),

  2920 	  Thm ("not_false",num_str @{thm not_false})

  2921 	  ],

  2922 	 scr = EmptyScr});

  2925 (*.simplifies expressions with numerals;

  2926    does NOT rearrange the term by AC-rewriting; thus terms with variables

  2927    need to have constants to be commuted together respectively.*)

  2928 val calculate_Rational = prep_rls (merge_rls "calculate_Rational"

  2929 	  (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  2930 	    erls = calc_rat_erls, srls = Erls,

  2931 	    calc = [], errpatts = [],

  2932 	    rules =

  2933 	      [Calc ("Fields.inverse_class.divide",eval_cancel "#divide_e"),

  2935 	       Thm ("minus_divide_left",num_str (@{thm minus_divide_left} RS @{thm sym})),

  2936 	         (*SYM - ?x / ?y = - (?x / ?y)  may come from subst*)

  2938 	       Thm ("rat_add",num_str @{thm rat_add}),

  2939 	         (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \

  2940 		           \a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*)

  2941 	       Thm ("rat_add1",num_str @{thm rat_add1}),

  2942 	         (*"[| a is_const; b is_const; c is_const |] ==> a / c + b / c = (a + b) / c"*)

  2943 	       Thm ("rat_add2",num_str @{thm rat_add2}),

  2944 	         (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> ?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)

  2945 	       Thm ("rat_add3",num_str @{thm rat_add3}),

  2946 	         (*"[| a is_const; b is_const; c is_const |] ==> a + b / c = (a * c) / c + b / c"\

  2947 		           .... is_const to be omitted here FIXME*)

  2949 	       Thm ("rat_mult",num_str @{thm rat_mult}),

  2950 	         (*a / b * (c / d) = a * c / (b * d)*)

  2951 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

  2952 	         (*?x * (?y / ?z) = ?x * ?y / ?z*)

  2953 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

  2954 	         (*?y / ?z * ?x = ?y * ?x / ?z*)

  2956 	       Thm ("real_divide_divide1",num_str @{thm real_divide_divide1}),

  2957 	         (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)

  2958 	       Thm ("divide_divide_eq_left",num_str @{thm divide_divide_eq_left}),

  2959 	         (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  2961 	       Thm ("rat_power", num_str @{thm rat_power}),

  2962 	         (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

  2964 	       Thm ("mult_cross",num_str @{thm mult_cross}),

  2965 	         (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*)

  2966 	       Thm ("mult_cross1",num_str @{thm mult_cross1}),

  2967 	         (*"   b ~= 0            ==> (a / b = c    ) = (a     = b * c)*)

  2968 	       Thm ("mult_cross2",num_str @{thm mult_cross2})

  2969 	         (*"           d ~= 0    ==> (a     = c / d) = (a * d =     c)*)

  2970 	       ], scr = EmptyScr})

  2971 	  calculate_Poly);

  2973 (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)

  2974 fun eval_is_expanded (thmid:string) _

  2975 		       (t as (Const("Rational.is'_expanded", _) $ arg)) thy =

  2976     if is_expanded arg

  2977     then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  2978 	         Trueprop $ (mk_equality (t, @{term True})))

  2979     else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  2980 	         Trueprop $ (mk_equality (t, @{term False})))

  2981   | eval_is_expanded _ _ _ _ = NONE;

  2983 val rational_erls =

  2984     merge_rls "rational_erls" calculate_Rational

  2985 	      (append_rls "is_expanded" Atools_erls

  2986 			  [Calc ("Rational.is'_expanded", eval_is_expanded "")

  2987 			   ]);

  2990 (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:

  2991  =================================================================

  2992  A[2] 'cancel_p': .

  2993  A[3] 'cancel': .

  2994  B[2] 'common_nominator_p': transforms summands in a term [2]

  2995          to fractions with the (least) common multiple as nominator.

  2996  B[3] 'norm_rational': normalizes arbitrary algebraic terms (without

  2997          radicals and transzendental functions) to one canceled fraction,

  2998 	 nominator and denominator in polynomial form.

  3000 In order to meet isac's requirements for interactive and stepwise calculation,

  3001 each 'reverse-rewerite-set' consists of an initialization for the interpreter

  3002 state and of 4 functions, each of which employs rewriting as much as possible.

  3003 The signature of these functions are the same in each 'reverse-rewrite-set'

  3004 respectively.*)

  3006 (* ************************************************************************* *)

  3008 local(*. cancel_p

  3009 ------------------------

  3010 cancels a single fraction consisting of two (uni- or multivariate)

  3011 polynomials WN0609???SK[2] into another such a fraction; examples:

  3013 	   a^2 + -1*b^2         a + b

  3014         -------------------- = ---------

  3015 	a^2 + -2*a*b + b^2     a + -1*b

  3017         a^2    a

  3018         --- = ---

  3019          a     1

  3021 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)

  3022 (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)

  3024 val {rules, rew_ord=(_,ro),...} =

  3025     rep_rls (assoc_rls "make_polynomial");

  3026 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3027            see rational.sml --- investigate rulesets for cancel_p ---*)

  3028 val {rules, rew_ord=(_,ro),...} =

  3029     rep_rls (assoc_rls "rev_rew_p");

  3031 (*.init_state = fn : term -> istate

  3032 initialzies the state of the script interpreter. The state is:

  3034 type rrlsstate =      (*state for reverse rewriting*)

  3035      (term *          (*the current formula*)

  3036       term *          (*the final term*)

  3037       rule list       (*'reverse rule list' (#)*)

  3038 	    list *    (*may be serveral, eg. in norm_rational*)

  3039       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3040        (term *        (*... rewrite with ...*)

  3041 	term list))   (*... assumptions*)

  3042 	  list);      (*derivation from given term to normalform

  3043 		       in reverse order with sym_thm;

  3044                        (#) could be extracted from here by (map #1)*).*)

  3045 (* val {rules, rew_ord=(_,ro),...} =

  3046        rep_rls (assoc_rls "rev_rew_p")        (*USE ALWAYS, SEE val cancel_p*);

  3047    val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*);

  3048    val t = t;

  3049    *)

  3050 fun init_state thy eval_rls ro t =

  3051     let val SOME (t',_) = factout_p_ thy t

  3052         val SOME (t'',asm) = cancel_p_ thy t

  3053         val der = reverse_deriv thy eval_rls rules ro NONE t'

  3054         val der = der @ [(Thm ("real_mult_div_cancel2",

  3055 			       num_str @{thm real_mult_div_cancel2}),

  3056 			  (t'',asm))]

  3057         val rs = (distinct_Thm o (map #1)) der

  3058 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3059 				      "sym_real_mult_0",

  3060 				      "sym_real_mult_1"

  3061 				      (*..insufficient,eg.make_Polynomial*)])rs

  3062     in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end;

  3064 (*.locate_rule = fn : rule list -> term -> rule

  3065 		      -> (rule * (term * term list) option) list.

  3066   checks a rule R for being a cancel-rule, and if it is,

  3067   then return the list of rules (+ the terms they are rewriting to)

  3068   which need to be applied before R should be applied.

  3069   precondition: the rule is applicable to the argument-term.

  3070 arguments:

  3071   rule list: the reverse rule list

  3072   -> term  : ... to which the rule shall be applied

  3073   -> rule  : ... to be applied to term

  3074 value:

  3075   -> (rule           : a rule rewriting to ...

  3076       * (term        : ... the resulting term ...

  3077          * term list): ... with the assumptions ( //#0).

  3078       ) list         : there may be several such rules;

  3079 		       the list is empty, if the rule has nothing to do

  3080 		       with cancelation.*)

  3081 (* val () = ();

  3082    *)

  3083 fun locate_rule thy eval_rls ro [rs] t r =

  3084     if (id_of_thm r) mem (map (id_of_thm)) rs

  3085     then let val ropt =

  3086 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3087 	 in case ropt of

  3088 		SOME ta => [(r, ta)]

  3089 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3090 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3091 			 []) end

  3092     else (tracing("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  3093   | locate_rule _ _ _ _ _ _ =

  3094     error ("locate_rule: doesnt match rev-sets in istate");

  3096 (*.next_rule = fn : rule list -> term -> rule option

  3097   for a given term return the next rules to be done for cancelling.

  3098 arguments:

  3099   rule list     : the reverse rule list

  3100   term          : the term for which ...

  3101 value:

  3102   -> rule option: ... this rule is appropriate for cancellation;

  3103 		  there may be no such rule (if the term is canceled already.*)

  3104 (* val thy = thy;

  3105    val Rrls {rew_ord=(_,ro),...} = cancel;

  3106    val ([rs],t) = (rss,f);

  3107    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3109    val (thy, [rs]) = (thy, revsets);

  3110    val Rrls {rew_ord=(_,ro),...} = cancel;

  3111    nex [rs] t;

  3112    *)

  3113 fun next_rule thy eval_rls ro [rs] t =

  3114     let val der = make_deriv thy eval_rls rs ro NONE t;

  3115     in case der of

  3116 (* val (_,r,_)::_ = der;

  3117    *)

  3118 	   (_,r,_)::_ => SOME r

  3119 	 | _ => NONE

  3120     end

  3121   | next_rule _ _ _ _ _ =

  3122     error ("next_rule: doesnt match rev-sets in istate");

  3124 (*.val attach_form = f : rule list -> term -> term

  3125 			 -> (rule * (term * term list)) list

  3126   checks an input term TI, if it may belong to a current cancellation, by

  3127   trying to derive it from the given term TG.

  3128 arguments:

  3129   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3130   -> term: TI, the next one input by the user

  3131 value:

  3132   -> (rule           : the rule to be applied in order to reach TI

  3133       * (term        : ... obtained by applying the rule ...

  3134          * term list): ... and the respective assumptions.

  3135       ) list         : there may be several such rules;

  3136                        the list is empty, if the users term does not belong

  3137 		       to a cancellation of the term last agreed upon.*)

  3138 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  3139     []:(rule * (term * term list)) list;

  3141 in

  3143 val cancel_p =

  3144     Rrls {id = "cancel_p", prepat=[],

  3145 	  rew_ord=("ord_make_polynomial",

  3146 		   ord_make_polynomial false thy),

  3147 	  erls = rational_erls,

  3148 	  calc = [("PLUS"    ,("Groups.plus_class.plus"        ,eval_binop "#add_")),

  3149 		  ("TIMES"   ,("Groups.times_class.times"        ,eval_binop "#mult_")),

  3150 		  ("DIVIDE" ,("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")),

  3151 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3152 	  errpatts = [],

  3153 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3154 		     normal_form = cancel_p_ thy,

  3155 		     locate_rule = locate_rule thy Atools_erls ro,

  3156 		     next_rule   = next_rule thy Atools_erls ro,

  3157 		     attach_form = attach_form}}

  3158 end;(*local*)

  3160 local(*.ad [2] 'common_nominator_p'

  3161 ---------------------------------

  3162 FIXME Beschreibung .*)

  3165 val {rules=rules,rew_ord=(_,ro),...} =

  3166     rep_rls (assoc_rls "make_polynomial");

  3167 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3168            see rational.sml --- investigate rulesets for cancel_p ---*)

  3169 val {rules, rew_ord=(_,ro),...} =

  3170     rep_rls (assoc_rls "rev_rew_p");

  3171 val thy = thy;

  3174 (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option

  3175   as defined above*)

  3177 (*.init_state = fn : term -> istate

  3178 initialzies the state of the interactive interpreter. The state is:

  3180 type rrlsstate =      (*state for reverse rewriting*)

  3181      (term *          (*the current formula*)

  3182       term *          (*the final term*)

  3183       rule list       (*'reverse rule list' (#)*)

  3184 	    list *    (*may be serveral, eg. in norm_rational*)

  3185       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3186        (term *        (*... rewrite with ...*)

  3187 	term list))   (*... assumptions*)

  3188 	  list);      (*derivation from given term to normalform

  3189 		       in reverse order with sym_thm;

  3190                        (#) could be extracted from here by (map #1)*).*)

  3191 fun init_state thy eval_rls ro t =

  3192     let val SOME (t',_) = common_nominator_p_ thy t;

  3193         val SOME (t'',asm) = add_fraction_p_ thy t;

  3194         val der = reverse_deriv thy eval_rls rules ro NONE t';

  3195         val der = der @ [(Thm ("real_mult_div_cancel2",

  3196 			       num_str @{thm real_mult_div_cancel2}),

  3197 			  (t'',asm))]

  3198         val rs = (distinct_Thm o (map #1)) der;

  3199 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3200 				      "sym_real_mult_0",

  3201 				      "sym_real_mult_1"]) rs;

  3202     in (t,t'',[rs(*here only _ONE_*)],der) end;

  3204 (* use"knowledge/Rational.ML";

  3205    *)

  3207 (*.locate_rule = fn : rule list -> term -> rule

  3208 		      -> (rule * (term * term list) option) list.

  3209   checks a rule R for being a cancel-rule, and if it is,

  3210   then return the list of rules (+ the terms they are rewriting to)

  3211   which need to be applied before R should be applied.

  3212   precondition: the rule is applicable to the argument-term.

  3213 arguments:

  3214   rule list: the reverse rule list

  3215   -> term  : ... to which the rule shall be applied

  3216   -> rule  : ... to be applied to term

  3217 value:

  3218   -> (rule           : a rule rewriting to ...

  3219       * (term        : ... the resulting term ...

  3220          * term list): ... with the assumptions ( //#0).

  3221       ) list         : there may be several such rules;

  3222 		       the list is empty, if the rule has nothing to do

  3223 		       with cancelation.*)

  3224 (* val () = ();

  3225    *)

  3226 fun locate_rule thy eval_rls ro [rs] t r =

  3227     if (id_of_thm r) mem (map (id_of_thm)) rs

  3228     then let val ropt =

  3229 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3230 	 in case ropt of

  3231 		SOME ta => [(r, ta)]

  3232 	      | NONE => (tracing("### locate_rule:  rewrite "^

  3233 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3234 			 []) end

  3235     else (tracing("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  3236   | locate_rule _ _ _ _ _ _ =

  3237     error ("locate_rule: doesnt match rev-sets in istate");

  3239 (*.next_rule = fn : rule list -> term -> rule option

  3240   for a given term return the next rules to be done for cancelling.

  3241 arguments:

  3242   rule list     : the reverse rule list

  3243   term          : the term for which ...

  3244 value:

  3245   -> rule option: ... this rule is appropriate for cancellation;

  3246 		  there may be no such rule (if the term is canceled already.*)

  3247 (* val thy = thy;

  3248    val Rrls {rew_ord=(_,ro),...} = cancel;

  3249    val ([rs],t) = (rss,f);

  3250    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3252    val (thy, [rs]) = (thy, revsets);

  3253    val Rrls {rew_ord=(_,ro),...} = cancel;

  3254    nex [rs] t;

  3255    *)

  3256 fun next_rule thy eval_rls ro [rs] t =

  3257     let val der = make_deriv thy eval_rls rs ro NONE t;

  3258     in case der of

  3259 (* val (_,r,_)::_ = der;

  3260    *)

  3261 	   (_,r,_)::_ => SOME r

  3262 	 | _ => NONE

  3263     end

  3264   | next_rule _ _ _ _ _ =

  3265     error ("next_rule: doesnt match rev-sets in istate");

  3267 (*.val attach_form = f : rule list -> term -> term

  3268 			 -> (rule * (term * term list)) list

  3269   checks an input term TI, if it may belong to a current cancellation, by

  3270   trying to derive it from the given term TG.

  3271 arguments:

  3272   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3273   -> term: TI, the next one input by the user

  3274 value:

  3275   -> (rule           : the rule to be applied in order to reach TI

  3276       * (term        : ... obtained by applying the rule ...

  3277          * term list): ... and the respective assumptions.

  3278       ) list         : there may be several such rules;

  3279                        the list is empty, if the users term does not belong

  3280 		       to a cancellation of the term last agreed upon.*)

  3281 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  3282     []:(rule * (term * term list)) list;

  3284 (* if each pat* matches with the current term, the Rrls is applied

  3285    (there are no preconditions to be checked, they are @{term True}) *)

  3286 val pat0 = parse_patt thy "?r/?s+?u/?v :: real";

  3287 val pat1 = parse_patt thy "?r/?s+?u    :: real";

  3288 val pat2 = parse_patt thy "?r   +?u/?v :: real";

  3289 val prepat = [([@{term True}], pat0),

  3290 	      ([@{term True}], pat1),

  3291 	      ([@{term True}], pat2)];

  3292 in

  3294 (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;

  3295   besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht

  3296   dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)

  3297 val common_nominator_p =

  3298     Rrls {id = "common_nominator_p", prepat=prepat,

  3299 	  rew_ord=("ord_make_polynomial",

  3300 		   ord_make_polynomial false thy),

  3301 	  erls = rational_erls,

  3302 	  calc = [("PLUS"    ,("Groups.plus_class.plus"        ,eval_binop "#add_")),

  3303 		  ("TIMES"   ,("Groups.times_class.times"        ,eval_binop "#mult_")),

  3304 		  ("DIVIDE" ,("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")),

  3305 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3306 	  errpatts = [],

  3307 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3308 		     normal_form = add_fraction_p_ thy, (*FIXME.WN0211*)

  3309 		     locate_rule = locate_rule thy Atools_erls ro,

  3310 		     next_rule   = next_rule thy Atools_erls ro,

  3311 		     attach_form = attach_form}}

  3312 end;(*local*)

  3314 end; (*struct*)

  3316 *}

  3317 ML {*

  3318 open RationalI;

  3319 (*##*)

  3321 (*.the expression contains + - * ^ / only ?.*)

  3322 fun is_ratpolyexp (Free _) = true

  3323   | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ Free _ $ Free _) = true

  3324   | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ Free _ $ Free _) = true

  3325   | is_ratpolyexp (Const ("Groups.times_class.times",_) $ Free _ $ Free _) = true

  3326   | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true

  3327   | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ Free _ $ Free _) = true

  3328   | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ t1 $ t2) =

  3329                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3330   | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) =

  3331                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3332   | is_ratpolyexp (Const ("Groups.times_class.times",_) $ t1 $ t2) =

  3333                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3334   | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =

  3335                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3336   | is_ratpolyexp (Const ("Fields.inverse_class.divide",_) $ t1 $ t2) =

  3337                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3338   | is_ratpolyexp _ = false;

  3340 (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)

  3341 fun eval_is_ratpolyexp (thmid:string) _

  3342 		       (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =

  3343     if is_ratpolyexp arg

  3344     then SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  3345 	         Trueprop $ (mk_equality (t, @{term True})))

  3346     else SOME (mk_thmid thmid "" (term_to_string''' thy arg) "",

  3347 	         Trueprop $ (mk_equality (t, @{term False})))

  3348   | eval_is_ratpolyexp _ _ _ _ = NONE;

  3350 (*("get_denominator", ("Rational.get_denominator", eval_get_denominator ""))*)

  3351 fun eval_get_denominator (thmid:string) _

  3352 		      (t as Const ("Rational.get_denominator", _) $

  3353               (Const ("Fields.inverse_class.divide", _) $ num $

  3354                 denom)) thy =

  3355       SOME (mk_thmid thmid "" (term_to_string''' thy denom) "",

  3356 	            Trueprop $ (mk_equality (t, denom)))

  3357   | eval_get_denominator _ _ _ _ = NONE;

  3359 (*("get_numerator", ("Rational.get_numerator", eval_get_numerator ""))*)

  3360 fun eval_get_numerator (thmid:string) _

  3361       (t as Const ("Rational.get_numerator", _) $

  3362           (Const ("Fields.inverse_class.divide", _) $num

  3363             $denom )) thy =

  3364     SOME (mk_thmid thmid "" (term_to_string''' thy num) "",

  3365 	    Trueprop $ (mk_equality (t, num)))

  3366   | eval_get_numerator _ _ _ _ = NONE;

  3368 (*-------------------18.3.03 --> struct <-----------vvv--*)

  3369 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)

  3371 (*WN100906 removed in favour of discard_minus = discard_minus_ formerly

  3372    discard binary minus, shift unary minus into -1*;

  3373    unary minus before numerals are put into the numeral by parsing;

  3374    contains absolute minimum of thms for context in norm_Rational

  3375 *** val discard_minus = prep_rls(

  3376   Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3377       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3378       rules =

  3379         [Thm ("real_diff_minus", num_str @{thm real_diff_minus}),

  3380 	           (*"a - b = a + -1 * b"*)

  3381 	         Thm ("sym_real_mult_minus1", num_str (@{thm real_mult_minus1} RS @{thm sym}))

  3382 	           (*- ?z = "-1 * ?z"*)],

  3383       scr = EmptyScr

  3384       }):rls;

  3385  *)

  3387 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)

  3388 val powers_erls = prep_rls(

  3389   Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3390       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3391       rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),

  3392 	       Calc ("Atools.is'_even",eval_is_even "#is_even_"),

  3393 	       Calc ("Orderings.ord_class.less",eval_equ "#less_"),

  3394 	       Thm ("not_false", num_str @{thm not_false}),

  3395 	       Thm ("not_true", num_str @{thm not_true}),

  3396 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

  3397 	       ],

  3398       scr = EmptyScr

  3399       }:rls);

  3400 (*.all powers over + distributed; atoms over * collected, other distributed

  3401    contains absolute minimum of thms for context in norm_Rational .*)

  3402 val powers = prep_rls(

  3403   Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3404       erls = powers_erls, srls = Erls, calc = [], errpatts = [],

  3405       rules = [Thm ("realpow_multI", num_str @{thm realpow_multI}),

  3406 	       (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)

  3407 	       Thm ("realpow_pow",num_str @{thm realpow_pow}),

  3408 	       (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)

  3409 	       Thm ("realpow_oneI",num_str @{thm realpow_oneI}),

  3410 	       (*"r ^^^ 1 = r"*)

  3411 	       Thm ("realpow_minus_even",num_str @{thm realpow_minus_even}),

  3412 	       (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)

  3413 	       Thm ("realpow_minus_odd",num_str @{thm realpow_minus_odd}),

  3414 	       (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)

  3416 	       (*----- collect atoms over * -----*)

  3417 	       Thm ("realpow_two_atom",num_str @{thm realpow_two_atom}),

  3418 	       (*"r is_atom ==> r * r = r ^^^ 2"*)

  3419 	       Thm ("realpow_plus_1",num_str @{thm realpow_plus_1}),

  3420 	       (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)

  3421 	       Thm ("realpow_addI_atom",num_str @{thm realpow_addI_atom}),

  3422 	       (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)

  3424 	       (*----- distribute none-atoms -----*)

  3425 	       Thm ("realpow_def_atom",num_str @{thm realpow_def_atom}),

  3426 	       (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)

  3427 	       Thm ("realpow_eq_oneI",num_str @{thm realpow_eq_oneI}),

  3428 	       (*"1 ^^^ n = 1"*)

  3429 	       Calc ("Groups.plus_class.plus",eval_binop "#add_")

  3430 	       ],

  3431       scr = EmptyScr

  3432       }:rls);

  3433 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3434 val rat_mult_divide = prep_rls(

  3435   Rls {id = "rat_mult_divide", preconds = [],

  3436        rew_ord = ("dummy_ord",dummy_ord),

  3437       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3438       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),

  3439 	       (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3440 	       Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

  3441 	       (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3442 	       otherwise inv.to a / b / c = ...*)

  3443 	       Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

  3444 	       (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much

  3445 		     and does not commute a / b * c ^^^ 2 !*)

  3447 	       Thm ("divide_divide_eq_right",

  3448                      num_str @{thm divide_divide_eq_right}),

  3449 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3450 	       Thm ("divide_divide_eq_left",

  3451                      num_str @{thm divide_divide_eq_left}),

  3452 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3453 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")

  3454 	       ],

  3455       scr = EmptyScr

  3456       }:rls);

  3458 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3459 val reduce_0_1_2 = prep_rls(

  3460   Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

  3461       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3462       rules = [(*Thm ("divide_1",num_str @{thm divide_1}),

  3463 		 "?x / 1 = ?x" unnecess.for normalform*)

  3464 	       Thm ("mult_1_left",num_str @{thm mult_1_left}),

  3465 	       (*"1 * z = z"*)

  3466 	       (*Thm ("real_mult_minus1",num_str @{thm real_mult_minus1}),

  3467 	       "-1 * z = - z"*)

  3468 	       (*Thm ("real_minus_mult_cancel",num_str @{thm real_minus_mult_cancel}),

  3469 	       "- ?x * - ?y = ?x * ?y"*)

  3471 	       Thm ("mult_zero_left",num_str @{thm mult_zero_left}),

  3472 	       (*"0 * z = 0"*)

  3473 	       Thm ("add_0_left",num_str @{thm add_0_left}),

  3474 	       (*"0 + z = z"*)

  3475 	       (*Thm ("right_minus",num_str @{thm right_minus}),

  3476 	       "?z + - ?z = 0"*)

  3478 	       Thm ("sym_real_mult_2",

  3479                      num_str (@{thm real_mult_2} RS @{thm sym})),

  3480 	       (*"z1 + z1 = 2 * z1"*)

  3481 	       Thm ("real_mult_2_assoc",num_str @{thm real_mult_2_assoc}),

  3482 	       (*"z1 + (z1 + k) = 2 * z1 + k"*)

  3484 	       Thm ("divide_zero_left",num_str @{thm divide_zero_left})

  3485 	       (*"0 / ?x = 0"*)

  3486 	       ], scr = EmptyScr}:rls);

  3488 (*erls for calculate_Rational;

  3489   make local with FIXX@ME result:term *term list WN0609???SKMG*)

  3490 val norm_rat_erls = prep_rls(

  3491   Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3492       erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3493       rules = [Calc ("Atools.is'_const",eval_const "#is_const_")

  3494 	       ], scr = EmptyScr}:rls);

  3496 (*.consists of rls containing the absolute minimum of thms.*)

  3497 (*040209: this version has been used by RL for his equations,

  3498 which is now replaced by MGs version below

  3499 vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)

  3500 val norm_Rational = prep_rls(

  3501   Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3502       erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3503       rules = [(*sequence given by operator precedence*)

  3504 	       Rls_ discard_minus,

  3505 	       Rls_ powers,

  3506 	       Rls_ rat_mult_divide,

  3507 	       Rls_ expand,

  3508 	       Rls_ reduce_0_1_2,

  3509 	       (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)

  3510 	       Rls_ order_add_mult,

  3511 	       Rls_ collect_numerals,

  3512 	       Rls_ add_fractions_p,

  3513 	       Rls_ cancel_p

  3514 	       ],

  3515       scr = EmptyScr}:rls);

  3517 val norm_Rational_parenthesized = prep_rls(

  3518   Seq {id = "norm_Rational_parenthesized", preconds = []:term list,

  3519        rew_ord = ("dummy_ord", dummy_ord),

  3520       erls = Atools_erls, srls = Erls,

  3521       calc = [], errpatts = [],

  3522       rules = [Rls_  norm_Rational, (*from RL -- not the latest one*)

  3523 	       Rls_ discard_parentheses

  3524 	       ],

  3525       scr = EmptyScr}:rls);

  3527 (*WN030318???SK: simplifies all but cancel and common_nominator*)

  3528 val simplify_rational =

  3529     merge_rls "simplify_rational" expand_binoms

  3530     (append_rls "divide" calculate_Rational

  3531 		[Thm ("divide_1",num_str @{thm divide_1}),

  3532 		 (*"?x / 1 = ?x"*)

  3533 		 Thm ("rat_mult",num_str @{thm rat_mult}),

  3534 		 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3535 		 Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}),

  3536 		 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3537 		 otherwise inv.to a / b / c = ...*)

  3538 		 Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}),

  3539 		 (*"?a / ?b * ?c = ?a * ?c / ?b"*)

  3540 		 Thm ("add_minus",num_str @{thm add_minus}),

  3541 		 (*"?a + ?b - ?b = ?a"*)

  3542 		 Thm ("add_minus1",num_str @{thm add_minus1}),

  3543 		 (*"?a - ?b + ?b = ?a"*)

  3544 		 Thm ("divide_minus1",num_str @{thm divide_minus1})

  3545 		 (*"?x / -1 = - ?x"*)

  3546 (*

  3547 ,

  3548 		 Thm ("",num_str @{thm })

  3549 *)

  3550 		 ]);

  3551 *}

  3552 ML {*

  3554 (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)

  3556 (* ------------------------------------------------------------------ *)

  3557 (*                  Simplifier für beliebige Buchterme                *)

  3558 (* ------------------------------------------------------------------ *)

  3559 (*----------------------- norm_Rational_mg ---------------------------*)

  3560 (*. description of the simplifier see MG-DA.p.56ff .*)

  3561 (* ------------------------------------------------------------------- *)

  3562 val common_nominator_p_rls = prep_rls(

  3563   Rls {id = "common_nominator_p_rls", preconds = [],

  3564        rew_ord = ("dummy_ord",dummy_ord),

  3565 	 erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3566 	 rules =

  3567 	 [Rls_ common_nominator_p

  3568 	  (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?

  3569 	    FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)

  3570 	  ],

  3571 	 scr = EmptyScr});

  3572 (* ------------------------------------------------------------------- *)

  3573 (* "Rls" causes repeated application of cancel_p to one and the same term *)

  3574 val cancel_p_rls = prep_rls(

  3575   Rls

  3576     {id = "cancel_p_rls", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

  3577     erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3578     rules =

  3579       [Rls_ cancel_p (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)

  3580       ],

  3581 	  scr = EmptyScr});

  3582 (* -------------------------------------------------------------------- *)

  3583 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3584     used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)

  3585 val rat_mult_poly = prep_rls(

  3586   Rls {id = "rat_mult_poly", preconds = [],

  3587        rew_ord = ("dummy_ord",dummy_ord),

  3588 	 erls =  append_rls "e_rls-is_polyexp" e_rls

  3589 	         [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3590 	 srls = Erls, calc = [], errpatts = [],

  3591 	 rules =

  3592 	 [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

  3593 	  (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3594 	  Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r})

  3595 	  (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3596 	  ],

  3597 	 scr = EmptyScr});

  3599 (* ------------------------------------------------------------------ *)

  3600 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3601     used in looping part norm_Rational_rls, see example DA-M02-main.p.60

  3602     .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls,

  3603     I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028

  3604     ... WN0609???MG.*)

  3605 val rat_mult_div_pow = prep_rls(

  3606   Rls {id = "rat_mult_div_pow", preconds = [],

  3607        rew_ord = ("dummy_ord",dummy_ord),

  3608        erls = e_rls,

  3609        (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls

  3610 			[Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3611          with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get

  3612 	 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.

  3613          thus we decided to go on with this flaw*)

  3614 		 srls = Erls, calc = [], errpatts = [],

  3615       rules = [Thm ("rat_mult",num_str @{thm rat_mult}),

  3616 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3617 	       Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

  3618 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3619 	       Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}),

  3620 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3622 	       Thm ("real_divide_divide1_mg",

  3623                      num_str @{thm real_divide_divide1_mg}),

  3624 	       (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)

  3625 	       Thm ("divide_divide_eq_right",

  3626                      num_str @{thm divide_divide_eq_right}),

  3627 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3628 	       Thm ("divide_divide_eq_left",

  3629                      num_str @{thm divide_divide_eq_left}),

  3630 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3631 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e"),

  3633 	       Thm ("rat_power", num_str @{thm rat_power})

  3634 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

  3635 	       ],

  3636       scr = EmptyScr}:rls);

  3637 (* ------------------------------------------------------------------ *)

  3638 val rat_reduce_1 = prep_rls(

  3639   Rls {id = "rat_reduce_1", preconds = [],

  3640        rew_ord = ("dummy_ord",dummy_ord),

  3641        erls = e_rls, srls = Erls, calc = [], errpatts = [],

  3642        rules = [Thm ("divide_1",num_str @{thm divide_1}),

  3643 		(*"?x / 1 = ?x"*)

  3644 		Thm ("mult_1_left",num_str @{thm mult_1_left})

  3645 		(*"1 * z = z"*)

  3646 		],

  3647        scr = EmptyScr}:rls);

  3648 (* ------------------------------------------------------------------ *)

  3649 (*. looping part of norm_Rational(*_mg*) .*)

  3650 val norm_Rational_rls = prep_rls(

  3651    Rls {id = "norm_Rational_rls", preconds = [],

  3652        rew_ord = ("dummy_ord",dummy_ord),

  3653        erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3654        rules = [Rls_ common_nominator_p_rls,

  3655 		Rls_ rat_mult_div_pow,

  3656 		Rls_ make_rat_poly_with_parentheses,

  3657 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)

  3658 		Rls_ rat_reduce_1

  3659 		],

  3660        scr = EmptyScr}:rls);

  3661 (* ------------------------------------------------------------------ *)

  3662 (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)

  3663  just be renaming:*)

  3664 val norm_Rational (*_mg*) = prep_rls(

  3665   Seq

  3666     {id = "norm_Rational"(*_mg*), preconds = [],

  3667     rew_ord = ("dummy_ord",dummy_ord),

  3668     erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

  3669     rules = [Rls_ discard_minus,

  3670       Rls_ rat_mult_poly,          (* removes double fractions like a/b/c    *)

  3671       Rls_ make_rat_poly_with_parentheses,           (*WN0510 also in(#)below*)

  3672       Rls_ cancel_p_rls,           (*FIXME.MG:cancel_p does NOT order sometim*)

  3673       Rls_ norm_Rational_rls,             (* the main rls, looping (#)       *)

  3674       Rls_ discard_parentheses1           (* mult only                       *)

  3675       ],

  3676     scr = EmptyScr}:rls);

  3677 "-------- rls norm_Rational --------------------------------------------------";

  3678 (* ------------------------------------------------------------------ *)

  3680 *}

  3681 ML {*

  3682 ruleset' := overwritelthy @{theory} (!ruleset',

  3683   [("calculate_Rational", calculate_Rational),

  3684    ("calc_rat_erls",calc_rat_erls),

  3685    ("rational_erls", rational_erls),

  3686    ("cancel_p", cancel_p),

  3687    ("common_nominator_p", common_nominator_p),

  3688    ("common_nominator_p_rls", common_nominator_p_rls),

  3689    (*WN120410 ("discard_minus", discard_minus), is defined in Poly.thy*)

  3690    ("powers_erls", powers_erls),

  3691    ("powers", powers),

  3692    ("rat_mult_divide", rat_mult_divide),

  3693    ("reduce_0_1_2", reduce_0_1_2),

  3694    ("rat_reduce_1", rat_reduce_1),

  3695    ("norm_rat_erls", norm_rat_erls),

  3696    ("norm_Rational", norm_Rational),

  3697    ("norm_Rational_rls", norm_Rational_rls),

  3698    ("norm_Rational_parenthesized", norm_Rational_parenthesized),

  3699    ("rat_mult_poly", rat_mult_poly),

  3700    ("rat_mult_div_pow", rat_mult_div_pow),

  3701    ("cancel_p_rls", cancel_p_rls)

  3702    ]);

  3704 calclist':= overwritel (!calclist',

  3705    [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))

  3706     ]);

  3708 (** problems **)

  3710 store_pbt

  3711  (prep_pbt thy "pbl_simp_rat" [] e_pblID

  3712  (["rational","simplification"],

  3713   [("#Given" ,["Term t_t"]),

  3714    ("#Where" ,["t_t is_ratpolyexp"]),

  3715    ("#Find"  ,["normalform n_n"])

  3716   ],

  3717   append_rls "e_rls" e_rls [(*for preds in where_*)],

  3718   SOME "Simplify t_t",

  3719   [["simplification","of_rationals"]]));

  3721 (** methods **)

  3723 (*WN061025 this methods script is copied from (auto-generated) script

  3724   of norm_Rational in order to ease repair on inform*)

  3725 store_met (prep_met thy "met_simp_rat" [] e_metID (["simplification","of_rationals"],

  3726 	  [("#Given" ,["Term t_t"]),

  3727 	   ("#Where" ,["t_t is_ratpolyexp"]),

  3728 	   ("#Find"  ,["normalform n_n"])],

  3729 	  {rew_ord'="tless_true", rls' = e_rls, calc = [], srls = e_rls,

  3730 	   prls = append_rls "simplification_of_rationals_prls" e_rls

  3731 				    [(*for preds in where_*) Calc ("Rational.is'_ratpolyexp", eval_is_ratpolyexp "")],

  3732 				   crls = e_rls, errpats = [],

  3733    nrls = norm_Rational_rls},

  3734    "Script SimplifyScript (t_t::real) =                             " ^

  3735    "  ((Try (Rewrite_Set discard_minus False) @@                   " ^

  3736    "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^

  3737    "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^

  3738    "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^

  3739    "    (Repeat                                                     " ^

  3740    "     ((Try (Rewrite_Set common_nominator_p_rls False) @@        " ^

  3741    "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^

  3742    "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^

  3743    "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^

  3744    "       Try (Rewrite_Set rat_reduce_1 False)))) @@               " ^

  3745    "    Try (Rewrite_Set discard_parentheses1 False))               " ^

  3746    "   t_t)"));

  3748 *}

  3749 ML {*"test"*}

  3750 end (* theory*)

author	Walther Neuper <neuper@ist.tugraz.at>
	Mon, 16 Sep 2013 10:46:51 +0200
changeset 52103	0d13f07d8e2a
parent 52101	c3f399ce32af
child 52105	2786cc9704c8
permissions	-rwxr-xr-x