(*.calculate in rationals: gcd, lcm, etc.

    (c) Stefan Karnel 2002

    Institute for Mathematics D and Institute for Software Technology, 

    TU-Graz SS 2002 

    Use is subject to license terms.

 use"IsacKnowledge/Rational.ML";

 use"Rational.ML";

 remove_thy"Rational";

 use_thy"IsacKnowledge/Isac";

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

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

   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].

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

 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

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

 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):BasisLibrary.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) =  BasisLibrary.Real.fromInt(x)*BasisLibrary.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:=BasisLibrary.Real.ceil(2.0 *   

 				  BasisLibrary.Real.fromInt(!d) *

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

 				  BasisLibrary.Real.fromInt(!d) * 

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

 					uv_mod_norm(!p2) / BasisLibrary.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)=print("[]\n")

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

   | mv_print_poly((x,y)::p1) = (print("("^BasisLibrary.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

   | uv_to_poly p1 = 

     let

 	val count=ref 0;

 	val help=ref p1;

 	val list=ref [];

     in

 	while length(!help)>0 do

 	    (

 	     if hd(!help)=0 then ()

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

 	     count:=(!count)+1;

 	     help:=tl(!help)

 	     );

 	    (!list)

     end;

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

 fun uv_gcd ([]:mv_poly) p2 = p2

   | uv_gcd p1 ([]:mv_poly) = p1

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

     let 

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

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

     in

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

     end

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

     let 

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

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

     in

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

     end 

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

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

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

   | mv_newton_help (pl:mv_poly list,k) = 

     let

 	val x=ref (rev(pl));

 	val t=ref [];

 	val y=ref [];

 	val n=ref 1;

 	val n1=ref[];

     in

 	(

 	 while length(!x)>1 do 

 	     (

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

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

 		   else n1:=[]; 

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

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

 	      x:=tl(!x)

 	      );

 	 (!y)

 	 )

     end;

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

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

   | mv_newton_add [x:mv_poly] t = x 

   | mv_newton_add (pl:mv_poly list) t = 

     let

 	val expos=ref [];

 	val pll=ref pl;

     in

 	(

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

 	     ( 

 	      pll:=tl(!pll)

 	      ); 

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

 	 mv_add(hd(pl),

 		mv_mul(

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

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

 		       LEX_

 		       ),

 		LEX_)

 	 )

     end;

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

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

 (*. this function returns [] .*)

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

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

   | mv_newton pl =

     let

 	val c=ref pl;

 	val c1=ref [];

 	val n=length(pl);

 	val k=ref 1;

 	val i=ref n;

 	val ppl=ref [];

     in

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

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

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

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

 	    (	 

 	     k:=(!k)+1; 

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

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

 	     ppl:= !c;

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

 	     );

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

 	c1:=rev(!c1);

 	ppl:= !c1;

 	mv_newton_add (!c1) 1

     end;

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

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

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

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

 fun mv_min_pp([]:mv_poly)=[]

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

     let

 	val y=ref xs;

 	val x=ref [];

     in

 	(

 	 x:=e;

 	 while length(!y)>0 do

 	     (

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

 	      y:=tl(!y)

 	      );

 	 !x

 	 )

     end;

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

 fun list_is_null [] = true 

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

 (* check if main variable is zero*)

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

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

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

   | mv_content(p1) = 

     let

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

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

 	val result=ref []; 

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

     in

 	(

 	 if length(!list)>1 then

 	     (

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

 		  (

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

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

 		   else list:=tl(!list) 

 		       );		  

 		  if length(!list)>0 then  

 		   ( 

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

 		    ) 

 		  else list:=(!result); 

 		  list:=mv_correct(!list,0);  

 		  (!list) 

 		  )

 	 else

 	     (

 	      mv_null3(!list) 

 	      )

 	     )

     end

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

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

   | mv_pp(p1) = let

 		    val cont=ref []; 

 		    val pp=ref[];

 		in

 		    cont:=mv_content(p1);

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

 		    if !pp=[] 

 			then raise error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")

 		    else (!pp)

 		end

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

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

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

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

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

      let

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

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

      in 

 	(

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

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

         )

     end 

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

 	(