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

   author: isac team

   Copyright (c) isac team 2002

   Use is subject to license terms.

   depends on Poly (and not on Atools), because 

   fractions with _normalized_ polynomials are canceled, added, etc.

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

*)

theory Rational 

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

begin

consts

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

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

  get_denominator :: "real => real"

  get_numerator   :: "real => real"

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

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

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

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

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

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

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

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

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

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

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

(*real_times_divide1_eq .. Isa02*) 

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

  real_times_divide_num:   "a is_const ==> 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

section {* Cancellation and addition of fractions *}

subsection {* Auxiliary functions and data *}

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

ML {*

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

    if is_numeral id 

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

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

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

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

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

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

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

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

fun monoms_of_term vs (t as Free _) =

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

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

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

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

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

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

    (monoms_of_term vs ms1) @ (monoms_of_term vs ms2)

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

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

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

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

* the coefficient may be somewhere

* he order and the parentheses within monomials are arbitrary

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

The function requires the free variables as strings already given, 

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

*)

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

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

      handle ERROR _ => NONE)

  | poly_of_term vs t =

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

      handle ERROR _ => NONE)

fun is_poly t =

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

  in 

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

  end

val is_expanded = is_poly

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

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

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

    [] @ term_of_es baseT expT vs es

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

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

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

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

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

    @ term_of_es baseT expT vs es

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

    let val es' = term_of_es baseT expT vs es

    in 

      if c = 1 

      then 

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

        then Free (isastr_of_int c, baseT)

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

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

    end

fun term_of_poly baseT expT vs p =

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

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

*}

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

      (c) Stefan Karnel 2002

      Institute for Mathematics D and Institute for Software Technology, 

      TU-Graz SS 2002 *}

text {* Remark on notions in the documentation below:

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

    [2] 'polynomial' normalform (Polynom)

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

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

    Instead of 

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

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

    we say: 

      'fraction in normalform [2]'

      'fraction in normalform [3]' 

    or

      'fraction [2]'

      'fraction [3]'.

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

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

*}

ML {* 

val thy = @{theory};

signature RATIONALI =

sig

  type mv_monom

  type mv_poly 

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

  val calculate_Rational : rls

  val calc_rat_erls:rls

  val cancel_p : rls   

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

  val common_nominator_p : rls              

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

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

			 (string * term) option                    

  val expanded2polynomial : term -> term option

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

  val is_expanded : term -> bool

  val is_polynomial : term -> bool

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

  val mv_lcm : mv_poly -> mv_poly -> mv_poly

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

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

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

                           roots into a simple and canceled fraction

                           with normalform [2].*)

*)

(*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME

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

                           roots into a simple and canceled fraction

                           with normalform [2].*)

*)

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

  val polynomial2expanded : term -> term option

  val rational_erls : 

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

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

end (* sig *)

(*.**************************************************************************

survey on the functions

~~~~~~~~~~~~~~~~~~~~~~~

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

--------------------:------------------+-------------------:-----------------

 factout_p_         :                  | factout_          :

 cancel_p_          :                  | cancel_           :

                    :cancel_p          |                   :cancel

--------------------:------------------+-------------------:-----------------

 common_nominator_p_:                  | common_nominator_ :

                    :common_nominator_p|                   :common_nominator

 add_fraction_p_    :                  | add_fraction_     :

--------------------:------------------+-------------------:-----------------

???SK                 :norm_rational_p   |                   :norm_rational

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

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

conversion from Isabelle-term to internal representation

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

... BITTE FORTSETZEN ...

polynomial2expanded = ...

expanded2polynomial = ...

remark: polynomial2expanded o expanded2polynomial = I, 

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

functions for greatest common divisor and canceling

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

################################################################################

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

##   Amine Chaieb, Robert Himmelmann, John Harrison.

################################################################################

mv_gcd

factout_

factout_p_

cancel_

cancel_p_

functions for least common multiple and addition of fractions

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

mv_lcm

common_nominator_

common_nominator_p_

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

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

functions for normalform of rationals

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

          schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?

norm_rational_

norm_expanded_rat_

**************************************************************************.*)

(*##*)

structure RationalI : RATIONALI = 

struct 

(*##*)

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

fun x mem [] = false

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

val is_expanded = is_poly

val is_polynomial = is_poly

fun mk_noteq_0 baseT t = 

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

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

fun mk_asms baseT ts =

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

  in map (mk_noteq_0 baseT) as' end

fun check_fraction t =

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

  in SOME (numerator, denominator) end

  handle Bind => NONE

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

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

fun factout_p_ (thy: theory) t =

  let val opt = check_fraction t

  in

    case opt of 

      NONE => NONE

    | SOME (numerator, denominator) =>

      let 

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

        val baseT = type_of numerator

        val expT = HOLogic.realT

      in

        case (poly_of_term vs numerator, poly_of_term vs denominator) of

          (SOME a, SOME b) =>

            let

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

              val es = replicate (length vs) 0 

            in

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

              then NONE

              else 

                let

                  val b't = term_of_poly baseT expT vs b'

                  val ct = term_of_poly baseT expT vs c

                  val t' = 

                    HOLogic.mk_binop "Fields.inverse_class.divide" 

                      (HOLogic.mk_binop "Groups.times_class.times"

                        (term_of_poly baseT expT vs a', ct),

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

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

                in SOME (t', asm) end

            end

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

      end

  end

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

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

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

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

  yields

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

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

    \<or> NONE *)

fun cancel_p_ (thy: theory) t =

  let val opt = check_fraction t

  in

    case opt of 

      NONE => NONE

    | SOME (numerator, denominator) =>

      let 

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

        val baseT = type_of numerator

        val expT = HOLogic.realT

      in

        case (poly_of_term vs numerator, poly_of_term vs denominator) of

          (SOME a, SOME b) =>

            let

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

              val es = replicate (length vs) 0 

            in

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

              then NONE

              else 

                let

                  val bt' = term_of_poly baseT expT vs b'

                  val ct = term_of_poly baseT expT vs c

                  val t' = 

                    HOLogic.mk_binop "Fields.inverse_class.divide" 

                      (term_of_poly baseT expT vs a', bt')

                  val asm = mk_asms baseT [bt']

                in SOME (t', asm) end

            end

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

      end

  end

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

fun norm_frac_sum 

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

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

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

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

  | norm_frac_sum 

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

      nofrac $ 

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

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

  | norm_frac_sum 

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

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

      nofrac)

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

  | norm_frac_sum _ = NONE  

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

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

fun common_nominator_p_ (thy: theory) t =

  let val opt = norm_frac_sum t

  in

    case opt of 

      NONE => NONE

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

      let 

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

        val baseT = type_of n1

        val expT = HOLogic.realT

      in

        case (poly_of_term vs d1, poly_of_term vs d2) of

          (SOME a, SOME b) =>

            let

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

              val d1' = term_of_poly baseT expT vs a'

              val d2' = term_of_poly baseT expT vs b'

              val c' = term_of_poly baseT expT vs c

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

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

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

                --------------------------------------------------------------------------*)

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

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

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

              val t' =

                HOLogic.mk_binop "Groups.plus_class.plus"

                  (HOLogic.mk_binop "Fields.inverse_class.divide"

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

                  HOLogic.mk_binop "Fields.inverse_class.divide" 

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

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

            in SOME (t', asm) end

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

      end

  end

(* add fractions

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

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

fun add_fraction_p_ (thy: theory) t =

  let val opt = norm_frac_sum t

  in

    case opt of 

      NONE => NONE

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

      let 

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

        val baseT = type_of n1

        val expT = HOLogic.realT

      in

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

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

            let

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

              val nomin = term_of_poly baseT expT vs 

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

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

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

              val asm = mk_asms baseT [denom]

            in SOME (t', asm) end

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

      end

  end

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

fun xs union [] = xs

  | [] union ys = ys

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

(*. gcd of integers .*)

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

fun gcd_int a b = if b=0 then a

		  else gcd_int b (a mod b);

(*. univariate polynomials (uv) .*)

(*. univariate polynomials are represented as a list 

    of the coefficent in reverse maximum degree order .*)

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

type uv_poly = int list;

(*. adds two uv polynomials .*)

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

  | uv_mod_add_poly (p1,[]) = p1

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

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

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

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

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

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

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

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

fun uv_mod_deg ([]:uv_poly) = 0  

  | uv_mod_deg p = length(p)-1;

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

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

fun uv_mod_mod2(x,p)=

    let

	val y=(x mod p);

    in

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

	    (

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

	     )

    end;

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

fun uv_mod_list_modp [] p = [] 

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

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

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

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

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

(*. only for uv_mod_divides .*)

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

   integer zu klein  *)

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

    error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

  | uv_mod_pdiv p1 [x] = 

    let

	val xs= Unsynchronized.ref  [];

    in

	if x<>0 then 

	    (

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

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

	     )

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

	([]:uv_poly,!xs:uv_poly)

    end

  | uv_mod_pdiv p1 p2 =  

    let

	val n= uv_mod_deg(p2);

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

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

	val p2'=(rev(p2));

	val lc2=hd(p2');

	val q= Unsynchronized.ref  [];

	val c= Unsynchronized.ref  0;

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

    in

	(

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

         then error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero") 

	 else

	     (

	      if (!m)<n then 

		  (

		   output:=([0],p1) 

		   ) 

	      else

		  (

		   while (!m)>=n do

		       (

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

			if !c<>0 then

			    (

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

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

			     m:=uv_mod_deg(!p1')

			     )

			else m:=0

			);

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

		   )

	      );

	     !output

	 )

    end;

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

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

    let

	val n=uv_mod_deg(p2);

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

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

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

	val lc2=hd(p2');

	val q= Unsynchronized.ref  [];

	val c= Unsynchronized.ref  0;

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

    in

	(

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

	 else

	     (

	      if (!m)<n then 

		  (

		   output:=([0],p1) 

		   ) 

	      else

		  (

		   while (!m)>=n do

		       (

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

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

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

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

			m:=(!m)-1

			);

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

		       (

			p1':=tl(!p1')

			);

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

		   )

	      );

	     !output:uv_poly * uv_poly

	 )

    end;

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

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

  | uv_mod_prest [] p2 = []:uv_poly

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

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

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

  | uv_mod_prestp [] p2 p= []:uv_poly 

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

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

fun uv_mod_cont ([]:uv_poly) = 0  

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

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

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

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

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

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

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

  | uv_mod_pp p =  

    let

	val c= Unsynchronized.ref  0;

    in

	(

	 c:=uv_mod_cont(p);

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

	 else uv_mod_div_list(p,!c)

	)

    end;

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

fun uv_mod_lc ([]:uv_poly) = 0 

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

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

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

    let

	val f = Unsynchronized.ref  [];

	val f'= Unsynchronized.ref  p2;

	val fi= Unsynchronized.ref  [];

    in

	( 

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

 	 while uv_mod_deg(!f')>0 do

	     (

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

	      if (!f')<>[] then 

		  (

		   fi:=(!f');

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

		   )

	      else ()

	      );

	      (!f)

	 )

    end;

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

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

  | uv_mod_gcd_modp p1 [] p= p1

  | uv_mod_gcd_modp p1 p2 p=

    let

	val p1'= Unsynchronized.ref [];

	val p2'= Unsynchronized.ref [];

	val pc= Unsynchronized.ref [];

	val g= Unsynchronized.ref  [];

	val d= Unsynchronized.ref  0;

	val prs= Unsynchronized.ref  [];

    in

	(

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

	     (

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

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

	      )

	 else 

	     (

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

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

	      );

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

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

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

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

	 else pc:=hd(!prs);

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

	 !g

	 )

    end;

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

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

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

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

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

fun uv_mod_add_qu [] = 0.0 

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

(*. calculates the euklidean norm .*)

fun uv_mod_norm ([]:uv_poly) = 0.0

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

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

fun uv_mod_multi a b = a * b;

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

fun uv_mod_prim(x,[])= false 

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

		else false

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

			then 

			    if uv_mod_prim(x,ys) then true 

			    else false

		    else false;

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

fun uv_mod_nextprime(g,p)= 

    let

	val list= Unsynchronized.ref  [2];

	val exit= Unsynchronized.ref  0;

	val i = Unsynchronized.ref 2

    in

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

	    (

	     if uv_mod_prim(!i,!list) then

		 (

		  if (p mod !i <> 0)

		      then

			  (

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

			   i:= (!i)+1

			   )

		  else i:=(!i)+1

		  )

	     else i:= (!i)+1

		 );

	    i:=(p+1);

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

	    (

	     if uv_mod_prim(!i,!list) then

		 (

		  if (g mod !i <> 0)

		      then

			  (

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

			   exit:= (!i)

			   )

		  else i:=(!i)+1

		      )

	     else i:= (!i)+1

		 ); 

	    !exit

    end;

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

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

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

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

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

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

(*. chinese remainder algorithm .*)

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

    let 

	val c= Unsynchronized.ref  0;

	val r1'= Unsynchronized.ref  0;

	val d= Unsynchronized.ref  0;

	val a= Unsynchronized.ref  0;

    in

	(

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

	     (

	      c:=(!c)+1

	      );

	 r1':= uv_mod_mod2(r1,m1);

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

	 !r1'+(!d)*m1    

	 )

    end;

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

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

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

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

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

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

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

    let 

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

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

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

	val temp= Unsynchronized.ref  [];

	val cp= Unsynchronized.ref  [];

	val qp= Unsynchronized.ref  [];

	val q= Unsynchronized.ref [];

	val pn= Unsynchronized.ref  0;

	val d= Unsynchronized.ref  0;

	val g1= Unsynchronized.ref  0;

	val p= Unsynchronized.ref  0;    

	val m= Unsynchronized.ref  0;

	val exit= Unsynchronized.ref  0;

	val i= Unsynchronized.ref  1;

    in

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

	else 

	    (

	     temp:= !p1;

	     p1:= !p2;

	     p2:= !temp

	     );

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

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

	p:=4;

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

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

	  Real.fromInt(!d) * 

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

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

	while (!exit=0) do  

	    (

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

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

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

		 (

		  i:=1;

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

		      (

		       i:=(!i)+1

		       );

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

		  )

	     else ();

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

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

	     pn:=(!p);

	     q:=(!qp);

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

		 (

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

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

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

 		      (

 		       i:=1;

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

 			   (

 			    i:=(!i)+1

		           );

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

 		      )

 		  else ();    

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

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

 		      (

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

 		       )

 		  else

 		      (

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

 			   (

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

			    pn:=(!pn) * !p;

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

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

		 	    )

		       else

			   (

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

				(

				 pn:= !p;

				 q:= !qp

				 )

			    else ()

			    )

		       )

		  );

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

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

	     );

	    uv_mod_smul_poly(!q,c):uv_poly

    end;

(*. multivariate polynomials .*)

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

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

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

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

(*. global variables .*)

(*. order indicators .*)

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

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

(*. datatypes for internal representation.*)

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

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

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

type mv_poly = mv_monom list; 

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

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

fun mv_mg_hlp([]) = EQUAL 

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

		    else if x>0 then GREATER

			 else mv_mg_hlp(list);

(*. adds a list of values .*)

fun mv_addlist([]) = 0

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

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

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

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

    if order=LEX_ then

	( 

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

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

	     )

    else

	if order=GGO_ then

	    ( 

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

	     else 

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

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

	     )

	else error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");

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

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

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

let 

    val temp= Unsynchronized.ref  EQUAL;

in

    if order=LEX_ then

	(

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

	     error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

	 else 

	     (

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

	      if !temp=EQUAL then 

		  ( if x1=x2 then EQUAL 

		    else if x1>x2 then GREATER

			 else LESS

			     )

	      else (!temp)

	      )

	     )

    else 

	if order=GGO_ then 

	    (

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

	      error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

	     else