wneuper/isa: src/Tools/isac/Knowledge/Rational.thy@525a28152a67

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

     2    author: isac team

     3    Copyright (c) isac team 2002

     4    Use is subject to license terms.

     6    depends on Poly (and not on Atools), because

     7    fractions with _normalized_ polynomials are canceled, added, etc.

     8 *)

    10 theory Rational imports Poly begin

    12 consts

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

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

    17 axioms (*.not contained in Isabelle2002,

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

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

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

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

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

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

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

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

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

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

    33 (*real_times_divide1_eq .. Isa02*)

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

    35   real_times_divide_num   "a is_const ==>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    66       (c) Stefan Karnel 2002

    67       Institute for Mathematics D and Institute for Software Technology,

    68       TU-Graz SS 2002 *}

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

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

    72     [2] 'polynomial' normalform (Polynom)

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

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

    75     Instead of

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

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

    78     we say:

    79       'fraction in normalform [2]'

    80       'fraction in normalform [3]'

    81     or

    82       'fraction [2]'

    83       'fraction [3]'.

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

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

    86 *}

    88 ML {*

    89 signature RATIONALI =

    90 sig

    91   type mv_monom

    92   type mv_poly

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

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

    95   val calculate_Rational : rls

    96   val calc_rat_erls:rls

    97   val cancel : rls

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

    99   val cancel_p : rls

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

   101   val common_nominator : rls

   102   val common_nominator_ : theory -> term -> (term * term list) option

   103   val common_nominator_p : rls

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

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

   106 			 (string * term) option

   107   val expanded2polynomial : term -> term option

   108   val factout_ : theory -> term -> (term * term list) option

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

   110   val is_expanded : term -> bool

   111   val is_polynomial : term -> bool

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

   114   val mv_lcm : mv_poly -> mv_poly -> mv_poly

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

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

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

   119                            roots into a simple and canceled fraction

   120                            with normalform [2].*)

   121 *)

   122 (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME

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

   124                            roots into a simple and canceled fraction

   125                            with normalform [2].*)

   126 *)

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

   128   val polynomial2expanded : term -> term option

   129   val rational_erls :

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

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

   133 end

   135 (*.**************************************************************************

   136 survey on the functions

   137 ~~~~~~~~~~~~~~~~~~~~~~~

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

   139 --------------------:------------------+-------------------:-----------------

   140  factout_p_         :                  | factout_          :

   141  cancel_p_          :                  | cancel_           :

   142                     :cancel_p          |                   :cancel

   143 --------------------:------------------+-------------------:-----------------

   144  common_nominator_p_:                  | common_nominator_ :

   145                     :common_nominator_p|                   :common_nominator

   146  add_fraction_p_    :                  | add_fraction_     :

   147 --------------------:------------------+-------------------:-----------------

   148 ???SK                 :norm_rational_p   |                   :norm_rational

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

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

   154 conversion from Isabelle-term to internal representation

   155 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   157 ... BITTE FORTSETZEN ...

   159 polynomial2expanded = ...

   160 expanded2polynomial = ...

   162 remark: polynomial2expanded o expanded2polynomial = I,

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

   165 functions for greatest common divisor and canceling

   166 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   167 mv_gcd

   168 factout_

   169 factout_p_

   170 cancel_

   171 cancel_p_

   173 functions for least common multiple and addition of fractions

   174 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   175 mv_lcm

   176 common_nominator_

   177 common_nominator_p_

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

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

   181 functions for normalform of rationals

   182 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

   184           schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ?

   186 norm_rational_

   187 norm_expanded_rat_

   189 **************************************************************************.*)

   192 (*##*)

   193 structure RationalI : RATIONALI =

   194 struct

   195 (*##*)

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

   198 fun x mem [] = false

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

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

   201 fun xs union [] = xs

   202   | [] union ys = ys

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

   205 (*. gcd of integers .*)

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

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

   208 		  else gcd_int b (a mod b);

   210 (*. univariate polynomials (uv) .*)

   211 (*. univariate polynomials are represented as a list

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

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

   214 type uv_poly = int list;

   216 (*. adds two uv polynomials .*)

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

   218   | uv_mod_add_poly (p1,[]) = p1

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

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

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

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

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

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

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

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

   230 fun uv_mod_deg ([]:uv_poly) = 0

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

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

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

   235 fun uv_mod_mod2(x,p)=

   236     let

   237 	val y=(x mod p);

   238     in

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

   240 	    (

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

   242 	     )

   243     end;

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

   246 fun uv_mod_list_modp [] p = []

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

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

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

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

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

   254 (*. only for uv_mod_divides .*)

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

   256    integer zu klein  *)

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

   258     raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero")

   259   | uv_mod_pdiv p1 [x] =

   260     let

   261 	val xs=ref [];

   262     in

   263 	if x<>0 then

   264 	    (

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

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

   267 	     )

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

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

   270     end

   271   | uv_mod_pdiv p1 p2 =

   272     let

   273 	val n= uv_mod_deg(p2);

   274 	val m= ref (uv_mod_deg(p1));

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

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

   277 	val lc2=hd(p2');

   278 	val q=ref [];

   279 	val c=ref 0;

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

   281     in

   282 	(

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

   284          then raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero")

   285 	 else

   286 	     (

   287 	      if (!m)<n then

   288 		  (

   289 		   output:=([0],p1)

   290 		   )

   291 	      else

   292 		  (

   293 		   while (!m)>=n do

   294 		       (

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

   296 			if !c<>0 then

   297 			    (

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

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

   300 			     m:=uv_mod_deg(!p1')

   301 			     )

   302 			else m:=0

   303 			);

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

   305 		   )

   306 	      );

   307 	     !output

   308 	 )

   309     end;

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

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

   313     let

   314 	val n=uv_mod_deg(p2);

   315 	val m=ref (uv_mod_deg(uv_mod_list_modp p1 p));

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

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

   318 	val lc2=hd(p2');

   319 	val q=ref [];

   320 	val c=ref 0;

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

   322     in

   323 	(

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

   325 	 else

   326 	     (

   327 	      if (!m)<n then

   328 		  (

   329 		   output:=([0],p1)

   330 		   )

   331 	      else

   332 		  (

   333 		   while (!m)>=n do

   334 		       (

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

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

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

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

   339 			m:=(!m)-1

   340 			);

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

   343 		       (

   344 			p1':=tl(!p1')

   345 			);

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

   348 		   )

   349 	      );

   350 	     !output:uv_poly * uv_poly

   351 	 )

   352     end;

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

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

   356   | uv_mod_prest [] p2 = []:uv_poly

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

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

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

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

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

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

   365 fun uv_mod_cont ([]:uv_poly) = 0

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

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

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

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

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

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

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

   375   | uv_mod_pp p =

   376     let

   377 	val c=ref 0;

   378     in

   379 	(

   380 	 c:=uv_mod_cont(p);

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

   383 	 else uv_mod_div_list(p,!c)

   384 	)

   385     end;

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

   388 fun uv_mod_lc ([]:uv_poly) = 0

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

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

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

   393     let

   394 	val f =ref [];

   395 	val f'=ref p2;

   396 	val fi=ref [];

   397     in

   398 	(

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

   400  	 while uv_mod_deg(!f')>0 do

   401 	     (

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

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

   404 		  (

   405 		   fi:=(!f');

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

   407 		   )

   408 	      else ()

   409 	      );

   410 	      (!f)

   412 	 )

   413     end;

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

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

   417   | uv_mod_gcd_modp p1 [] p= p1

   418   | uv_mod_gcd_modp p1 p2 p=

   419     let

   420 	val p1'=ref[];

   421 	val p2'=ref[];

   422 	val pc=ref[];

   423 	val g=ref [];

   424 	val d=ref 0;

   425 	val prs=ref [];

   426     in

   427 	(

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

   429 	     (

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

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

   432 	      )

   433 	 else

   434 	     (

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

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

   437 	      );

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

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

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

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

   444 	 else pc:=hd(!prs);

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

   447 	 !g

   448 	 )

   449     end;

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

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

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

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

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

   458 fun uv_mod_add_qu [] = 0.0

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

   461 (*. calculates the euklidean norm .*)

   462 fun uv_mod_norm ([]:uv_poly) = 0.0

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

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

   466 fun uv_mod_multi a b = a * b;

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

   469 fun uv_mod_prim(x,[])= false

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

   471 		else false

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

   473 			then

   474 			    if uv_mod_prim(x,ys) then true

   475 			    else false

   476 		    else false;

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

   479 fun uv_mod_nextprime(g,p)=

   480     let

   481 	val list=ref [2];

   482 	val exit=ref 0;

   483 	val i = ref 2

   484     in

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

   486 	    (

   487 	     if uv_mod_prim(!i,!list) then

   488 		 (

   489 		  if (p mod !i <> 0)

   490 		      then

   491 			  (

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

   493 			   i:= (!i)+1

   494 			   )

   495 		  else i:=(!i)+1

   496 		  )

   497 	     else i:= (!i)+1

   498 		 );

   499 	    i:=(p+1);

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

   501 	    (

   502 	     if uv_mod_prim(!i,!list) then

   503 		 (

   504 		  if (g mod !i <> 0)

   505 		      then

   506 			  (

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

   508 			   exit:= (!i)

   509 			   )

   510 		  else i:=(!i)+1

   511 		      )

   512 	     else i:= (!i)+1

   513 		 );

   514 	    !exit

   515     end;

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

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

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

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

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

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

   525 (*. chinese remainder algorithm .*)

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

   527     let

   528 	val c=ref 0;

   529 	val r1'=ref 0;

   530 	val d=ref 0;

   531 	val a=ref 0;

   532     in

   533 	(

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

   535 	     (

   536 	      c:=(!c)+1

   537 	      );

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

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

   540 	 !r1'+(!d)*m1

   541 	 )

   542     end;

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

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

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

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

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

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

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

   552     let

   553 	val p1=ref (uv_mod_pp(p1'));

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

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

   556 	val temp=ref [];

   557 	val cp=ref [];

   558 	val qp=ref [];

   559 	val q=ref[];

   560 	val pn=ref 0;

   561 	val d=ref 0;

   562 	val g1=ref 0;

   563 	val p=ref 0;

   564 	val m=ref 0;

   565 	val exit=ref 0;

   566 	val i=ref 1;

   567     in

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

   569 	else

   570 	    (

   571 	     temp:= !p1;

   572 	     p1:= !p2;

   573 	     p2:= !temp

   574 	     );

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

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

   579 	p:=4;

   581 	m:=BasisLibrary.Real.ceil(2.0 *

   582 				  BasisLibrary.Real.fromInt(!d) *

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

   584 				  BasisLibrary.Real.fromInt(!d) *

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

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

   588 	while (!exit=0) do

   589 	    (

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

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

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

   593 		 (

   594 		  i:=1;

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

   596 		      (

   597 		       i:=(!i)+1

   598 		       );

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

   600 		  )

   601 	     else ();

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

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

   607 	     pn:=(!p);

   608 	     q:=(!qp);

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

   611 		 (

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

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

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

   615  		      (

   616  		       i:=1;

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

   618  			   (

   619  			    i:=(!i)+1

   620 		           );

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

   622  		      )

   623  		  else ();

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

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

   627  		      (

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

   629  		       )

   630  		  else

   631  		      (

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

   633  			   (

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

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

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

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

   638 		 	    )

   639 		       else

   640 			   (

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

   642 				(

   643 				 pn:= !p;

   644 				 q:= !qp

   645 				 )

   646 			    else ()

   647 			    )

   648 		       )

   649 		  );

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

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

   652 	     );

   653 	    uv_mod_smul_poly(!q,c):uv_poly

   654     end;

   656 (*. multivariate polynomials .*)

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

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

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

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

   662 (*. global variables .*)

   663 (*. order indicators .*)

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

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

   667 (*. datatypes for internal representation.*)

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

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

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

   672 type mv_poly = mv_monom list;

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

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

   676 fun mv_mg_hlp([]) = EQUAL

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

   678 		    else if x>0 then GREATER

   679 			 else mv_mg_hlp(list);

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

   682 fun mv_addlist([]) = 0

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

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

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

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

   688     if order=LEX_ then

   689 	(

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

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

   692 	     )

   693     else

   694 	if order=GGO_ then

   695 	    (

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

   697 	     else

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

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

   700 	     )

   701 	else raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order");

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

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

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

   706 let

   707     val temp=ref EQUAL;

   708 in

   709     if order=LEX_ then

   710 	(

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

   712 	     raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   713 	 else

   714 	     (

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

   716 	      if !temp=EQUAL then

   717 		  ( if x1=x2 then EQUAL

   718 		    else if x1>x2 then GREATER

   719 			 else LESS

   720 			     )

   721 	      else (!temp)

   722 	      )

   723 	     )

   724     else

   725 	if order=GGO_ then

   726 	    (

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

   728 	      raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error")

   729 	     else

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

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

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

   733 		     )

   734 	else raise error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order")

   735 end;

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

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

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

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

   742 (*. leading power product .*)

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

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

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

   747 (*. leading monomial .*)

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

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

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

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

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

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

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

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

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

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

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

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

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

   765   | mv_lc(p1,order)  =

   766     let

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

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

   769 	val lc=ref [];

   770     in

   771 	(

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

   773 	     (

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

   775 	      p1o:=tl(!p1o)

   776 	      );

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

   778 	 mv_rev_to(!lc)

   779 	 )

   780     end;

   782 (*. compares two powerproducts .*)

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

   785 (*. help function for mv_add .*)

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

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

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

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

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

   791   | mv_madd(p1,p2,order) =

   792     (

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

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

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

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

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

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

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

   800 	      )

   802 (*. adds two multivariate polynomials .*)

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

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

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

   807 (*. monom multiplication .*)

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

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

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

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

   814 (*. zeros a list .*)

   815 fun mv_null2([])=[]

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

   818 (*. multiplies two multivariate polynomials .*)

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

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

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

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

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

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

   826 fun mv_getmax([])=0

   827   | mv_getmax(x::p1)= let

   828 		       val m=mv_getmax(p1);

   829 		   in

   830 		       if m>x then m

   831 		       else x

   832 		   end;

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

   834 fun mv_grad([]:mv_poly) = 0

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

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

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

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

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

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

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

   845 fun mv_is_negativ([])=false

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

   848 (*. division of monomials .*)

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

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

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

   852     let

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

   854 	val pp=ref [];

   855     in

   856 	(

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

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

   859 	     if #1(p2)<>0 then

   860 		 (

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

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

   863 		  else (!c,!pp)

   864 		      )

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

   866 		 )

   867     end;

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

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

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

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

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

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

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

   878   | mv_division(f,g,order)=

   879     let

   880 	val r=ref [];

   881 	val q=ref [];

   882 	val g'=ref [];

   883 	val k=ref 0;

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

   885 	val exit=ref 0;

   886     in

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

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

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

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

   891 	else

   892 	    (

   893 	     exit:=0;

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

   895 		 (

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

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

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

   899 		      (

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

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

   902 		       )

   903 		  else exit:=1;

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

   905 		  else (exit:=1)

   906 		  );

   907 		 (rev(!q),!r)

   908 		 )

   909     end;

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

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

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

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

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

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

   919 (*. multivariate case .*)

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

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

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

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

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

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

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

   930 fun mv_deg([]:mv_poly) = 0

   931   | mv_deg(p1)=

   932     let

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

   934     in

   935 	if length(p1')=0 then 0

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

   937     end;

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

   940 fun mv_deg2([]:mv_poly) = 0

   941   | mv_deg2(p1)=

   942     let

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

   944     in

   945 	if length(p1')=0 then 0

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

   947     end;

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

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

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

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

   954 fun uv_content2([]:mv_poly) = 0

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

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

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

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

   960     let

   961 	val count=ref 0;

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

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

   964 	val list=ref [];

   965     in

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

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

   968 	     else

   969 		 (

   970 		  count:=0;

   971 		  while (!count)<=max do

   972 		      (

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

   974 			   (

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

   976 			    help:=tl(!help)

   977 			    )

   978 		       else

   979 			   (

   980 			    list:= 0::(!list)

   981 			    );

   982 		       count := (!count) + 1

   983 		       );

   984 		      (!list)

   985 		      )

   986     end;

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

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

   990   | uv_to_poly p1 =

   991     let

   992 	val count=ref 0;

   993 	val help=ref p1;

   994 	val list=ref [];

   995     in

   996 	while length(!help)>0 do

   997 	    (

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

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

  1000 	     count:=(!count)+1;

  1001 	     help:=tl(!help)

  1002 	     );

  1003 	    (!list)

  1004     end;

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

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

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

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

  1010     let

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

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

  1013     in

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

  1015     end

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

  1017     let

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

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

  1020     in

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

  1022     end

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

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

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

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

  1028     let

  1029 	val x=ref (rev(pl));

  1030 	val t=ref [];

  1031 	val y=ref [];

  1032 	val n=ref 1;

  1033 	val n1=ref[];

  1034     in

  1035 	(

  1036 	 while length(!x)>1 do

  1037 	     (

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

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

  1040 		   else n1:=[];

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

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

  1043 	      x:=tl(!x)

  1044 	      );

  1045 	 (!y)

  1046 	 )

  1047     end;

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

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

  1051   | mv_newton_add [x:mv_poly] t = x

  1052   | mv_newton_add (pl:mv_poly list) t =

  1053     let

  1054 	val expos=ref [];

  1055 	val pll=ref pl;

  1056     in

  1057 	(

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

  1060 	     (

  1061 	      pll:=tl(!pll)

  1062 	      );

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

  1064 	 mv_add(hd(pl),

  1065 		mv_mul(

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

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

  1068 		       LEX_

  1069 		       ),

  1070 		LEX_)

  1071 	 )

  1072     end;

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

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

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

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

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

  1079   | mv_newton pl =

  1080     let

  1081 	val c=ref pl;

  1082 	val c1=ref [];

  1083 	val n=length(pl);

  1084 	val k=ref 1;

  1085 	val i=ref n;

  1086 	val ppl=ref [];

  1087     in

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

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

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

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

  1092 	    (

  1093 	     k:=(!k)+1;

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

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

  1096 	     ppl:= !c;

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

  1098 	     );

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

  1100 	c1:=rev(!c1);

  1101 	ppl:= !c1;

  1102 	mv_newton_add (!c1) 1

  1103     end;

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

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

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

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

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

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

  1112     let

  1113 	val y=ref xs;

  1114 	val x=ref [];

  1115     in

  1116 	(

  1117 	 x:=e;

  1118 	 while length(!y)>0 do

  1119 	     (

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

  1121 	      y:=tl(!y)

  1122 	      );

  1123 	 !x

  1124 	 )

  1125     end;

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

  1128 fun list_is_null [] = true

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

  1131 (* check if main variable is zero*)

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

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

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

  1136   | mv_content(p1) =

  1137     let

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

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

  1140 	val result=ref [];

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

  1142     in

  1143 	(

  1144 	 if length(!list)>1 then

  1145 	     (

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

  1147 		  (

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

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

  1151 		   else list:=tl(!list)

  1153 		       );

  1154 		  if length(!list)>0 then

  1155 		   (

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

  1157 		    )

  1158 		  else list:=(!result);

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

  1160 		  (!list)

  1161 		  )

  1162 	 else

  1163 	     (

  1164 	      mv_null3(!list)

  1165 	      )

  1166 	     )

  1167     end

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

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

  1171   | mv_pp(p1) = let

  1172 		    val cont=ref [];

  1173 		    val pp=ref[];

  1174 		in

  1175 		    cont:=mv_content(p1);

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

  1177 		    if !pp=[]

  1178 			then raise error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ")

  1179 		    else (!pp)

  1180 		end

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

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

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

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

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

  1187      let

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

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

  1190      in

  1191 	(

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

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

  1194         )

  1195     end

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

  1197 	(

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

  1199 	)

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

  1201 	(

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

  1203         )

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

  1205     let

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

  1207 	val cont =

  1208 		  (

  1209                    if main_zero(mv_content(p1')) andalso

  1210                      (main_zero(mv_content(p2'))) then

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

  1212                    else

  1213                      mv_gcd (mv_content(p1')) (mv_content(p2'))

  1214                   );

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

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

  1217 	val gcd=ref [];

  1218 	val candidate=ref [];

  1219 	val interpolation_list=ref [];

  1220 	val delta=ref [];

  1221         val p1r = ref [];

  1222         val p2r = ref [];

  1223         val p1r' = ref [];

  1224         val p2r' = ref [];

  1225 	val factor=ref [];

  1226 	val r=ref 0;

  1227 	val gcd_r=ref [];

  1228 	val d=ref 0;

  1229 	val exit=ref 0;

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

  1231     in

  1232 	(

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

  1234 	     (

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

  1236 	      )

  1237 	 else

  1238 	     (

  1239 	      while !exit=0 do

  1240 		  (

  1241 		   r:=(!r)+1;

  1242                    p1r := mv_lc(p1,LEX_);

  1243 		   p2r := mv_lc(p2,LEX_);

  1244                    if main_zero(!p1r) andalso

  1245                       main_zero(!p2r)

  1246                    then

  1247                        (

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

  1249                        )

  1250                    else

  1251                        (

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

  1253                        );

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

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

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

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

  1258                    then

  1259                        (

  1260 		       )

  1261 		   else

  1262 		       (

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

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

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

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

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

  1268 			if (!d < !current_degree) then

  1269 			    (

  1270 			     current_degree:= !d;

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

  1272 			     )

  1273 			else

  1274 			    (

  1275 			     if (!d = !current_degree) then

  1276 				 (

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

  1278 				  )

  1279 			     else ()

  1280 				 )

  1281 			    );

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

  1283 			  (

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

  1285 			   if !candidate=[] then ()

  1286 			   else

  1287 			       (

  1288 				candidate:=mv_pp(!candidate);

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

  1290 				    (

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

  1292 				     exit:=1

  1293 				     )

  1294 				else ()

  1295 				    );

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

  1297 			       )

  1298 		      else ()

  1299 			  )

  1300 	     );

  1301 	     (!gcd):mv_poly

  1302 	     )

  1303     end;

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

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

  1308     (

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

  1310      );

  1312 (*. gets the variables (strings) of a term .*)

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

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

  1316 fun count_neg ([]:mv_poly) = 0

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

  1318 			  else count_neg xs;

  1320 (*. help function for is_polynomial

  1321     checks the order of the operators .*)

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

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

  1324   | test_polynomial (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false

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

  1326   | test_polynomial (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

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

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

  1329   | test_polynomial _ v = false;

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

  1332 fun is_polynomial t = test_polynomial t " ";

  1334 (*. help function for is_expanded

  1335     checks the order of the operators .*)

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

  1337   | test_exp (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false

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

  1339   | test_exp (t as Const ("op +",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

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

  1341   | test_exp (t as Const ("op -",_) $ t1 $ t2) v = if v="*" orelse v="^" then false

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

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

  1344   | test_exp  _ v = false;

  1347 (*. help function for check_coeff:

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

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

  1350     let

  1351 	val x=ref NONE;

  1352 	val len=ref 0;

  1353 	val vl=ref [];

  1354 	val vh=ref [];

  1355 	val i=ref 0;

  1356     in

  1357 	if is_numeral str then

  1358 	    (

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

  1360 		 )

  1361 	else (* variable *)

  1362 	    (

  1363 	     len:=length(v);

  1364 	     vh:=v;

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

  1366 		 (

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

  1368 		      (

  1369 		       vl:=1::(!vl)

  1370 		       )

  1371 		  else

  1372 		      (

  1373 		       vl:=0::(!vl)

  1374 		       );

  1375 		      vh:=tl(!vh);

  1376 		      i:=(!i)+1

  1377 		      );

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

  1379 	    )

  1380     end

  1381   | term2coef' (Const ("op *",_) $ t1 $ t2) v :mv_poly option=

  1382     let

  1383 	val t1pp=ref [];

  1384 	val t2pp=ref [];

  1385 	val t1c=ref 0;

  1386 	val t2c=ref 0;

  1387     in

  1388 	(

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

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

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

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

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

  1396 	 )

  1397     end

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

  1399     let

  1400 	val x=ref NONE;

  1401 	val len=ref 0;

  1402 	val vl=ref [];

  1403 	val vh=ref [];

  1404 	val vtemp=ref [];

  1405 	val i=ref 0;

  1406     in

  1407     (

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

  1409 	 (

  1410 	  len:=length(v);

  1411 	  vh:=v;

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

  1414 	      (

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

  1416 		   (

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

  1418 		    )

  1419 	       else

  1420 		   (

  1421 		    vl:=0::(!vl)

  1422 		    );

  1423 		   vh:=tl(!vh);

  1424 		   i:=(!i)+1

  1425 		   );

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

  1427 	      )

  1428      else raise error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term")

  1429 	 )

  1430     end

  1431   | term2coef' (Const ("op +",_) $ t1 $ t2) v :mv_poly option=

  1432     (

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

  1434 	 )

  1435   | term2coef' (Const ("op -",_) $ t1 $ t2) v :mv_poly option=

  1436     (

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

  1438 	 )

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

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

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

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

  1444     else false;

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

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

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

  1450 fun mk_monom v' p vs =

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

  1452     in map (conv p) vs end;

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

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

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

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

  1458     if is_numeral str

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

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

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

  1463     let

  1464 	val x=ref NONE;

  1465 	val len=ref 0;

  1466 	val vl=ref [];

  1467 	val vh=ref [];

  1468 	val i=ref 0;

  1469     in

  1470 	if is_numeral str then

  1471 	    (

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

  1473 		 )

  1474 	else (* variable *)

  1475 	    (

  1476 	     len:=length v;

  1477 	     vh:= v;

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

  1479 		 (

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

  1481 		      (

  1482 		       vl:=1::(!vl)

  1483 		       )

  1484 		  else

  1485 		      (

  1486 		       vl:=0::(!vl)

  1487 		       );

  1488 		      vh:=tl(!vh);

  1489 		      i:=(!i)+1

  1490 		      );

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

  1492 	    )

  1493     end

  1494   | term2poly' (Const ("op *",_) $ t1 $ t2) v :mv_poly option=

  1495     let

  1496 	val t1pp=ref [];

  1497 	val t2pp=ref [];

  1498 	val t1c=ref 0;

  1499 	val t2c=ref 0;

  1500     in

  1501 	(

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

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

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

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

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

  1508 	 handle _ => NONE

  1510 	 )

  1511     end

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

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

  1514     let

  1515 	val x=ref NONE;

  1516 	val len=ref 0;

  1517 	val vl=ref [];

  1518 	val vh=ref [];

  1519 	val vtemp=ref [];

  1520 	val i=ref 0;

  1521     in

  1522     (

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

  1524 	 (

  1525 	  len:=length(v);

  1526 	  vh:=v;

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

  1529 	      (

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

  1531 		   (

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

  1533 		    )

  1534 	       else

  1535 		   (

  1536 		    vl:=0::(!vl)

  1537 		    );

  1538 		   vh:=tl(!vh);

  1539 		   i:=(!i)+1

  1540 		   );

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

  1542 	      )

  1543      else raise error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term")

  1544 	 )

  1545     end

  1546   | term2poly' (Const ("op +",_) $ t1 $ t2) v :mv_poly option =

  1547     (

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

  1549 	 )

  1550   | term2poly' (Const ("op -",_) $ t1 $ t2) v :mv_poly option =

  1551     (

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

  1553 	 )

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

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

  1557     term2poly

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

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

  1560     			       DIE DU MIR LETZTES MAL GEGEBEN HAST*)

  1561     	 mv_monom list        (*internal representation           *)

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

  1563 .*)

  1564 fun term2poly (t:term) v =

  1565      if is_polynomial t then term2poly' t v

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

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

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

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

  1572 fun expanded2poly (t:term) v =

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

  1574     (*else raise error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*);

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

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

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

  1580 fun powerproduct2term(xs,v) =

  1581     let

  1582 	val xss=ref xs;

  1583 	val vv=ref v;

  1584     in

  1585 	(

  1586 	 while hd(!xss)=0 do

  1587 	     (

  1588 	      xss:=tl(!xss);

  1589 	      vv:=tl(!vv)

  1590 	      );

  1592 	 if list_is_null(tl(!xss)) then

  1593 	     (

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

  1595 	      else

  1596 		  (

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

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

  1599 		   )

  1600 	      )

  1601 	 else

  1602 	     (

  1603 	      if hd(!xss)=1 then

  1604 		  (

  1605 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

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

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

  1608 		   )

  1609 	      else

  1610 		  (

  1611 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1612 		   (

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

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

  1615 		    ) $

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

  1617 		   )

  1618 	      )

  1619 	 )

  1620     end;

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

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

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

  1625     else

  1626 	(

  1627 	 if list_is_null(e) then

  1628 	     (

  1629 	      Free(str_of_int c,HOLogic.realT)

  1630 	      )

  1631 	 else

  1632 	     (

  1633 	      if c=1 then

  1634 		  (

  1635 		   powerproduct2term(e,v)

  1636 		   )

  1637 	      else

  1638 		  (

  1639 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1640 		   Free(str_of_int c,HOLogic.realT)  $

  1641 		   powerproduct2term(e,v)

  1642 		   )

  1643 		  )

  1644 	     );*)

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

  1648     if list_is_null is

  1649     then

  1650 	if i >= 0

  1651 	then Free (str_of_int i, HOLogic.realT)

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

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

  1654     else

  1655 	if i > 0

  1656 	then Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $

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

  1658 		   powerproduct2term(is, v)

  1659 	else Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $

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

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

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

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

  1664     if list_is_null is

  1665     then Free (str_of_int i, HOLogic.realT)

  1666     else if i = 1

  1667     then powerproduct2term (is, vs)

  1668     else Const ("op *", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $

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

  1670 	       powerproduct2term (is, vs);

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

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

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

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

  1676     Const("op +", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $

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

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

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

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

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

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

  1685     else

  1686 	(

  1687 	 if list_is_null(e) then

  1688 	     (

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

  1690 	      )

  1691 	 else

  1692 	     (

  1693 	      if abs(c)=1 then

  1694 		  (

  1695 		   powerproduct2term(e,v)

  1696 		   )

  1697 	      else

  1698 		  (

  1699 		   Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

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

  1701 		   powerproduct2term(e,v)

  1702 		   )

  1703 		  )

  1704 	     );

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

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

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

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

  1710     if c1<0 then

  1711 	Const("op -",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1712 	exp2term'(others,vars) $

  1713 	(

  1714 	 monom2term2((c1,e1),vars)

  1715 	 )

  1716     else

  1717 	Const("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1718 	exp2term'(others,vars) $

  1719 	(

  1720 	 monom2term2((c1,e1),vars)

  1721 	 );

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

  1724     a term in polynomial representation .*)

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

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

  1728 fun polynomial2expanded t =

  1729     (let

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

  1731     in

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

  1733     end) handle _ => NONE;

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

  1736 fun expanded2polynomial t =

  1737     (let

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

  1739     in

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

  1741     end) handle _ => NONE;

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

  1745 fun step_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) =

  1746     let

  1747 	val p1' = ref [];

  1748 	val p2' = ref [];

  1749 	val p3  = ref []

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

  1751     in

  1752 	(

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

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

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

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

  1757 	     (

  1758 	      Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2

  1759 	      )

  1760 	 else

  1761 	     (

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

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

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

  1767 	      (

  1768 	       Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1769 	       $

  1770 	       (

  1771 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1772 		poly2term(!p1',vars) $

  1773 		poly2term(!p3,vars)

  1774 		)

  1775 	       $

  1776 	       (

  1777 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1778 		poly2term(!p2',vars) $

  1779 		poly2term(!p3,vars)

  1780 		)

  1781 	       )

  1782 	      else

  1783 	      (

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

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

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

  1787 	       (

  1788 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1789 		$

  1790 		(

  1791 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1792 		 poly2term(!p1',vars) $

  1793 		 poly2term(!p3,vars)

  1794 		 )

  1795 		$

  1796 		(

  1797 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1798 		 poly2term(!p2',vars) $

  1799 		 poly2term(!p3,vars)

  1800 		 )

  1801 		)

  1802 	       )

  1803 	      )

  1804 	     )

  1805     end

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

  1809 (*. same as step_cancel, this time for expanded forms (input+output) .*)

  1810 fun step_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) =

  1811     let

  1812 	val p1' = ref [];

  1813 	val p2' = ref [];

  1814 	val p3  = ref []

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

  1816     in

  1817 	(

  1818          p1':= sort (mv_geq LEX_) (the (expanded2poly p1 vars ));

  1819        	 p2':= sort (mv_geq LEX_) (the (expanded2poly p2 vars ));

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

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

  1822 	     (

  1823 	      Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2

  1824 	      )

  1825 	 else

  1826 	     (

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

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

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

  1832 	      (

  1833 	       Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1834 	       $

  1835 	       (

  1836 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1837 		poly2expanded(!p1',vars) $

  1838 		poly2expanded(!p3,vars)

  1839 		)

  1840 	       $

  1841 	       (

  1842 		Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1843 		poly2expanded(!p2',vars) $

  1844 		poly2expanded(!p3,vars)

  1845 		)

  1846 	       )

  1847 	      else

  1848 	      (

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

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

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

  1852 	       (

  1853 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1854 		$

  1855 		(

  1856 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1857 		 poly2expanded(!p1',vars) $

  1858 		 poly2expanded(!p3,vars)

  1859 		 )

  1860 		$

  1861 		(

  1862 		 Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  1863 		 poly2expanded(!p2',vars) $

  1864 		 poly2expanded(!p3,vars)

  1865 		 )

  1866 		)

  1867 	       )

  1868 	      )

  1869 	     )

  1870     end

  1871 | step_cancel_expanded _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction");

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

  1874 fun direct_cancel (t as Const ("HOL.divide",_) $ p1 $ p2) =

  1875     let

  1876 	val p1' = ref [];

  1877 	val p2' = ref [];

  1878 	val p3  = ref []

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

  1880     in

  1881 	(

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

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

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

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

  1887 	     (

  1888 	      (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])

  1889 	      )

  1890 	 else

  1891 	     (

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

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

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

  1895 	      (

  1896 	       (

  1897 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1898 		$

  1899 		(

  1900 		 poly2term((!p1'),vars)

  1901 		 )

  1902 		$

  1903 		(

  1904 		 poly2term((!p2'),vars)

  1905 		 )

  1906 		)

  1907 	       ,

  1908 	       if mv_grad(!p3)>0 then

  1909 		   [

  1910 		    (

  1911 		     Const ("Not",[bool]--->bool) $

  1912 		     (

  1913 		      Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $

  1914 		      poly2term((!p3),vars) $

  1915 		      Free("0",HOLogic.realT)

  1916 		      )

  1917 		     )

  1918 		    ]

  1919 	       else

  1920 		   []

  1921 		   )

  1922 	      else

  1923 		  (

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

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

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

  1927 		       (

  1928 			Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1929 			$

  1930 			(

  1931 			 poly2term((!p1'),vars)

  1932 			 )

  1933 			$

  1934 			(

  1935 			 poly2term((!p2'),vars)

  1936 			 )

  1937 			,

  1938 			if mv_grad(!p3)>0 then

  1939 			    [

  1940 			     (

  1941 			      Const ("Not",[bool]--->bool) $

  1942 			      (

  1943 			       Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $

  1944 			       poly2term((!p3),vars) $

  1945 			       Free("0",HOLogic.realT)

  1946 			       )

  1947 			      )

  1948 			     ]

  1949 			else

  1950 			    []

  1951 			    )

  1952 		       )

  1953 		  )

  1954 	     )

  1955     end

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

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

  1959 fun direct_cancel_expanded (t as Const ("HOL.divide",_) $ p1 $ p2) =

  1960     let

  1961 	val p1' = ref [];

  1962 	val p2' = ref [];

  1963 	val p3  = ref []

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

  1965     in

  1966 	(

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

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

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

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

  1972 	     (

  1973 	      (Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[])

  1974 	      )

  1975 	 else

  1976 	     (

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

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

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

  1980 	      (

  1981 	       (

  1982 		Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  1983 		$

  1984 		(

  1985 		 poly2expanded((!p1'),vars)

  1986 		 )

  1987 		$

  1988 		(

  1989 		 poly2expanded((!p2'),vars)

  1990 		 )

  1991 		)

  1992 	       ,

  1993 	       if mv_grad(!p3)>0 then

  1994 		   [

  1995 		    (

  1996 		     Const ("Not",[bool]--->bool) $

  1997 		     (

  1998 		      Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $

  1999 		      poly2expanded((!p3),vars) $

  2000 		      Free("0",HOLogic.realT)

  2001 		      )

  2002 		     )

  2003 		    ]

  2004 	       else

  2005 		   []

  2006 		   )

  2007 	      else

  2008 		  (

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

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

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

  2012 		       (

  2013 			Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2014 			$

  2015 			(

  2016 			 poly2expanded((!p1'),vars)

  2017 			 )

  2018 			$

  2019 			(

  2020 			 poly2expanded((!p2'),vars)

  2021 			 )

  2022 			,

  2023 			if mv_grad(!p3)>0 then

  2024 			    [

  2025 			     (

  2026 			      Const ("Not",[bool]--->bool) $

  2027 			      (

  2028 			       Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $

  2029 			       poly2expanded((!p3),vars) $

  2030 			       Free("0",HOLogic.realT)

  2031 			       )

  2032 			      )

  2033 			     ]

  2034 			else

  2035 			    []

  2036 			    )

  2037 		       )

  2038 		  )

  2039 	     )

  2040     end

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

  2044 (*. adds two fractions .*)

  2045 fun add_fract ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =

  2046     let

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

  2048 	val t11'=ref (the(term2poly t11 vars));

  2049 val _= writeln"### add_fract: done t11"

  2050 	val t12'=ref (the(term2poly t12 vars));

  2051 val _= writeln"### add_fract: done t12"

  2052 	val t21'=ref (the(term2poly t21 vars));

  2053 val _= writeln"### add_fract: done t21"

  2054 	val t22'=ref (the(term2poly t22 vars));

  2055 val _= writeln"### add_fract: done t22"

  2056 	val den=ref [];

  2057 	val nom=ref [];

  2058 	val m1=ref [];

  2059 	val m2=ref [];

  2060     in

  2062 	(

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

  2064 writeln"### add_fract: done sort mv_lcm";

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

  2066 writeln"### add_fract: done sort mv_division t12";

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

  2068 writeln"### add_fract: done sort mv_division t22";

  2069 	 nom :=sort (mv_geq LEX_)

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

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

  2072 				       LEX_),

  2073 				LEX_));

  2074 writeln"### add_fract: done sort mv_add";

  2075 	 (

  2076 	  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2077 	  $

  2078 	  (

  2079 	   poly2term((!nom),vars)

  2080 	   )

  2081 	  $

  2082 	  (

  2083 	   poly2term((!den),vars)

  2084 	   )

  2085 	  )

  2086 	 )

  2087     end

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

  2090 (*. adds two expanded fractions .*)

  2091 fun add_fract_exp ((Const("HOL.divide",_) $ t11 $ t12),(Const("HOL.divide",_) $ t21 $ t22)) =

  2092     let

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

  2094 	val t11'=ref (the(expanded2poly t11 vars));

  2095 	val t12'=ref (the(expanded2poly t12 vars));

  2096 	val t21'=ref (the(expanded2poly t21 vars));

  2097 	val t22'=ref (the(expanded2poly t22 vars));

  2098 	val den=ref [];

  2099 	val nom=ref [];

  2100 	val m1=ref [];

  2101 	val m2=ref [];

  2102     in

  2104 	(

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

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

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

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

  2109 	 (

  2110 	  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2111 	  $

  2112 	  (

  2113 	   poly2expanded((!nom),vars)

  2114 	   )

  2115 	  $

  2116 	  (

  2117 	   poly2expanded((!den),vars)

  2118 	   )

  2119 	  )

  2120 	 )

  2121     end

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

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

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

  2126   | add_list_of_fractions [x]= direct_cancel x

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

  2128     let

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

  2130 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(x)";

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

  2132 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(y)";

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

  2134 val _= writeln"### add_list_of_fractions xs: has done add_list_of_fraction xs";

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

  2136 val _= writeln"### add_list_of_fractions xs: has done direct_cancel(t3a)";

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

  2138     in

  2139 	(writeln"### add_list_of_fractions in";

  2140 	 (

  2141 	 (t4a,rest)

  2142 	 )

  2143 	 )

  2144     end;

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

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

  2148   | add_list_of_fractions_exp [x]= direct_cancel_expanded x

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

  2150     let

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

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

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

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

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

  2156     in

  2157 	(

  2158 	 (t4a,rest)

  2159 	 )

  2160     end;

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

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

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

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

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

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

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

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

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

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

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

  2176   | termlist2denominators xs =

  2177     let

  2178 	val xxs=ref xs;

  2179 	val var=ref [];

  2180     in

  2181 	var:=[];

  2182 	while length(!xxs)>0 do

  2183 	    (

  2184 	     let

  2185 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);

  2186 	     in

  2187 		 (

  2188 		  xxs:=tl(!xxs);

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

  2190 		  )

  2191 	     end

  2192 	     );

  2193 	    (t2d(xs,!var),!var)

  2194     end;

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

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

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

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

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

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

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

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

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

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

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

  2210   | termlist2denominators xs =

  2211     let

  2212 	val xxs=ref xs;

  2213 	val var=ref [];

  2214     in

  2215 	var:=[];

  2216 	while length(!xxs)>0 do

  2217 	    (

  2218 	     let

  2219 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);

  2220 	     in

  2221 		 (

  2222 		  xxs:=tl(!xxs);

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

  2224 		  )

  2225 	     end

  2226 	     );

  2227 	    (t2d(xs,!var),!var)

  2228     end;

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

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

  2232   | termlist2denominators_exp xs =

  2233     let

  2234 	val xxs=ref xs;

  2235 	val var=ref [];

  2236     in

  2237 	var:=[];

  2238 	while length(!xxs)>0 do

  2239 	    (

  2240 	     let

  2241 		 val (t as Const ("HOL.divide",_) $ p1x $ p2x)=hd(!xxs);

  2242 	     in

  2243 		 (

  2244 		  xxs:=tl(!xxs);

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

  2246 		  )

  2247 	     end

  2248 	     );

  2249 	    (t2d_exp(xs,!var),!var)

  2250     end;

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

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

  2254     let

  2255 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;

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

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

  2258 	val p1var=get_vars(p1');

  2259     in

  2260 	if length(xs)>0 then

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

  2262 		(

  2263 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2264 		 $

  2265 		 (

  2266 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2267 		  $

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

  2269 		  $

  2270 		  den

  2271 		  )

  2272 		 $

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

  2274 		,

  2275 		[]

  2276 		)

  2277 	    else

  2278 		(

  2279 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2280 		 $

  2281 		 (

  2282 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2283 		  $

  2284 		  (

  2285 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2286 		   poly2term(the (term2poly p1' p1var),p1var) $

  2287 		   poly2term(p3,var)

  2288 		   )

  2289 		  $

  2290 		  (

  2291 		   den

  2292 		   )

  2293 		  )

  2294 		 $

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

  2296 		,

  2297 		[]

  2298 		)

  2299 	else

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

  2301 		(

  2302 		 (

  2303 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2304 		  $

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

  2306 		  $

  2307 		  den

  2308 		  )

  2309 		 ,

  2310 		 []

  2311 		 )

  2312 	     else

  2313 		 (

  2314 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2315 		  $

  2316 		  (

  2317 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2318 		   poly2term(the (term2poly p1' p1var),p1var) $

  2319 		   poly2term(p3,var)

  2320 		   )

  2321 		  $

  2322 		  den

  2323 		  ,

  2324 		  []

  2325 		  )

  2326     end;

  2328 (*. same as com_den, this time for expanded forms .*)

  2329 fun com_den_exp(x::xs,denom,den,var)=

  2330     let

  2331 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;

  2332 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

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

  2334 	val p1var=get_vars(p1');

  2335     in

  2336 	if length(xs)>0 then

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

  2338 		(

  2339 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2340 		 $

  2341 		 (

  2342 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2343 		  $

  2344 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2345 		  $

  2346 		  den

  2347 		  )

  2348 		 $

  2349 		 #1(com_den_exp(xs,denom,den,var))

  2350 		,

  2351 		[]

  2352 		)

  2353 	    else

  2354 		(

  2355 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2356 		 $

  2357 		 (

  2358 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2359 		  $

  2360 		  (

  2361 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2362 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2363 		   poly2expanded(p3,var)

  2364 		   )

  2365 		  $

  2366 		  (

  2367 		   den

  2368 		   )

  2369 		  )

  2370 		 $

  2371 		 #1(com_den_exp(xs,denom,den,var))

  2372 		,

  2373 		[]

  2374 		)

  2375 	else

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

  2377 		(

  2378 		 (

  2379 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2380 		  $

  2381 		  poly2expanded(the(expanded2poly p1' p1var),p1var)

  2382 		  $

  2383 		  den

  2384 		  )

  2385 		 ,

  2386 		 []

  2387 		 )

  2388 	     else

  2389 		 (

  2390 		  Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT)

  2391 		  $

  2392 		  (

  2393 		   Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2394 		   poly2expanded(the(expanded2poly p1' p1var),p1var) $

  2395 		   poly2expanded(p3,var)

  2396 		   )

  2397 		  $

  2398 		  den

  2399 		  ,

  2400 		  []

  2401 		  )

  2402     end;

  2404 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2405 -------------------------------------------------------------

  2406 (* WN0210???SK brauch ma des überhaupt *)

  2407 fun com_den2(x::xs,denom,den,var)=

  2408     let

  2409 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;

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

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

  2412 	val p1var=get_vars(p1');

  2413     in

  2414 	if length(xs)>0 then

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

  2416 		(

  2417 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2418 		 poly2term(the(term2poly p1' p1var),p1var) $

  2419 		 com_den2(xs,denom,den,var)

  2420 		)

  2421 	    else

  2422 		(

  2423 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2424 		 (

  2425 		   let

  2426 		       val p3'=poly2term(p3,var);

  2427 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2428 		   in

  2429 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2430 		   end

  2431 		  ) $

  2432 		 com_den2(xs,denom,den,var)

  2433 		)

  2434 	else

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

  2436 		(

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

  2438 		 )

  2439 	     else

  2440 		 (

  2441 		   let

  2442 		       val p3'=poly2term(p3,var);

  2443 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2444 		   in

  2445 		       poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars)

  2446 		   end

  2447 		  )

  2448     end;

  2450 (* WN0210???SK brauch ma des überhaupt *)

  2451 fun com_den_exp2(x::xs,denom,den,var)=

  2452     let

  2453 	val (t as Const ("HOL.divide",_) $ p1' $ p2')=x;

  2454 	val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var));

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

  2456 	val p1var=get_vars p1';

  2457     in

  2458 	if length(xs)>0 then

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

  2460 		(

  2461 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2462 		 poly2expanded(the (expanded2poly p1' p1var),p1var) $

  2463 		 com_den_exp2(xs,denom,den,var)

  2464 		)

  2465 	    else

  2466 		(

  2467 		 Const ("op +",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2468 		 (

  2469 		   let

  2470 		       val p3'=poly2expanded(p3,var);

  2471 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2472 		   in

  2473 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2474 		   end

  2475 		  ) $

  2476 		 com_den_exp2(xs,denom,den,var)

  2477 		)

  2478 	else

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

  2480 		(

  2481 		 poly2expanded(the (expanded2poly p1' p1var),p1var)

  2482 		 )

  2483 	     else

  2484 		 (

  2485 		   let

  2486 		       val p3'=poly2expanded(p3,var);

  2487 		       val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3');

  2488 		   in

  2489 		       poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars)

  2490 		   end

  2491 		  )

  2492     end;

  2493 ---------------------------------------------------------*)

  2496 (*. searches for an element y of a list ys, which has an gcd not 1 with x .*)

  2497 fun exists_gcd (x,[]) = false

  2498   | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then  exists_gcd (x,ys)

  2499 			   else true;

  2501 (*. divides each element of the list xs with y .*)

  2502 fun list_div ([],y) = []

  2503   | list_div (x::xs,y) =

  2504     let

  2505 	val (d,r)=mv_division(x,y,LEX_);

  2506     in

  2507 	if r=[] then

  2508 	    d::list_div(xs,y)

  2509 	else x::list_div(xs,y)

  2510     end;

  2512 (*. checks if x is in the list ys .*)

  2513 fun in_list (x,[]) = false

  2514   | in_list (x,y::ys) = if x=y then true

  2515 			else in_list(x,ys);

  2517 (*. deletes all equal elements of the list xs .*)

  2518 fun kill_equal [] = []

  2519   | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs)

  2520 			 else x::kill_equal(xs);

  2522 (*. searches for new factors .*)

  2523 fun new_factors [] = []

  2524   | new_factors (list:mv_poly list):mv_poly list =

  2525     let

  2526 	val l = kill_equal list;

  2527 	val len = length(l);

  2528     in

  2529 	if len>=2 then

  2530 	    (

  2531 	     let

  2532 		 val x::y::xs=l;

  2533 		 val gcd=mv_gcd x y;

  2534 	     in

  2535 		 if gcd=[(1,mv_null2(#2(hd(x))))] then

  2536 		     (

  2537 		      if exists_gcd(x,xs) then new_factors (y::xs @ [x])

  2538 		      else x::new_factors(y::xs)

  2539 	             )

  2540 		 else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd)))

  2541 	     end

  2542 	     )

  2543 	else

  2544 	    if len=1 then [hd(l)]

  2545 	    else []

  2546     end;

  2548 (*. gets the factors of a list .*)

  2549 fun get_factors x = new_factors x;

  2551 (*. multiplies the elements of the list .*)

  2552 fun multi_list [] = []

  2553   | multi_list (x::xs) = if xs=[] then x

  2554 			 else mv_mul(x,multi_list xs,LEX_);

  2556 (*. makes a term out of the elements of the list (polynomial representation) .*)

  2557 fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2558   | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars)

  2559 			       else

  2560 				   (

  2561 				    Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2562 				    poly2term(sort (mv_geq LEX_) (x),vars) $

  2563 				    make_term(xs,vars)

  2564 				    );

  2566 (*. factorizes the denominator (polynomial representation) .*)

  2567 fun factorize_den (l,den,vars) =

  2568     let

  2569 	val factor_list=kill_equal( (get_factors l));

  2570 	val mlist=multi_list(factor_list);

  2571 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2572     in

  2573 	if rest=[] then

  2574 	    (

  2575 	     if last=[(1,mv_null2(vars))] then make_term(factor_list,vars)

  2576 	     else make_term(last::factor_list,vars)

  2577 	     )

  2578 	else raise error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division")

  2579     end;

  2581 (*. makes a term out of the elements of the list (expanded polynomial representation) .*)

  2582 fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT)

  2583   | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars)

  2584 			       else

  2585 				   (

  2586 				    Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2587 				    poly2expanded(sort (mv_geq LEX_) (x),vars) $

  2588 				    make_exp(xs,vars)

  2589 				    );

  2591 (*. factorizes the denominator (expanded polynomial representation) .*)

  2592 fun factorize_den_exp (l,den,vars) =

  2593     let

  2594 	val factor_list=kill_equal( (get_factors l));

  2595 	val mlist=multi_list(factor_list);

  2596 	val (last,rest)=mv_division(den,multi_list(factor_list),LEX_);

  2597     in

  2598 	if rest=[] then

  2599 	    (

  2600 	     if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars)

  2601 	     else make_exp(last::factor_list,vars)

  2602 	     )

  2603 	else raise error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division")

  2604     end;

  2606 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2607 (*. the nominators and denominators with the correct factor .*)

  2608 (*. (polynomial representation) .*)

  2609 fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list)

  2610   | step_add_list_of_fractions [x]= raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add")

  2611   | step_add_list_of_fractions (xs) =

  2612     let

  2613         val den_list=termlist2denominators (xs); (* list of denominators *)

  2614 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2615 	val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2616     in

  2617 	com_den(xs,denom,den,var)

  2618     end;

  2620 (*. calculates the common denominator of all elements of the list and multiplies .*)

  2621 (*. the nominators and denominators with the correct factor .*)

  2622 (*. (expanded polynomial representation) .*)

  2623 fun step_add_list_of_fractions_exp []  = (Free("0",HOLogic.realT),[]:term list)

  2624   | step_add_list_of_fractions_exp [x] = raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add")

  2625   | step_add_list_of_fractions_exp (xs)=

  2626     let

  2627         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2628 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2629 	val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *)

  2630     in

  2631 	com_den_exp(xs,denom,den,var)

  2632     end;

  2634 (* wird aktuell nicht mehr gebraucht, bei rückänderung schon

  2635 -------------------------------------------------------------

  2636 (* WN0210???SK brauch ma des überhaupt *)

  2637 fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list)

  2638   | step_add_list_of_fractions2 [x]=(x,[])

  2639   | step_add_list_of_fractions2 (xs) =

  2640     let

  2641         val den_list=termlist2denominators (xs); (* list of denominators *)

  2642 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2643 	val den=factorize_den(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2644     in

  2645 	(

  2646 	 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2647 	 com_den2(xs,denom, poly2term(denom,var)(*den*),var) $

  2648 	 poly2term(denom,var)

  2649 	,

  2650 	[]

  2651 	)

  2652     end;

  2654 (* WN0210???SK brauch ma des überhaupt *)

  2655 fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list)

  2656   | step_add_list_of_fractions2_exp [x]=(x,[])

  2657   | step_add_list_of_fractions2_exp (xs) =

  2658     let

  2659         val den_list=termlist2denominators_exp (xs); (* list of denominators *)

  2660 	val (denom,var)=calc_lcm(den_list);      (* common denominator *)

  2661 	val den=factorize_den_exp(#1(den_list),denom,var);  (* faktorisierter Nenner !!! *)

  2662     in

  2663 	(

  2664 	 Const ("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2665 	 com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $

  2666 	 poly2expanded(denom,var)

  2667 	,

  2668 	[]

  2669 	)

  2670     end;

  2671 ---------------------------------------------- *)

  2674 (*. converts a term, which contains severel terms seperated by +, into a list of these terms .*)

  2675 fun term2list (t as (Const("HOL.divide",_) $ _ $ _)) = [t]

  2676   | term2list (t as (Const("Atools.pow",_) $ _ $ _)) =

  2677     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2678 	  t $ Free("1",HOLogic.realT)

  2679      ]

  2680   | term2list (t as (Free(_,_))) =

  2681     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2682 	  t $  Free("1",HOLogic.realT)

  2683      ]

  2684   | term2list (t as (Const("op *",_) $ _ $ _)) =

  2685     [Const("HOL.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $

  2686 	  t $ Free("1",HOLogic.realT)

  2687      ]

  2688   | term2list (Const("op +",_) $ t1 $ t2) = term2list(t1) @ term2list(t2)

  2689   | term2list (Const("op -",_) $ t1 $ t2) =

  2690     raise error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet")

  2691   | term2list _ = raise error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term");

  2693 (*.factors out the gcd of nominator and denominator:

  2694    a/b = (a' * gcd)/(b' * gcd),  a,b,gcd  are poly[2].*)

  2695 fun factout_p_  (thy:theory) t = SOME (step_cancel t,[]:term list);

  2696 fun factout_ (thy:theory) t = SOME (step_cancel_expanded t,[]:term list);

  2698 (*.cancels a single fraction with normalform [2]

  2699    resulting in a canceled fraction [2], see factout_ .*)

  2700 fun cancel_p_ (thy:theory) t = (*WN.2.6.03 no rewrite -> NONE !*)

  2701     (let val (t',asm) = direct_cancel(*_expanded ... corrected MG.21.8.03*) t

  2702      in if t = t' then NONE else SOME (t',asm)

  2703      end) handle _ => NONE;

  2704 (*.the same as above with normalform [3]

  2705   val cancel_ :

  2706       theory ->        (*10.02 unused                                    *)

  2707       term -> 	       (*fraction in normalform [3]                      *)

  2708       (term * 	       (*fraction in normalform [3]                      *)

  2709        term list)      (*casual asumptions in normalform [3]             *)

  2710 	  option       (*NONE: the function is not applicable            *).*)

  2711 fun cancel_ (thy:theory) t = SOME (direct_cancel_expanded t) handle _ => NONE;

  2713 (*.transforms sums of at least 2 fractions [3] to

  2714    sums with the least common multiple as nominator.*)

  2715 fun common_nominator_p_ (thy:theory) t =

  2716 ((*writeln("### common_nominator_p_ called");*)

  2717     SOME (step_add_list_of_fractions(term2list(t))) handle _ => NONE

  2718 );

  2719 fun common_nominator_ (thy:theory) t =

  2720     SOME (step_add_list_of_fractions_exp(term2list(t))) handle _ => NONE;

  2722 (*.add 2 or more fractions

  2723 val add_fraction_p_ :

  2724       theory ->        (*10.02 unused                                    *)

  2725       term -> 	       (*2 or more fractions with normalform [2]         *)

  2726       (term * 	       (*one fraction with normalform [2]                *)

  2727        term list)      (*casual assumptions in normalform [2] WN0210???SK  *)

  2728 	  option       (*NONE: the function is not applicable            *).*)

  2729 fun add_fraction_p_ (thy:theory) t =

  2730 (writeln("### add_fraction_p_ called");

  2731     (let val ts = term2list t

  2732      in if 1 < length ts

  2733 	then SOME (add_list_of_fractions ts)

  2734 	else NONE (*raise error ("RATIONALS_ADD_EXCEPTION: nothing to add")*)

  2735      end) handle _ => NONE

  2736 );

  2737 (*.same as add_fraction_p_ but with normalform [3].*)

  2738 (*SOME (step_add_list_of_fractions2(term2list(t))); *)

  2739 fun add_fraction_ (thy:theory) t =

  2740     if length(term2list(t))>1

  2741     then SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE

  2742     else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)

  2743 	NONE;

  2744 fun add_fraction_ (thy:theory) t =

  2745     (if 1 < length (term2list t)

  2746      then SOME (add_list_of_fractions_exp (term2list t))

  2747      else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*)

  2748 	 NONE) handle _ => NONE;

  2750 (*SOME (step_add_list_of_fractions2_exp(term2list(t))); *)

  2752 (*. brings the term into a normal form .*)

  2753 fun norm_rational_ (thy:theory) t =

  2754     SOME (add_list_of_fractions(term2list(t))) handle _ => NONE;

  2755 fun norm_expanded_rat_ (thy:theory) t =

  2756     SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE;

  2759 (*.evaluates conditions in calculate_Rational.*)

  2760 (*make local with FIXX@ME result:term *term list*)

  2761 val calc_rat_erls = prep_rls(

  2762   Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  2763 	 erls = e_rls, srls = Erls, calc = [], (*asm_thm = [], *)

  2764 	 rules =

  2765 	 [Calc ("op =",eval_equal "#equal_"),

  2766 	  Calc ("Atools.is'_const",eval_const "#is_const_"),

  2767 	  Thm ("not_true",num_str not_true),

  2768 	  Thm ("not_false",num_str not_false)

  2769 	  ],

  2770 	 scr = EmptyScr});

  2773 (*.simplifies expressions with numerals;

  2774    does NOT rearrange the term by AC-rewriting; thus terms with variables

  2775    need to have constants to be commuted together respectively.*)

  2776 val calculate_Rational = prep_rls(

  2777     merge_rls "calculate_Rational"

  2778 	(Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  2779 	      erls = calc_rat_erls, srls = Erls, (*asm_thm = [],*)

  2780 	      calc = [],

  2781 	      rules =

  2782 	      [Calc ("HOL.divide"  ,eval_cancel "#divide_"),

  2784 	       Thm ("sym_real_minus_divide_eq",

  2785 		    num_str (real_minus_divide_eq RS sym)),

  2786 	       (*SYM - ?x / ?y = - (?x / ?y)  may come from subst*)

  2788 	       Thm ("rat_add",num_str rat_add),

  2789 	       (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \

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

  2791 	       Thm ("rat_add1",num_str rat_add1),

  2792 	       (*"[| a is_const; b is_const; c is_const |] ==> \

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

  2794 	       Thm ("rat_add2",num_str rat_add2),

  2795 	       (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> \

  2796 		 \?a / ?c + ?b = (?a + ?b * ?c) / ?c"*)

  2797 	       Thm ("rat_add3",num_str rat_add3),

  2798 	       (*"[| a is_const; b is_const; c is_const |] ==> \

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

  2800 		 \.... is_const to be omitted here FIXME*)

  2802 	       Thm ("rat_mult",num_str rat_mult),

  2803 	       (*a / b * (c / d) = a * c / (b * d)*)

  2804 	       Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),

  2805 	       (*?x * (?y / ?z) = ?x * ?y / ?z*)

  2806 	       Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),

  2807 	       (*?y / ?z * ?x = ?y * ?x / ?z*)

  2809 	       Thm ("real_divide_divide1",num_str real_divide_divide1),

  2810 	       (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*)

  2811 	       Thm ("real_divide_divide2_eq",num_str real_divide_divide2_eq),

  2812 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  2814 	       Thm ("rat_power", num_str rat_power),

  2815 	       (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

  2817 	       Thm ("mult_cross",num_str mult_cross),

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

  2819 	       Thm ("mult_cross1",num_str mult_cross1),

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

  2821 	       Thm ("mult_cross2",num_str mult_cross2)

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

  2823 	       ], scr = EmptyScr})

  2824 	calculate_Poly);

  2827 (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*)

  2828 fun eval_is_expanded (thmid:string) _

  2829 		       (t as (Const("Rational.is'_expanded", _) $ arg)) thy =

  2830     if is_expanded arg

  2831     then SOME (mk_thmid thmid ""

  2832 			((Syntax.string_of_term (thy2ctxt thy)) arg) "",

  2833 	       Trueprop $ (mk_equality (t, HOLogic.true_const)))

  2834     else SOME (mk_thmid thmid ""

  2835 			((Syntax.string_of_term (thy2ctxt thy)) arg) "",

  2836 	       Trueprop $ (mk_equality (t, HOLogic.false_const)))

  2837   | eval_is_expanded _ _ _ _ = NONE;

  2839 val rational_erls =

  2840     merge_rls "rational_erls" calculate_Rational

  2841 	      (append_rls "is_expanded" Atools_erls

  2842 			  [Calc ("Rational.is'_expanded", eval_is_expanded "")

  2843 			   ]);

  2847 (*.3 'reverse-rewrite-sets' for symbolic computation on rationals:

  2848  =================================================================

  2849  A[2] 'cancel_p': .

  2850  A[3] 'cancel': .

  2851  B[2] 'common_nominator_p': transforms summands in a term [2]

  2852          to fractions with the (least) common multiple as nominator.

  2853  B[3] 'norm_rational': normalizes arbitrary algebraic terms (without

  2854          radicals and transzendental functions) to one canceled fraction,

  2855 	 nominator and denominator in polynomial form.

  2857 In order to meet isac's requirements for interactive and stepwise calculation,

  2858 each 'reverse-rewerite-set' consists of an initialization for the interpreter

  2859 state and of 4 functions, each of which employs rewriting as much as possible.

  2860 The signature of these functions are the same in each 'reverse-rewrite-set'

  2861 respectively.*)

  2863 (* ************************************************************************* *)

  2866 local(*. cancel_p

  2867 ------------------------

  2868 cancels a single fraction consisting of two (uni- or multivariate)

  2869 polynomials WN0609???SK[2] into another such a fraction; examples:

  2871 	   a^2 + -1*b^2         a + b

  2872         -------------------- = ---------

  2873 	a^2 + -2*a*b + b^2     a + -1*b

  2875         a^2    a

  2876         --- = ---

  2877          a     1

  2879 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)

  2880 (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)

  2882 val {rules, rew_ord=(_,ro),...} =

  2883     rep_rls (assoc_rls "make_polynomial");

  2884 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  2885            see rational.sml --- investigate rulesets for cancel_p ---*)

  2886 val {rules, rew_ord=(_,ro),...} =

  2887     rep_rls (assoc_rls "rev_rew_p");

  2889 val thy = Rational.thy;

  2891 (*.init_state = fn : term -> istate

  2892 initialzies the state of the script interpreter. The state is:

  2894 type rrlsstate =      (*state for reverse rewriting*)

  2895      (term *          (*the current formula*)

  2896       term *          (*the final term*)

  2897       rule list       (*'reverse rule list' (#)*)

  2898 	    list *    (*may be serveral, eg. in norm_rational*)

  2899       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  2900        (term *        (*... rewrite with ...*)

  2901 	term list))   (*... assumptions*)

  2902 	  list);      (*derivation from given term to normalform

  2903 		       in reverse order with sym_thm;

  2904                        (#) could be extracted from here by (map #1)*).*)

  2905 (* val {rules, rew_ord=(_,ro),...} =

  2906        rep_rls (assoc_rls "rev_rew_p")        (*USE ALWAYS, SEE val cancel_p*);

  2907    val (thy, eval_rls, ro) =(Rational.thy, Atools_erls, ro) (*..val cancel_p*);

  2908    val t = t;

  2909    *)

  2910 fun init_state thy eval_rls ro t =

  2911     let val SOME (t',_) = factout_p_ thy t

  2912         val SOME (t'',asm) = cancel_p_ thy t

  2913         val der = reverse_deriv thy eval_rls rules ro NONE t'

  2914         val der = der @ [(Thm ("real_mult_div_cancel2",

  2915 			       num_str real_mult_div_cancel2),

  2916 			  (t'',asm))]

  2917         val rs = (distinct_Thm o (map #1)) der

  2918 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  2919 				      "sym_real_mult_0",

  2920 				      "sym_real_mult_1"

  2921 				      (*..insufficient,eg.make_Polynomial*)])rs

  2922     in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end;

  2924 (*.locate_rule = fn : rule list -> term -> rule

  2925 		      -> (rule * (term * term list) option) list.

  2926   checks a rule R for being a cancel-rule, and if it is,

  2927   then return the list of rules (+ the terms they are rewriting to)

  2928   which need to be applied before R should be applied.

  2929   precondition: the rule is applicable to the argument-term.

  2930 arguments:

  2931   rule list: the reverse rule list

  2932   -> term  : ... to which the rule shall be applied

  2933   -> rule  : ... to be applied to term

  2934 value:

  2935   -> (rule           : a rule rewriting to ...

  2936       * (term        : ... the resulting term ...

  2937          * term list): ... with the assumptions ( //#0).

  2938       ) list         : there may be several such rules;

  2939 		       the list is empty, if the rule has nothing to do

  2940 		       with cancelation.*)

  2941 (* val () = ();

  2942    *)

  2943 fun locate_rule thy eval_rls ro [rs] t r =

  2944     if (id_of_thm r) mem (map (id_of_thm)) rs

  2945     then let val ropt =

  2946 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  2947 	 in case ropt of

  2948 		SOME ta => [(r, ta)]

  2949 	      | NONE => (writeln("### locate_rule:  rewrite "^

  2950 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  2951 			 []) end

  2952     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  2953   | locate_rule _ _ _ _ _ _ =

  2954     raise error ("locate_rule: doesnt match rev-sets in istate");

  2956 (*.next_rule = fn : rule list -> term -> rule option

  2957   for a given term return the next rules to be done for cancelling.

  2958 arguments:

  2959   rule list     : the reverse rule list

  2960   term          : the term for which ...

  2961 value:

  2962   -> rule option: ... this rule is appropriate for cancellation;

  2963 		  there may be no such rule (if the term is canceled already.*)

  2964 (* val thy = Rational.thy;

  2965    val Rrls {rew_ord=(_,ro),...} = cancel;

  2966    val ([rs],t) = (rss,f);

  2967    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  2969    val (thy, [rs]) = (Rational.thy, revsets);

  2970    val Rrls {rew_ord=(_,ro),...} = cancel;

  2971    nex [rs] t;

  2972    *)

  2973 fun next_rule thy eval_rls ro [rs] t =

  2974     let val der = make_deriv thy eval_rls rs ro NONE t;

  2975     in case der of

  2976 (* val (_,r,_)::_ = der;

  2977    *)

  2978 	   (_,r,_)::_ => SOME r

  2979 	 | _ => NONE

  2980     end

  2981   | next_rule _ _ _ _ _ =

  2982     raise error ("next_rule: doesnt match rev-sets in istate");

  2984 (*.val attach_form = f : rule list -> term -> term

  2985 			 -> (rule * (term * term list)) list

  2986   checks an input term TI, if it may belong to a current cancellation, by

  2987   trying to derive it from the given term TG.

  2988 arguments:

  2989   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  2990   -> term: TI, the next one input by the user

  2991 value:

  2992   -> (rule           : the rule to be applied in order to reach TI

  2993       * (term        : ... obtained by applying the rule ...

  2994          * term list): ... and the respective assumptions.

  2995       ) list         : there may be several such rules;

  2996                        the list is empty, if the users term does not belong

  2997 		       to a cancellation of the term last agreed upon.*)

  2998 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  2999     []:(rule * (term * term list)) list;

  3001 in

  3003 val cancel_p =

  3004     Rrls {id = "cancel_p", prepat=[],

  3005 	  rew_ord=("ord_make_polynomial",

  3006 		   ord_make_polynomial false Rational.thy),

  3007 	  erls = rational_erls,

  3008 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),

  3009 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),

  3010 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),

  3011 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3012 	  (*asm_thm=[("real_mult_div_cancel2","")],*)

  3013 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3014 		     normal_form = cancel_p_ thy,

  3015 		     locate_rule = locate_rule thy Atools_erls ro,

  3016 		     next_rule   = next_rule thy Atools_erls ro,

  3017 		     attach_form = attach_form}}

  3018 end;(*local*)

  3021 local(*.ad (1) 'cancel'

  3022 ------------------------------

  3023 cancels a single fraction consisting of two (uni- or multivariate)

  3024 polynomials WN0609???SK[3] into another such a fraction; examples:

  3026 	   a^2 - b^2           a + b

  3027         -------------------- = ---------

  3028 	a^2 - 2*a*b + b^2      a - *b

  3030 Remark: the reverse ruleset does _NOT_ work properly with other input !.*)

  3031 (*WN 24.8.02: wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*)

  3033 (*

  3034 val SOME (Rls {rules=rules,rew_ord=(_,ro),...}) =

  3035     assoc'(!ruleset',"expand_binoms");

  3036 *)

  3037 val {rules=rules,rew_ord=(_,ro),...} =

  3038     rep_rls (assoc_rls "expand_binoms");

  3039 val thy = Rational.thy;

  3041 fun init_state thy eval_rls ro t =

  3042     let val SOME (t',_) = factout_ thy t;

  3043         val SOME (t'',asm) = cancel_ thy t;

  3044         val der = reverse_deriv thy eval_rls rules ro NONE t';

  3045         val der = der @ [(Thm ("real_mult_div_cancel2",

  3046 			       num_str real_mult_div_cancel2),

  3047 			  (t'',asm))]

  3048         val rs = map #1 der;

  3049     in (t,t'',[rs],der) end;

  3051 fun locate_rule thy eval_rls ro [rs] t r =

  3052     if (id_of_thm r) mem (map (id_of_thm)) rs

  3053     then let val ropt =

  3054 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3055 	 in case ropt of

  3056 		SOME ta => [(r, ta)]

  3057 	      | NONE => (writeln("### locate_rule:  rewrite "^

  3058 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3059 			 []) end

  3060     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  3061   | locate_rule _ _ _ _ _ _ =

  3062     raise error ("locate_rule: doesnt match rev-sets in istate");

  3064 fun next_rule thy eval_rls ro [rs] t =

  3065     let val der = make_deriv thy eval_rls rs ro NONE t;

  3066     in case der of

  3067 (* val (_,r,_)::_ = der;

  3068    *)

  3069 	   (_,r,_)::_ => SOME r

  3070 	 | _ => NONE

  3071     end

  3072   | next_rule _ _ _ _ _ =

  3073     raise error ("next_rule: doesnt match rev-sets in istate");

  3075 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  3076     []:(rule * (term * term list)) list;

  3078 val pat = (term_of o the o (parse thy)) "?r/?s";

  3079 val pre1 = (term_of o the o (parse thy)) "?r is_expanded";

  3080 val pre2 = (term_of o the o (parse thy)) "?s is_expanded";

  3081 val prepat = [([pre1, pre2], pat)];

  3083 in

  3086 val cancel =

  3087     Rrls {id = "cancel", prepat=prepat,

  3088 	  rew_ord=("ord_make_polynomial",

  3089 		   ord_make_polynomial false Rational.thy),

  3090 	  erls = rational_erls,

  3091 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),

  3092 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),

  3093 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),

  3094 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3095 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3096 		     normal_form = cancel_ thy,

  3097 		     locate_rule = locate_rule thy Atools_erls ro,

  3098 		     next_rule   = next_rule thy Atools_erls ro,

  3099 		     attach_form = attach_form}}

  3100 end;(*local*)

  3104 local(*.ad [2] 'common_nominator_p'

  3105 ---------------------------------

  3106 FIXME Beschreibung .*)

  3109 val {rules=rules,rew_ord=(_,ro),...} =

  3110     rep_rls (assoc_rls "make_polynomial");

  3111 (*WN060829 ... make_deriv does not terminate with 1st expl above,

  3112            see rational.sml --- investigate rulesets for cancel_p ---*)

  3113 val {rules, rew_ord=(_,ro),...} =

  3114     rep_rls (assoc_rls "rev_rew_p");

  3115 val thy = Rational.thy;

  3118 (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option

  3119   as defined above*)

  3121 (*.init_state = fn : term -> istate

  3122 initialzies the state of the interactive interpreter. The state is:

  3124 type rrlsstate =      (*state for reverse rewriting*)

  3125      (term *          (*the current formula*)

  3126       term *          (*the final term*)

  3127       rule list       (*'reverse rule list' (#)*)

  3128 	    list *    (*may be serveral, eg. in norm_rational*)

  3129       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3130        (term *        (*... rewrite with ...*)

  3131 	term list))   (*... assumptions*)

  3132 	  list);      (*derivation from given term to normalform

  3133 		       in reverse order with sym_thm;

  3134                        (#) could be extracted from here by (map #1)*).*)

  3135 fun init_state thy eval_rls ro t =

  3136     let val SOME (t',_) = common_nominator_p_ thy t;

  3137         val SOME (t'',asm) = add_fraction_p_ thy t;

  3138         val der = reverse_deriv thy eval_rls rules ro NONE t';

  3139         val der = der @ [(Thm ("real_mult_div_cancel2",

  3140 			       num_str real_mult_div_cancel2),

  3141 			  (t'',asm))]

  3142         val rs = (distinct_Thm o (map #1)) der;

  3143 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3144 				      "sym_real_mult_0",

  3145 				      "sym_real_mult_1"]) rs;

  3146     in (t,t'',[rs(*here only _ONE_*)],der) end;

  3148 (* use"knowledge/Rational.ML";

  3149    *)

  3151 (*.locate_rule = fn : rule list -> term -> rule

  3152 		      -> (rule * (term * term list) option) list.

  3153   checks a rule R for being a cancel-rule, and if it is,

  3154   then return the list of rules (+ the terms they are rewriting to)

  3155   which need to be applied before R should be applied.

  3156   precondition: the rule is applicable to the argument-term.

  3157 arguments:

  3158   rule list: the reverse rule list

  3159   -> term  : ... to which the rule shall be applied

  3160   -> rule  : ... to be applied to term

  3161 value:

  3162   -> (rule           : a rule rewriting to ...

  3163       * (term        : ... the resulting term ...

  3164          * term list): ... with the assumptions ( //#0).

  3165       ) list         : there may be several such rules;

  3166 		       the list is empty, if the rule has nothing to do

  3167 		       with cancelation.*)

  3168 (* val () = ();

  3169    *)

  3170 fun locate_rule thy eval_rls ro [rs] t r =

  3171     if (id_of_thm r) mem (map (id_of_thm)) rs

  3172     then let val ropt =

  3173 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3174 	 in case ropt of

  3175 		SOME ta => [(r, ta)]

  3176 	      | NONE => (writeln("### locate_rule:  rewrite "^

  3177 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3178 			 []) end

  3179     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  3180   | locate_rule _ _ _ _ _ _ =

  3181     raise error ("locate_rule: doesnt match rev-sets in istate");

  3183 (*.next_rule = fn : rule list -> term -> rule option

  3184   for a given term return the next rules to be done for cancelling.

  3185 arguments:

  3186   rule list     : the reverse rule list

  3187   term          : the term for which ...

  3188 value:

  3189   -> rule option: ... this rule is appropriate for cancellation;

  3190 		  there may be no such rule (if the term is canceled already.*)

  3191 (* val thy = Rational.thy;

  3192    val Rrls {rew_ord=(_,ro),...} = cancel;

  3193    val ([rs],t) = (rss,f);

  3194    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3196    val (thy, [rs]) = (Rational.thy, revsets);

  3197    val Rrls {rew_ord=(_,ro),...} = cancel;

  3198    nex [rs] t;

  3199    *)

  3200 fun next_rule thy eval_rls ro [rs] t =

  3201     let val der = make_deriv thy eval_rls rs ro NONE t;

  3202     in case der of

  3203 (* val (_,r,_)::_ = der;

  3204    *)

  3205 	   (_,r,_)::_ => SOME r

  3206 	 | _ => NONE

  3207     end

  3208   | next_rule _ _ _ _ _ =

  3209     raise error ("next_rule: doesnt match rev-sets in istate");

  3211 (*.val attach_form = f : rule list -> term -> term

  3212 			 -> (rule * (term * term list)) list

  3213   checks an input term TI, if it may belong to a current cancellation, by

  3214   trying to derive it from the given term TG.

  3215 arguments:

  3216   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3217   -> term: TI, the next one input by the user

  3218 value:

  3219   -> (rule           : the rule to be applied in order to reach TI

  3220       * (term        : ... obtained by applying the rule ...

  3221          * term list): ... and the respective assumptions.

  3222       ) list         : there may be several such rules;

  3223                        the list is empty, if the users term does not belong

  3224 		       to a cancellation of the term last agreed upon.*)

  3225 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  3226     []:(rule * (term * term list)) list;

  3228 val pat0 = (term_of o the o (parse thy)) "?r/?s+?u/?v";

  3229 val pat1 = (term_of o the o (parse thy)) "?r/?s+?u   ";

  3230 val pat2 = (term_of o the o (parse thy)) "?r   +?u/?v";

  3231 val prepat = [([HOLogic.true_const], pat0),

  3232 	      ([HOLogic.true_const], pat1),

  3233 	      ([HOLogic.true_const], pat2)];

  3235 in

  3237 (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt;

  3238   besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht

  3239   dh. wie common_nominator_p_, aber auf 1 Bruchstrich*)

  3240 val common_nominator_p =

  3241     Rrls {id = "common_nominator_p", prepat=prepat,

  3242 	  rew_ord=("ord_make_polynomial",

  3243 		   ord_make_polynomial false Rational.thy),

  3244 	  erls = rational_erls,

  3245 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),

  3246 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),

  3247 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),

  3248 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3249 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3250 		     normal_form = add_fraction_p_ thy,(*FIXME.WN0211*)

  3251 		     locate_rule = locate_rule thy Atools_erls ro,

  3252 		     next_rule   = next_rule thy Atools_erls ro,

  3253 		     attach_form = attach_form}}

  3254 end;(*local*)

  3257 local(*.ad [2] 'common_nominator'

  3258 ---------------------------------

  3259 FIXME Beschreibung .*)

  3262 val {rules=rules,rew_ord=(_,ro),...} =

  3263     rep_rls (assoc_rls "make_polynomial");

  3264 val thy = Rational.thy;

  3267 (*.common_nominator_ = fn : theory -> term -> (term * term list) option

  3268   as defined above*)

  3270 (*.init_state = fn : term -> istate

  3271 initialzies the state of the interactive interpreter. The state is:

  3273 type rrlsstate =      (*state for reverse rewriting*)

  3274      (term *          (*the current formula*)

  3275       term *          (*the final term*)

  3276       rule list       (*'reverse rule list' (#)*)

  3277 	    list *    (*may be serveral, eg. in norm_rational*)

  3278       (rule *         (*Thm (+ Thm generated from Calc) resulting in ...*)

  3279        (term *        (*... rewrite with ...*)

  3280 	term list))   (*... assumptions*)

  3281 	  list);      (*derivation from given term to normalform

  3282 		       in reverse order with sym_thm;

  3283                        (#) could be extracted from here by (map #1)*).*)

  3284 fun init_state thy eval_rls ro t =

  3285     let val SOME (t',_) = common_nominator_ thy t;

  3286         val SOME (t'',asm) = add_fraction_ thy t;

  3287         val der = reverse_deriv thy eval_rls rules ro NONE t';

  3288         val der = der @ [(Thm ("real_mult_div_cancel2",

  3289 			       num_str real_mult_div_cancel2),

  3290 			  (t'',asm))]

  3291         val rs = (distinct_Thm o (map #1)) der;

  3292 	val rs = filter_out (eq_Thms ["sym_real_add_zero_left",

  3293 				      "sym_real_mult_0",

  3294 				      "sym_real_mult_1"]) rs;

  3295     in (t,t'',[rs(*here only _ONE_*)],der) end;

  3297 (* use"knowledge/Rational.ML";

  3298    *)

  3300 (*.locate_rule = fn : rule list -> term -> rule

  3301 		      -> (rule * (term * term list) option) list.

  3302   checks a rule R for being a cancel-rule, and if it is,

  3303   then return the list of rules (+ the terms they are rewriting to)

  3304   which need to be applied before R should be applied.

  3305   precondition: the rule is applicable to the argument-term.

  3306 arguments:

  3307   rule list: the reverse rule list

  3308   -> term  : ... to which the rule shall be applied

  3309   -> rule  : ... to be applied to term

  3310 value:

  3311   -> (rule           : a rule rewriting to ...

  3312       * (term        : ... the resulting term ...

  3313          * term list): ... with the assumptions ( //#0).

  3314       ) list         : there may be several such rules;

  3315 		       the list is empty, if the rule has nothing to do

  3316 		       with cancelation.*)

  3317 (* val () = ();

  3318    *)

  3319 fun locate_rule thy eval_rls ro [rs] t r =

  3320     if (id_of_thm r) mem (map (id_of_thm)) rs

  3321     then let val ropt =

  3322 		 rewrite_ thy ro eval_rls true (thm_of_thm r) t;

  3323 	 in case ropt of

  3324 		SOME ta => [(r, ta)]

  3325 	      | NONE => (writeln("### locate_rule:  rewrite "^

  3326 				 (id_of_thm r)^" "^(term2str t)^" = NONE");

  3327 			 []) end

  3328     else (writeln("### locate_rule:  "^(id_of_thm r)^" not mem rrls");[])

  3329   | locate_rule _ _ _ _ _ _ =

  3330     raise error ("locate_rule: doesnt match rev-sets in istate");

  3332 (*.next_rule = fn : rule list -> term -> rule option

  3333   for a given term return the next rules to be done for cancelling.

  3334 arguments:

  3335   rule list     : the reverse rule list

  3336   term          : the term for which ...

  3337 value:

  3338   -> rule option: ... this rule is appropriate for cancellation;

  3339 		  there may be no such rule (if the term is canceled already.*)

  3340 (* val thy = Rational.thy;

  3341    val Rrls {rew_ord=(_,ro),...} = cancel;

  3342    val ([rs],t) = (rss,f);

  3343    next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*)

  3345    val (thy, [rs]) = (Rational.thy, revsets);

  3346    val Rrls {rew_ord=(_,ro),...} = cancel_p;

  3347    nex [rs] t;

  3348    *)

  3349 fun next_rule thy eval_rls ro [rs] t =

  3350     let val der = make_deriv thy eval_rls rs ro NONE t;

  3351     in case der of

  3352 (* val (_,r,_)::_ = der;

  3353    *)

  3354 	   (_,r,_)::_ => SOME r

  3355 	 | _ => NONE

  3356     end

  3357   | next_rule _ _ _ _ _ =

  3358     raise error ("next_rule: doesnt match rev-sets in istate");

  3360 (*.val attach_form = f : rule list -> term -> term

  3361 			 -> (rule * (term * term list)) list

  3362   checks an input term TI, if it may belong to a current cancellation, by

  3363   trying to derive it from the given term TG.

  3364 arguments:

  3365   term   : TG, the last one in the cancellation agreed upon by user + math-eng

  3366   -> term: TI, the next one input by the user

  3367 value:

  3368   -> (rule           : the rule to be applied in order to reach TI

  3369       * (term        : ... obtained by applying the rule ...

  3370          * term list): ... and the respective assumptions.

  3371       ) list         : there may be several such rules;

  3372                        the list is empty, if the users term does not belong

  3373 		       to a cancellation of the term last agreed upon.*)

  3374 fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*)

  3375     []:(rule * (term * term list)) list;

  3377 val pat0 =  (term_of o the o (parse thy)) "?r/?s+?u/?v";

  3378 val pat01 = (term_of o the o (parse thy)) "?r/?s-?u/?v";

  3379 val pat1 =  (term_of o the o (parse thy)) "?r/?s+?u   ";

  3380 val pat11 = (term_of o the o (parse thy)) "?r/?s-?u   ";

  3381 val pat2 =  (term_of o the o (parse thy)) "?r   +?u/?v";

  3382 val pat21 = (term_of o the o (parse thy)) "?r   -?u/?v";

  3383 val prepat = [([HOLogic.true_const], pat0),

  3384 	      ([HOLogic.true_const], pat01),

  3385 	      ([HOLogic.true_const], pat1),

  3386 	      ([HOLogic.true_const], pat11),

  3387 	      ([HOLogic.true_const], pat2),

  3388 	      ([HOLogic.true_const], pat21)];

  3391 in

  3393 val common_nominator =

  3394     Rrls {id = "common_nominator", prepat=prepat,

  3395 	  rew_ord=("ord_make_polynomial",

  3396 		   ord_make_polynomial false Rational.thy),

  3397 	  erls = rational_erls,

  3398 	  calc = [("PLUS"    ,("op +"        ,eval_binop "#add_")),

  3399 		  ("TIMES"   ,("op *"        ,eval_binop "#mult_")),

  3400 		  ("DIVIDE" ,("HOL.divide"  ,eval_cancel "#divide_")),

  3401 		  ("POWER"  ,("Atools.pow"  ,eval_binop "#power_"))],

  3402 	  (*asm_thm=[("real_mult_div_cancel2","")],*)

  3403 	  scr=Rfuns {init_state  = init_state thy Atools_erls ro,

  3404 		     normal_form = add_fraction_ (*NOT common_nominator_*) thy,

  3405 		     locate_rule = locate_rule thy Atools_erls ro,

  3406 		     next_rule   = next_rule thy Atools_erls ro,

  3407 		     attach_form = attach_form}}

  3409 end;(*local*)

  3412 (*##*)

  3413 end;(*struct*)

  3415 open RationalI;

  3416 (*##*)

  3418 (*.the expression contains + - * ^ / only ?.*)

  3419 fun is_ratpolyexp (Free _) = true

  3420   | is_ratpolyexp (Const ("op +",_) $ Free _ $ Free _) = true

  3421   | is_ratpolyexp (Const ("op -",_) $ Free _ $ Free _) = true

  3422   | is_ratpolyexp (Const ("op *",_) $ Free _ $ Free _) = true

  3423   | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true

  3424   | is_ratpolyexp (Const ("HOL.divide",_) $ Free _ $ Free _) = true

  3425   | is_ratpolyexp (Const ("op +",_) $ t1 $ t2) =

  3426                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3427   | is_ratpolyexp (Const ("op -",_) $ t1 $ t2) =

  3428                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3429   | is_ratpolyexp (Const ("op *",_) $ t1 $ t2) =

  3430                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3431   | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) =

  3432                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3433   | is_ratpolyexp (Const ("HOL.divide",_) $ t1 $ t2) =

  3434                ((is_ratpolyexp t1) andalso (is_ratpolyexp t2))

  3435   | is_ratpolyexp _ = false;

  3437 (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*)

  3438 fun eval_is_ratpolyexp (thmid:string) _

  3439 		       (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy =

  3440     if is_ratpolyexp arg

  3441     then SOME (mk_thmid thmid ""

  3442 			((Syntax.string_of_term (thy2ctxt thy)) arg) "",

  3443 	       Trueprop $ (mk_equality (t, HOLogic.true_const)))

  3444     else SOME (mk_thmid thmid ""

  3445 			((Syntax.string_of_term (thy2ctxt thy)) arg) "",

  3446 	       Trueprop $ (mk_equality (t, HOLogic.false_const)))

  3447   | eval_is_ratpolyexp _ _ _ _ = NONE;

  3451 (*-------------------18.3.03 --> struct <-----------vvv--*)

  3452 val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*)

  3454 (*.discard binary minus, shift unary minus into -1*;

  3455    unary minus before numerals are put into the numeral by parsing;

  3456    contains absolute minimum of thms for context in norm_Rational .*)

  3457 val discard_minus = prep_rls(

  3458   Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3459       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)

  3460       rules = [Thm ("real_diff_minus", num_str real_diff_minus),

  3461 	       (*"a - b = a + -1 * b"*)

  3462 	       Thm ("sym_real_mult_minus1",num_str (real_mult_minus1 RS sym))

  3463 	       (*- ?z = "-1 * ?z"*)

  3464 	       ],

  3465       scr = Script ((term_of o the o (parse thy))

  3466       "empty_script")

  3467       }):rls;

  3468 (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*)

  3469 val powers_erls = prep_rls(

  3470   Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3471       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)

  3472       rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"),

  3473 	       Calc ("Atools.is'_even",eval_is_even "#is_even_"),

  3474 	       Calc ("op <",eval_equ "#less_"),

  3475 	       Thm ("not_false", not_false),

  3476 	       Thm ("not_true", not_true),

  3477 	       Calc ("op +",eval_binop "#add_")

  3478 	       ],

  3479       scr = Script ((term_of o the o (parse thy))

  3480       "empty_script")

  3481       }:rls);

  3482 (*.all powers over + distributed; atoms over * collected, other distributed

  3483    contains absolute minimum of thms for context in norm_Rational .*)

  3484 val powers = prep_rls(

  3485   Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3486       erls = powers_erls, srls = Erls, calc = [], (*asm_thm = [],*)

  3487       rules = [Thm ("realpow_multI", num_str realpow_multI),

  3488 	       (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*)

  3489 	       Thm ("realpow_pow",num_str realpow_pow),

  3490 	       (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*)

  3491 	       Thm ("realpow_oneI",num_str realpow_oneI),

  3492 	       (*"r ^^^ 1 = r"*)

  3493 	       Thm ("realpow_minus_even",num_str realpow_minus_even),

  3494 	       (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*)

  3495 	       Thm ("realpow_minus_odd",num_str realpow_minus_odd),

  3496 	       (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*)

  3498 	       (*----- collect atoms over * -----*)

  3499 	       Thm ("realpow_two_atom",num_str realpow_two_atom),

  3500 	       (*"r is_atom ==> r * r = r ^^^ 2"*)

  3501 	       Thm ("realpow_plus_1",num_str realpow_plus_1),

  3502 	       (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*)

  3503 	       Thm ("realpow_addI_atom",num_str realpow_addI_atom),

  3504 	       (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*)

  3506 	       (*----- distribute none-atoms -----*)

  3507 	       Thm ("realpow_def_atom",num_str realpow_def_atom),

  3508 	       (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*)

  3509 	       Thm ("realpow_eq_oneI",num_str realpow_eq_oneI),

  3510 	       (*"1 ^^^ n = 1"*)

  3511 	       Calc ("op +",eval_binop "#add_")

  3512 	       ],

  3513       scr = Script ((term_of o the o (parse thy))

  3514       "empty_script")

  3515       }:rls);

  3516 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3517 val rat_mult_divide = prep_rls(

  3518   Rls {id = "rat_mult_divide", preconds = [],

  3519        rew_ord = ("dummy_ord",dummy_ord),

  3520       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)

  3521       rules = [Thm ("rat_mult",num_str rat_mult),

  3522 	       (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3523 	       Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),

  3524 	       (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3525 	       otherwise inv.to a / b / c = ...*)

  3526 	       Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),

  3527 	       (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much

  3528 		     and does not commute a / b * c ^^^ 2 !*)

  3530 	       Thm ("real_divide_divide1_eq", real_divide_divide1_eq),

  3531 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3532 	       Thm ("real_divide_divide2_eq", real_divide_divide2_eq),

  3533 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3534 	       Calc ("HOL.divide"  ,eval_cancel "#divide_")

  3535 	       ],

  3536       scr = Script ((term_of o the o (parse thy)) "empty_script")

  3537       }:rls);

  3538 (*.contains absolute minimum of thms for context in norm_Rational.*)

  3539 val reduce_0_1_2 = prep_rls(

  3540   Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord),

  3541       erls = e_rls,srls = Erls,calc = [],(*asm_thm = [],*)

  3542       rules = [(*Thm ("real_divide_1",num_str real_divide_1),

  3543 		 "?x / 1 = ?x" unnecess.for normalform*)

  3544 	       Thm ("real_mult_1",num_str real_mult_1),

  3545 	       (*"1 * z = z"*)

  3546 	       (*Thm ("real_mult_minus1",num_str real_mult_minus1),

  3547 	       "-1 * z = - z"*)

  3548 	       (*Thm ("real_minus_mult_cancel",num_str real_minus_mult_cancel),

  3549 	       "- ?x * - ?y = ?x * ?y"*)

  3551 	       Thm ("real_mult_0",num_str real_mult_0),

  3552 	       (*"0 * z = 0"*)

  3553 	       Thm ("real_add_zero_left",num_str real_add_zero_left),

  3554 	       (*"0 + z = z"*)

  3555 	       (*Thm ("real_add_minus",num_str real_add_minus),

  3556 	       "?z + - ?z = 0"*)

  3558 	       Thm ("sym_real_mult_2",num_str (real_mult_2 RS sym)),

  3559 	       (*"z1 + z1 = 2 * z1"*)

  3560 	       Thm ("real_mult_2_assoc",num_str real_mult_2_assoc),

  3561 	       (*"z1 + (z1 + k) = 2 * z1 + k"*)

  3563 	       Thm ("real_0_divide",num_str real_0_divide)

  3564 	       (*"0 / ?x = 0"*)

  3565 	       ], scr = EmptyScr}:rls);

  3567 (*erls for calculate_Rational;

  3568   make local with FIXX@ME result:term *term list WN0609???SKMG*)

  3569 val norm_rat_erls = prep_rls(

  3570   Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3571       erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*)

  3572       rules = [Calc ("Atools.is'_const",eval_const "#is_const_")

  3573 	       ],

  3574       scr = Script ((term_of o the o (parse thy))

  3575       "empty_script")

  3576       }:rls);

  3577 (*.consists of rls containing the absolute minimum of thms.*)

  3578 (*040209: this version has been used by RL for his equations,

  3579 which is now replaced by MGs version below

  3580 vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)

  3581 val norm_Rational = prep_rls(

  3582   Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord),

  3583       erls = norm_rat_erls, srls = Erls, calc = [], (*asm_thm = [],*)

  3584       rules = [(*sequence given by operator precedence*)

  3585 	       Rls_ discard_minus,

  3586 	       Rls_ powers,

  3587 	       Rls_ rat_mult_divide,

  3588 	       Rls_ expand,

  3589 	       Rls_ reduce_0_1_2,

  3590 	       (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*)

  3591 	       Rls_ order_add_mult,

  3592 	       Rls_ collect_numerals,

  3593 	       Rls_ add_fractions_p,

  3594 	       Rls_ cancel_p

  3595 	       ],

  3596       scr = Script ((term_of o the o (parse thy))

  3597       "empty_script")

  3598       }:rls);

  3599 val norm_Rational_parenthesized = prep_rls(

  3600   Seq {id = "norm_Rational_parenthesized", preconds = []:term list,

  3601        rew_ord = ("dummy_ord", dummy_ord),

  3602       erls = Atools_erls, srls = Erls,

  3603       calc = [], (*asm_thm = [],*)

  3604       rules = [Rls_  norm_Rational, (*from RL -- not the latest one*)

  3605 	       Rls_ discard_parentheses

  3606 	       ],

  3607       scr = EmptyScr

  3608       }:rls);

  3611 (*-------------------18.3.03 --> struct <-----------^^^--*)

  3615 theory' := overwritel (!theory', [("Rational.thy",Rational.thy)]);

  3618 (*WN030318???SK: simplifies all but cancel and common_nominator*)

  3619 val simplify_rational =

  3620     merge_rls "simplify_rational" expand_binoms

  3621     (append_rls "divide" calculate_Rational

  3622 		[Thm ("real_divide_1",num_str real_divide_1),

  3623 		 (*"?x / 1 = ?x"*)

  3624 		 Thm ("rat_mult",num_str rat_mult),

  3625 		 (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3626 		 Thm ("real_times_divide1_eq",num_str real_times_divide1_eq),

  3627 		 (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2],

  3628 		 otherwise inv.to a / b / c = ...*)

  3629 		 Thm ("real_times_divide2_eq",num_str real_times_divide2_eq),

  3630 		 (*"?a / ?b * ?c = ?a * ?c / ?b"*)

  3631 		 Thm ("add_minus",num_str add_minus),

  3632 		 (*"?a + ?b - ?b = ?a"*)

  3633 		 Thm ("add_minus1",num_str add_minus1),

  3634 		 (*"?a - ?b + ?b = ?a"*)

  3635 		 Thm ("real_divide_minus1",num_str real_divide_minus1)

  3636 		 (*"?x / -1 = - ?x"*)

  3637 (*

  3638 ,

  3639 		 Thm ("",num_str )

  3640 *)

  3641 		 ]);

  3643 (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*)

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

  3646 (*                  Simplifier für beliebige Buchterme                *)

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

  3648 (*----------------------- norm_Rational_mg ---------------------------*)

  3649 (*. description of the simplifier see MG-DA.p.56ff .*)

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

  3651 val common_nominator_p_rls = prep_rls(

  3652   Rls {id = "common_nominator_p_rls", preconds = [],

  3653        rew_ord = ("dummy_ord",dummy_ord),

  3654 	 erls = e_rls, srls = Erls, calc = [],

  3655 	 rules =

  3656 	 [Rls_ common_nominator_p

  3657 	  (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ?

  3658 	    FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*)

  3659 	  ],

  3660 	 scr = EmptyScr});

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

  3662 val cancel_p_rls = prep_rls(

  3663   Rls {id = "cancel_p_rls", preconds = [],

  3664        rew_ord = ("dummy_ord",dummy_ord),

  3665 	 erls = e_rls, srls = Erls, calc = [],

  3666 	 rules =

  3667 	 [Rls_ cancel_p

  3668 	  (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*)

  3669 	  ],

  3670 	 scr = EmptyScr});

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

  3672 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3673     used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*)

  3674 val rat_mult_poly = prep_rls(

  3675   Rls {id = "rat_mult_poly", preconds = [],

  3676        rew_ord = ("dummy_ord",dummy_ord),

  3677 	 erls =  append_rls "e_rls-is_polyexp" e_rls

  3678 	         [Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3679 	 srls = Erls, calc = [],

  3680 	 rules =

  3681 	 [Thm ("rat_mult_poly_l",num_str rat_mult_poly_l),

  3682 	  (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3683 	  Thm ("rat_mult_poly_r",num_str rat_mult_poly_r)

  3684 	  (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3685 	  ],

  3686 	 scr = EmptyScr});

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

  3688 (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions;

  3689     used in looping part norm_Rational_rls, see example DA-M02-main.p.60

  3690     .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls,

  3691     I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028

  3692     ... WN0609???MG.*)

  3693 val rat_mult_div_pow = prep_rls(

  3694   Rls {id = "rat_mult_div_pow", preconds = [],

  3695        rew_ord = ("dummy_ord",dummy_ord),

  3696        erls = e_rls,

  3697        (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls

  3698 			[Calc ("Poly.is'_polyexp", eval_is_polyexp "")],

  3699          with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get

  3700 	 error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc.

  3701          thus we decided to go on with this flaw*)

  3702 		 srls = Erls, calc = [],

  3703       rules = [Thm ("rat_mult",num_str rat_mult),

  3704 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

  3705 	       Thm ("rat_mult_poly_l",num_str rat_mult_poly_l),

  3706 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

  3707 	       Thm ("rat_mult_poly_r",num_str rat_mult_poly_r),

  3708 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

  3710 	       Thm ("real_divide_divide1_mg", real_divide_divide1_mg),

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

  3712 	       Thm ("real_divide_divide1_eq", real_divide_divide1_eq),

  3713 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

  3714 	       Thm ("real_divide_divide2_eq", real_divide_divide2_eq),

  3715 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

  3716 	       Calc ("HOL.divide"  ,eval_cancel "#divide_"),

  3718 	       Thm ("rat_power", num_str rat_power)

  3719 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

  3720 	       ],

  3721       scr = Script ((term_of o the o (parse thy)) "empty_script")

  3722       }:rls);

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

  3724 val rat_reduce_1 = prep_rls(

  3725   Rls {id = "rat_reduce_1", preconds = [],

  3726        rew_ord = ("dummy_ord",dummy_ord),

  3727        erls = e_rls, srls = Erls, calc = [],

  3728        rules = [Thm ("real_divide_1",num_str real_divide_1),

  3729 		(*"?x / 1 = ?x"*)

  3730 		Thm ("real_mult_1",num_str real_mult_1)

  3731 		(*"1 * z = z"*)

  3732 		],

  3733        scr = Script ((term_of o the o (parse thy)) "empty_script")

  3734        }:rls);

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

  3736 (*. looping part of norm_Rational(*_mg*) .*)

  3737 val norm_Rational_rls = prep_rls(

  3738    Rls {id = "norm_Rational_rls", preconds = [],

  3739        rew_ord = ("dummy_ord",dummy_ord),

  3740        erls = norm_rat_erls, srls = Erls, calc = [],

  3741        rules = [Rls_ common_nominator_p_rls,

  3742 		Rls_ rat_mult_div_pow,

  3743 		Rls_ make_rat_poly_with_parentheses,

  3744 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)

  3745 		Rls_ rat_reduce_1

  3746 		],

  3747        scr = Script ((term_of o the o (parse thy)) "empty_script")

  3748        }:rls);

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

  3750 (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG)

  3751  just be renaming:*)

  3752 val norm_Rational(*_mg*) = prep_rls(

  3753    Seq {id = "norm_Rational"(*_mg*), preconds = [],

  3754        rew_ord = ("dummy_ord",dummy_ord),

  3755        erls = norm_rat_erls, srls = Erls, calc = [],

  3756        rules = [Rls_ discard_minus_,

  3757 		Rls_ rat_mult_poly,(* removes double fractions like a/b/c    *)

  3758 		Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)

  3759 		Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)

  3760 		Rls_ norm_Rational_rls,   (* the main rls, looping (#)       *)

  3761 		Rls_ discard_parentheses_ (* mult only                       *)

  3762 		],

  3763        scr = Script ((term_of o the o (parse thy)) "empty_script")

  3764        }:rls);

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

  3768 ruleset' := overwritelthy thy (!ruleset',

  3769   [("calculate_Rational", calculate_Rational),

  3770    ("calc_rat_erls",calc_rat_erls),

  3771    ("rational_erls", rational_erls),

  3772    ("cancel_p", cancel_p),

  3773    ("cancel", cancel),

  3774    ("common_nominator_p", common_nominator_p),

  3775    ("common_nominator_p_rls", common_nominator_p_rls),

  3776    ("common_nominator"  , common_nominator),

  3777    ("discard_minus", discard_minus),

  3778    ("powers_erls", powers_erls),

  3779    ("powers", powers),

  3780    ("rat_mult_divide", rat_mult_divide),

  3781    ("reduce_0_1_2", reduce_0_1_2),

  3782    ("rat_reduce_1", rat_reduce_1),

  3783    ("norm_rat_erls", norm_rat_erls),

  3784    ("norm_Rational", norm_Rational),

  3785    ("norm_Rational_rls", norm_Rational_rls),

  3786    ("norm_Rational_parenthesized", norm_Rational_parenthesized),

  3787    ("rat_mult_poly", rat_mult_poly),

  3788    ("rat_mult_div_pow", rat_mult_div_pow),

  3789    ("cancel_p_rls", cancel_p_rls)

  3790    ]);

  3792 calclist':= overwritel (!calclist',

  3793    [("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))

  3794     ]);

  3796 (** problems **)

  3798 store_pbt

  3799  (prep_pbt Rational.thy "pbl_simp_rat" [] e_pblID

  3800  (["rational","simplification"],

  3801   [("#Given" ,["term t_"]),

  3802    ("#Where" ,["t_ is_ratpolyexp"]),

  3803    ("#Find"  ,["normalform n_"])

  3804   ],

  3805   append_rls "e_rls" e_rls [(*for preds in where_*)],

  3806   SOME "Simplify t_",

  3807   [["simplification","of_rationals"]]));

  3809 (** methods **)

  3811 (*WN061025 this methods script is copied from (auto-generated) script

  3812   of norm_Rational in order to ease repair on inform*)

  3813 store_met

  3814     (prep_met Rational.thy "met_simp_rat" [] e_metID

  3815 	      (["simplification","of_rationals"],

  3816 	       [("#Given" ,["term t_"]),

  3817 		("#Where" ,["t_ is_ratpolyexp"]),

  3818 		("#Find"  ,["normalform n_"])

  3819 		],

  3820 	       {rew_ord'="tless_true",

  3821 		rls' = e_rls,

  3822 		calc = [], srls = e_rls,

  3823 		prls = append_rls "simplification_of_rationals_prls" e_rls

  3824 				[(*for preds in where_*)

  3825 				 Calc ("Rational.is'_ratpolyexp",

  3826 				       eval_is_ratpolyexp "")],

  3827 		crls = e_rls, nrls = norm_Rational_rls},

  3828 "Script SimplifyScript (t_::real) =                              " ^

  3829 "  ((Try (Rewrite_Set discard_minus_ False) @@                   " ^

  3830 "    Try (Rewrite_Set rat_mult_poly False) @@                    " ^

  3831 "    Try (Rewrite_Set make_rat_poly_with_parentheses False) @@   " ^

  3832 "    Try (Rewrite_Set cancel_p_rls False) @@                     " ^

  3833 "    (Repeat                                                     " ^

  3834 "     ((Try (Rewrite_Set common_nominator_p_rls False) @@        " ^

  3835 "       Try (Rewrite_Set rat_mult_div_pow False) @@              " ^

  3836 "       Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^

  3837 "       Try (Rewrite_Set cancel_p_rls False) @@                  " ^

  3838 "       Try (Rewrite_Set rat_reduce_1 False)))) @@               " ^

  3839 "    Try (Rewrite_Set discard_parentheses_ False))               " ^

  3840 "    t_)"

  3841 	       ));

  3842 *}

  3844 end

author	Walther Neuper <neuper@ist.tugraz.at>
	Thu, 26 Aug 2010 18:15:30 +0200
branch	isac-update-Isa09-2
changeset 37950	525a28152a67
parent 37947	22235e4dbe5f
child 37953	369b3012f6f6
permissions	-rw-r--r--