(* 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.

*)

theory Rational imports Poly begin

consts

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

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

axioms (*.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)"

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

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

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

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

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

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

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

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

(*real_times_divide1_eq .. Isa02*) 

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

  real_times_divide_num:   "a is_const ==> 

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

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

(*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)"

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

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

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

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

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

  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"

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

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

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

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

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

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

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

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

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

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

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

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

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

val ------------------------------------------------------+ = "11111";

signature RATIONALI =

sig

  type mv_monom

  type mv_poly 

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

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

  val calculate_Rational : rls

  val calc_rat_erls:rls

  val cancel : rls

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

  val cancel_p : rls   

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

  val common_nominator : rls              

  val common_nominator_ : 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_ : theory -> term -> (term * term list) 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

val ------------------------------------------------------++ = "22222";

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

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

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

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;

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

    raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

  | uv_mod_pdiv p1 [x] = 

    let

	val xs=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 raise 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= ref (uv_mod_deg(p1));

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

	val p2'=(rev(p2));

	val lc2=hd(p2');

	val q=ref [];

	val c=ref 0;

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

    in

	(

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

         then raise 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=ref (uv_mod_deg(uv_mod_list_modp p1 p));

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

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

	val lc2=hd(p2');

	val q=ref [];

	val c=ref 0;

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

    in

	(

	 if (!m)=0 orelse p2=[0] then raise 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) = raise 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= raise 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) = raise 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=ref 0;

    in

	(

	 c:=uv_mod_cont(p);

	 if !c=0 then raise 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 =ref [];

	val f'=ref p2;

	val fi=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'=ref[];

	val p2'=ref[];

	val pc=ref[];

	val g=ref [];

	val d=ref 0;

	val prs=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=ref [2];

	val exit=ref 0;

	val i = 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= raise 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) = raise 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=ref 0;

	val r1'=ref 0;

	val d=ref 0;

	val a=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) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x1")

  | uv_mod_cra_2 (x1,[],m1,m2) = raise 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=ref (uv_mod_pp(p1'));

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

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

	val temp=ref [];

	val cp=ref [];

	val qp=ref [];

	val q=ref[];

	val pn=ref 0;

	val d=ref 0;

	val g1=ref 0;

	val p=ref 0;    

	val m=ref 0;

	val exit=ref 0;

	val i=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 raise 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 raise 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 raise 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=ref EQUAL;

in

    if order=LEX_ then

	(

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

	     raise 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 

	      raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

	     else 

		 if mv_addlist(x)=mv_addlist(y) then 

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

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

		     )

	else raise error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")

end;

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

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

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

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

(*. leading power product .*)

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

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

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

(*. leading monomial .*)

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

  | mv_lm([x],order) = x 

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

(*. leading coefficient in term order .*)

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

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

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

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

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

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

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

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

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

  | mv_lc(p1,order)  = 

    let

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

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

	val lc=ref [];

    in

	(

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

	     (

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

	      p1o:=tl(!p1o)

	      );

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

	 mv_rev_to(!lc)

	 )

    end;

(*. compares two powerproducts .*)

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

(*. help function for mv_add .*)

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

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

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

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

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

  | mv_madd(p1,p2,order) = 

    (

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

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

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

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

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

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

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

	      )

(*. adds two multivariate polynomials .*)

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

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

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

(*. monom multiplication .*)

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

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

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

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

(*. zeros a list .*)

fun mv_null2([])=[]

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

(*. multiplies two multivariate polynomials .*)

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

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

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

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

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

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

fun mv_getmax([])=0

  | mv_getmax(x::p1)= let 

		       val m=mv_getmax(p1);

		   in

		       if m>x then m

		       else x

		   end;

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

fun mv_grad([]:mv_poly) = 0 

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

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

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

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

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

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

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

fun mv_is_negativ([])=false

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

(*. division of monomials .*)

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

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

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

    let

	val c=ref (#1(p2));

	val pp=ref [];

    in 

	(

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

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

	     if #1(p2)<>0 then 

		 (

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

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

		  else (!c,!pp) 

		      )

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

		 )

    end;

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

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

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

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

(*. division of two multivariate polynomials .*) 

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

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

  | mv_division(f,g,order)=

    let 

	val r=ref [];

	val q=ref [];

	val g'=ref [];

	val k=ref 0;

	val m=ref (0,[0]);

	val exit=ref 0;

    in

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

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

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

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

	else

	    (

	     exit:=0;

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

		 (

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

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

		  if #1(!m)<>0 then

		      ( 

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

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

		       )

		  else exit:=1;

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

		  else (exit:=1)

		  );

		 (rev(!q),!r)

		 )

    end;

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

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

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

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

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

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

(*. multivariate case .*)

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

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

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

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

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

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

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

fun mv_deg([]:mv_poly) = 0  

  | mv_deg(p1)=

    let

	val p1'=mv_shorten(p1,LEX_);

    in

	if length(p1')=0 then 0 

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

    end;

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

fun mv_deg2([]:mv_poly) = 0

  | mv_deg2(p1)=

    let

	val p1'=mv_shorten(p1,LEX_);

    in

	if length(p1')=0 then 0 

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

    end;

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

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

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

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

fun uv_content2([]:mv_poly) = 0

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

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

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

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

    let

	val count=ref 0;

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

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

	val list=ref [];

    in

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

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

	     else

		 (

		  count:=0;

		  while (!count)<=max do

		      (

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

			   (

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

			    help:=tl(!help) 

			    )

		       else 

			   (

			    list:= 0::(!list)

			    );

		       count := (!count) + 1

		       );

		      (!list)

		      )

    end;

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

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